Example 3: Kinetic energy dissipation in a pipe

In this example, we show how to minimize the shape functional

\[\mathcal{J}(\Omega) = \int_\Omega \nu \nabla \mathbf{u} : \nabla \mathbf{u} \,\mathrm{d}\mathbf{x}\,,\]

where \(\mathbf{u}:\mathbb{R}^d\to\mathbb{R^d}\,, d = 2,3,\) is the velocity of an incompressible fluid and \(\nu\) is the fluid viscosity. The fluid velocity \(\mathbf{u}\) and the fluid pressure \(p:\mathbb{R}^d\to\mathbb{R}\) satisfy the incompressible Navier-Stokes equations

\begin{align*} -\nu \Delta \mathbf{u} + \mathbf{u} \nabla \mathbf{u} + \nabla p &= 0 & \text{in } \Omega\,,\\ \operatorname{div} \mathbf{u} &= 0& \text{in } \Omega\,,\\ \mathbf{u} &= \mathbf{g} &\text{on } \partial\Omega\setminus \Gamma\,,\\ p\mathbf{n} - \nu \nabla u\cdot \mathbf{n} & = 0 &\text{on } \Gamma\,. \end{align*}

Here, \(\mathbf{g}\) is given by a Poiseuille flow at the inlet and is zero on the walls of the pipe. The letter \(\Gamma\) denotes the outlet.

In addition to the PDE-contstraint, we enforce a volume constraint: the volume of the optimized domain should be equal to the volume of the initial domain.

Initial domain

pipe2D_initial pipe3D_initial

Initial pipe design in 2D (with mesh) and 3D, and magnitude of fluid velocity.

For the 2D-example, the geometry of the initial domain is described in the following gmsh script (this has been tested with gmsh v4.1.2).

 1h = 1.; //width
 2H = 5.; //height
 3
 4//lower points
 5Point(1) = { 0, 0, 0, 1.};
 6Point(2) = { 2, 0, 0, 1.}; //first spline point
 7Point(3) = { 4, 0, 0, 1.};
 8Point(4) = { 8, 6, 0, 1.};
 9Point(5) = {10, 6, 0, 1.};
10Point(6) = {12, 6, 0, 1.}; //last spline point
11Point(7) = {15, 6, 0, 1.};
12//upper points
13Point( 8) = {15, 7, 0, 1.};
14Point( 9) = {12, 7, 0, 1.}; //first spline point
15Point(10) = {10, 7, 0, 1.};
16Point(11) = { 8, 7, 0, 1.};
17Point(12) = { 4, 1, 0, 1.};
18Point(13) = { 2, 1, 0, 1.}; //last spline point
19Point(14) = { 0, 1, 0, 1.};
20
21//edges
22Line(1) = {1, 2};
23BSpline(2) = {2, 3, 4, 5, 6};
24Line(3) = {6, 7};
25Line(4) = {7, 8};
26Line(5) = {8, 9};
27BSpline(6) = {9, 10, 11, 12, 13};
28Line(7) = {13, 14};
29Line(8) = {14,  1};
30
31//boundary and physical curves
32Curve Loop(9) = {1, 2, 3, 4, 5, 6, 7, 8};
33Physical Curve("Inflow", 10) = {8};
34Physical Curve("Outflow", 11) = {4};
35Physical Curve("WallFixed", 12) = {1, 3, 5, 7};
36Physical Curve("WallFree", 13) = {2, 6};
37
38
39//domain and physical surface
40Plane Surface(1) = {9};
41Physical Surface("Pipe", 2) = {1};

The mesh can be generated typing gmsh -2 -clscale 0.1 -format msh2 -o pipe.msh pipe2d.geo in the terminal.

For the 3D-example, the geometry of the initial domain is described in the following gmsh script.

 1// Gmsh project created on Tue Jan 22 11:40:52 2019
 2SetFactory("OpenCASCADE");
 3Circle(1) = {0, 0, 0, 0.5, 0, 2*Pi};
 4Line Loop(2) = {1};
 5Plane Surface(1) = {2};
 6Point( 5) = {0, 0.0,  0.00, 1.0};
 7Point( 6) = {0, 0.0,  2.00, 1.0}; //first spline point
 8Point( 7) = {0, 0.0,  3.00, 1.0};
 9Point( 8) = {0, 0.0,  4.00, 1.0};
10Point( 9) = {0, 0.2,  4.75, 1.0};
11Point(10) = {0, 5.0,  6.00, 1.0};
12Point(11) = {0, 5.0, 10.00, 1.0};
13Point(12) = {0, 5.0, 12.00, 1.0}; //last spline point
14Point(13) = {0, 5.0, 15.00, 1.0};
15
16Bezier(10) = {6, 7, 8, 9, 10, 11, 12};
17Line(15) = {5, 6};
18Line(16) = {12, 13};
19Wire(1) = {15, 10, 16};
20Extrude { Surface{1}; } Using Wire {1}
21Delete{ Surface{1}; }
22
23// When using fancy commands like "extrude", it is not quite obvious
24// what the numbering of the generated lines, surfaces and volumes is
25// going to be. To figure this out, we open the .geo file in the gmsh
26// gui and use the Tools -> Visibility menu to find the numbers for
27// each entity.
28
29Physical Surface("Inflow", 10) = {2};
30Physical Surface("Outflow", 11) = {6};
31Physical Surface("WallFixed", 12) = {3, 5};
32Physical Surface("WallFree", 13) = {4};
33//Physical Surface("Inflow") = {2};
34//Physical Surface("Outflow") = {6};
35//Physical Surface("WallFixed") = {3, 5};
36//Physical Surface("WallFree") = {4};
37Physical Volume("PhysVol") = {1};

The mesh can be generated typing gmsh -3 -clscale 0.2 -format msh2 -o pipe.msh pipe3d.geo in the terminal (this has been tested with gmsh v4.1.2).

Implementing the PDE constraint

We implement the boundary value problem that acts as PDE constraint in a python module named PDEconstraint_pipe.py. In the code, we highlight the lines which characterize the weak formulation of this boundary value problem.

Note

The Dirichlet boundary data \(\mathbf{g}\) depends on the dimension \(d\) (see Lines 36-42)

 1import firedrake as fd
 2from fireshape import PdeConstraint
 3
 4
 5class NavierStokesSolver(PdeConstraint):
 6    """Incompressible Navier-Stokes as PDE constraint."""
 7
 8    def __init__(self, mesh_m, viscosity):
 9        super().__init__()
10        self.mesh_m = mesh_m
11        self.failed_to_solve = False  # when self.solver.solve() fail
12
13        # Setup problem, Taylor-Hood finite elements
14        self.V = fd.VectorFunctionSpace(self.mesh_m, "CG", 2) \
15            * fd.FunctionSpace(self.mesh_m, "CG", 1)
16
17        # Preallocate solution variables for state equation
18        self.solution = fd.Function(self.V, name="State")
19        self.testfunction = fd.TestFunction(self.V)
20
21        # Define viscosity parameter
22        self.viscosity = viscosity
23
24        # Weak form of incompressible Navier-Stokes equations
25        z = self.solution
26        u, p = fd.split(z)
27        test = self.testfunction
28        v, q = fd.split(test)
29        nu = self.viscosity  # shorten notation
30        self.F = nu*fd.inner(fd.grad(u), fd.grad(v))*fd.dx - p*fd.div(v)*fd.dx\
31            + fd.inner(fd.dot(fd.grad(u), u), v)*fd.dx + fd.div(u)*q*fd.dx
32
33        # Dirichlet Boundary conditions
34        X = fd.SpatialCoordinate(self.mesh_m)
35        dim = self.mesh_m.topological_dimension()
36        if dim == 2:
37            uin = 4 * fd.as_vector([(1-X[1])*X[1], 0])
38        elif dim == 3:
39            rsq = X[0]**2+X[1]**2  # squared radius = 0.5**2 = 1/4
40            uin = fd.as_vector([0, 0, 1-4*rsq])
41        else:
42            raise NotImplementedError
43        self.bcs = [fd.DirichletBC(self.V.sub(0), 0., [12, 13]),
44                    fd.DirichletBC(self.V.sub(0), uin, [10])]
45
46        # PDE-solver parameters
47        self.nsp = None
48        self.params = {
49            "snes_max_it": 10, "mat_type": "aij", "pc_type": "lu",
50            "pc_factor_mat_solver_type": "superlu_dist",
51            # "snes_monitor": None, "ksp_monitor": None,
52        }
53
54    def solve(self):
55        super().solve()
56        self.failed_to_solve = False
57        u_old = self.solution.copy(deepcopy=True)
58        try:
59            fd.solve(self.F == 0, self.solution, bcs=self.bcs,
60                     solver_parameters=self.params)
61        except fd.ConvergenceError:
62            self.failed_to_solve = True
63            self.solution.assign(u_old)
64
65
66if __name__ == "__main__":
67    mesh = fd.Mesh("pipe.msh")
68    if mesh.topological_dimension() == 2:  # in 2D
69        viscosity = fd.Constant(1./400.)
70    elif mesh.topological_dimension() == 3:  # in 3D
71        viscosity = fd.Constant(1/10.)  # simpler problem in 3D
72    else:
73        raise NotImplementedError
74    e = NavierStokesSolver(mesh, viscosity)
75    e.solve()
76    print(e.failed_to_solve)
77    out = fd.File("temp_PDEConstrained_u.pvd")
78    out.write(e.solution.split()[0])

Note

The Navier-Stokes solver may fail to converge if too big an optimization step occurs in the optimization process. To address this issue, we use a trust-region algorithm as optimization solver, and we make the functional \(\mathcal{J}\) return NaN whenever the state solver fails. This way, the trust-region method will notice that there is no improvement if the state solver fails and will thus reduce the trust-region radius.

_images/pipe2D_fail.png

If the trust-region radius is too large, the algorithm may try to evaluate the objective functional on a domain that is not feasible. In this case, the PDE-solver fails, the domain is rejected, and the trust-region radius is reduced.

Implementing the shape functional

We implement the shape functional \(\mathcal{J}\) in a python file named objective_pipe.py. In the code, we highlight the lines which characterize \(\mathcal{J}\).

 1import firedrake as fd
 2from fireshape import ShapeObjective
 3from PDEconstraint_pipe import NavierStokesSolver
 4import numpy as np
 5
 6
 7class PipeObjective(ShapeObjective):
 8    """L2 tracking functional for Poisson problem."""
 9
10    def __init__(self, pde_solver: NavierStokesSolver, *args, **kwargs):
11        super().__init__(*args, **kwargs)
12        self.pde_solver = pde_solver
13
14    def value_form(self):
15        """Evaluate misfit functional."""
16        nu = self.pde_solver.viscosity
17
18        if self.pde_solver.failed_to_solve:  # return NaNs if state solve fails
19            return np.nan * fd.dx(self.pde_solver.mesh_m)
20        else:
21            z = self.pde_solver.solution
22            u, p = fd.split(z)
23            return nu * fd.inner(fd.grad(u), fd.grad(u)) * fd.dx

Setting up and solving the problem

We set up the problem in the script main_pipe.py.

To set up the problem, we need to:

  • load the mesh of the initial guess (Line 9),

  • choose the discretization of the control (Line 10, Lagrangian finite elements of degree 1),

  • choose the metric of the control space (Line 11, \(H^1\)-seminorm with homogeneous Dirichlet boundary conditions on fixed boundaries),

  • initialize the PDE contraint on the physical mesh mesh_m (Line 15-21), choosing different viscosity parameters depending on the physical dimension dim

  • specify to save the function \(\mathbf{u}\) after each iteration in the file solution/u.pvd by setting the function cb appropriately (Lines 24-28),

  • initialize the shape functional (Line 31), and the reduced shape functional (Line 32),

  • add a regularization term to improve the mesh quality in the updated domains (Lines 35-36),

  • specify the volume equality constraint (Lines 39-42)

  • create a ROL optimization prolem (Lines 45-58), and solve it (Line 60). Note that the volume equality constraint is imposed in Line 58.

 1import firedrake as fd
 2import fireshape as fs
 3import fireshape.zoo as fsz
 4import ROL
 5from PDEconstraint_pipe import NavierStokesSolver
 6from objective_pipe import PipeObjective
 7
 8# setup problem
 9mesh = fd.Mesh("pipe.msh")
10Q = fs.FeControlSpace(mesh)
11inner = fs.LaplaceInnerProduct(Q, fixed_bids=[10, 11, 12])
12q = fs.ControlVector(Q, inner)
13
14# setup PDE constraint
15if mesh.topological_dimension() == 2:  # in 2D
16    viscosity = fd.Constant(1./400.)
17elif mesh.topological_dimension() == 3:  # in 3D
18    viscosity = fd.Constant(1/10.)  # simpler problem in 3D
19else:
20    raise NotImplementedError
21e = NavierStokesSolver(Q.mesh_m, viscosity)
22
23# save state variable evolution in file u2.pvd or u3.pvd
24if mesh.topological_dimension() == 2:  # in 2D
25    out = fd.File("solution/u2D.pvd")
26elif mesh.topological_dimension() == 3:  # in 3D
27    out = fd.File("solution/u3D.pvd")
28
29
30def cb():
31    return out.write(e.solution.split()[0])
32
33
34# create PDEconstrained objective functional
35J_ = PipeObjective(e, Q, cb=cb)
36J = fs.ReducedObjective(J_, e)
37
38# add regularization to improve mesh quality
39Jq = fsz.MoYoSpectralConstraint(10, fd.Constant(0.5), Q)
40J = J + Jq
41
42# Set up volume constraint
43vol = fsz.VolumeFunctional(Q)
44initial_vol = vol.value(q, None)
45econ = fs.EqualityConstraint([vol], target_value=[initial_vol])
46emul = ROL.StdVector(1)
47
48# ROL parameters
49params_dict = {
50    'General': {'Print Verbosity': 0,  # set to 1 to understand output
51                'Secant': {'Type': 'Limited-Memory BFGS',
52                           'Maximum Storage': 10}},
53    'Step': {'Type': 'Augmented Lagrangian',
54             'Augmented Lagrangian':
55             {'Subproblem Step Type': 'Trust Region',
56              'Print Intermediate Optimization History': False,
57              'Subproblem Iteration Limit': 10}},
58    'Status Test': {'Gradient Tolerance': 1e-2,
59                    'Step Tolerance': 1e-2,
60                    'Constraint Tolerance': 1e-1,
61                    'Iteration Limit': 10}}
62params = ROL.ParameterList(params_dict, "Parameters")
63problem = ROL.OptimizationProblem(J, q, econ=econ, emul=emul)
64solver = ROL.OptimizationSolver(problem, params)
65solver.solve()

Note

This problem can also be solved using Bsplines to discretize the control. For instance, one could replace Line 10-11 with

bbox = [(1.5, 12.), (0, 6.)]
orders = [4, 4]
levels = [4, 3]
Q = fs.BsplineControlSpace(mesh, bbox, orders, levels, boundary_regularities=[2, 0])
inner = fs.H1InnerProduct(Q)

In this case, the control is discretized using tensorized cubic (order = [4, 4]) Bsplines (roughly \(2^4\) in the \(x\)-direction \(\times\, 2^3\) in the \(y\)-direction; levels = [4, 3]). These Bsplines lie in the box with lower left corner \((1.5, 0)\) and upper right corner \((12., 6.)\) (bbox = [(1.5, 12.), (0, 6.)]). With boundary_regularities = [2, 0] we prescribe that the transformation vanishes for \(x=1.5\) and \(x=12\) with \(C^1\)-regularity, but it does not necessarily vanish for \(y=0\) and \(y=6\). In light of this, we do not need to specify fixed_bids in the inner product.

Using Bsplines to discretize the control leads to similar results.

Result

pipe2D_final pipe3D_final

Optimized pipe design in 2D (left) and 3D (right).

For the 2D-example, typing python3 main_pipe.py in the terminal returns:

Augmented Lagrangian Solver
Subproblem Solver: Trust Region
  iter  fval           cnorm          gLnorm         snorm          penalty   feasTol   optTol    #fval   #grad   #cval   subIter
  0     4.390787e-01   0.000000e+00   4.690997e-01                  1.00e+01  1.26e-01  1.00e-02
  1     2.823114e-01   6.377197e-01   7.192107e-02   9.086354e-01   1.00e+02  1.26e-01  1.00e-01  15      11      22      10
  2     3.266521e-01   9.002146e-03   2.942106e-01   4.741487e-01   1.00e+02  7.94e-02  1.00e-03  23      15      32      5
  3     3.256262e-01   3.716964e-05   9.015756e-02   7.546472e-02   1.00e+02  5.01e-02  1.00e-04  36      25      53      10
  4     3.251405e-01   6.469895e-06   5.961345e-02   4.412538e-02   1.00e+02  3.16e-02  1.00e-04  49      35      74      10
  5     3.249782e-01   4.207885e-04   6.196224e-02   6.207929e-02   1.00e+02  2.00e-02  1.00e-04  62      45      95      10
  6     3.248956e-01   1.972863e-05   8.589613e-03   9.381617e-03   1.00e+02  1.26e-02  1.00e-04  75      55      116     10
Optimization Terminated with Status: Converged

Note

To store the terminal output in a txt file, use the bash command python3 main_pipe.py >> output.txt.

We can inspect the result by opening the file u.pvd with ParaView. We see that the difference between the volume of the initial guess and of the retrieved optimized design is roughly \(2\cdot 10^{-5}\).

For the 3D-example, typing python3 main_pipe.py in the terminal returns:

Augmented Lagrangian Solver
Subproblem Solver: Trust Region
  iter  fval           cnorm          gLnorm         snorm          penalty   feasTol   optTol    #fval   #grad   #cval   subIter
  0     1.211321e+01   0.000000e+00   1.000000e+00                  1.50e+01  1.00e-01  1.00e-02
  1     8.681223e+00   1.932438e+00   1.998650e-01   1.058693e+00   1.50e+02  8.72e-02  6.65e-02  15      13      24      10
  2     1.111633e+01   2.419067e-02   3.687839e-01   1.333067e+00   1.50e+02  5.28e-02  4.42e-04  23      19      36      5
  3     1.114521e+01   1.079619e-03   3.528500e-01   1.243047e-01   1.50e+02  3.20e-02  1.00e-04  36      29      57      10
  4     1.113536e+01   7.521605e-04   6.820438e-02   1.229308e-01   1.50e+02  1.94e-02  1.00e-04  49      39      78      10
  5     1.113373e+01   5.221684e-04   1.130481e-01   9.808759e-02   1.50e+02  1.17e-02  1.00e-04  62      49      99      10
  6     1.113360e+01   3.785105e-04   5.083579e-02   7.107051e-02   1.50e+02  7.11e-03  1.00e-04  75      59      120     10
  7     1.113168e+01   2.478858e-04   8.290102e-02   4.790457e-02   1.50e+02  4.30e-03  1.00e-04  88      69      141     10
  8     1.113127e+01   6.115143e-05   6.973287e-02   4.889087e-02   1.50e+02  2.61e-03  1.00e-04  101     79      162     10
  9     1.113080e+01   1.956210e-04   3.754803e-02   3.813478e-02   1.50e+02  1.58e-03  1.00e-04  114     89      183     10
  10    1.112986e+01   5.181766e-06   3.406590e-02   2.990831e-02   1.50e+02  1.00e-03  1.00e-04  127     99      204     10
Optimization Terminated with Status: Iteration Limit Exceeded

We can inspect the result by opening the file u.pvd with ParaView. We see that the difference between the volume of the initial guess and of the retrieved optimized design is roughly \(5\cdot 10^{-6}\).