Example 2: L2-tracking

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

\[\mathcal{J}(\Omega) = \int_\Omega \big(u(\mathbf{x}) - u_t(\mathbf{x})\big)^2 \,\mathrm{d}\mathbf{x}\,.\]

where \(u:\mathbb{R}^2\to\mathbb{R}\) is the solution to the scalar boundary value problem

\[-\Delta u = 4 \quad \text{in }\Omega\,, \qquad u = 0 \quad \text{on } \partial\Omega\]

and \(u_t:\mathbb{R}^2\to\mathbb{R}\) is a target function. In particular, we consider

\[u_t(x,y) = 0.36 - (x - 0.5)^2 + (y - 0.5)^2 - 0.5\,.\]

Beside the empty set, the domain that minimizes this shape functional is a disc of radius \(0.6\) centered at \((0.5,0.5)\).

Implementing the PDE constraint

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

import firedrake as fd
from fireshape import PdeConstraint


class PoissonSolver(PdeConstraint):
    """A Poisson BVP with hom DirBC as PDE constraint."""
    def __init__(self, mesh_m):
        super().__init__()

        # Setup problem
        V = fd.FunctionSpace(mesh_m, "CG", 1)

        # Weak form of Poisson problem
        u = fd.Function(V, name="State")
        v = fd.TestFunction(V)
        f = fd.Constant(4.)
        self.F = (fd.inner(fd.grad(u), fd.grad(v)) - f * v) * fd.dx
        self.bcs = fd.DirichletBC(V, 0., "on_boundary")

        # PDE-solver parameters
        self.params = {
            "ksp_type": "cg",
            "mat_type": "aij",
            "pc_type": "hypre",
            "pc_factor_mat_solver_package": "boomerang",
            "ksp_rtol": 1e-11,
            "ksp_atol": 1e-11,
            "ksp_stol": 1e-15,
        }

        self.solution = u

        # problem = fd.NonlinearVariationalProblem(
        #     self.F, self.solution, bcs=self.bcs)
        # self.solver = fd.NonlinearVariationalSolver(
        #     problem, solver_parameters=self.params)

    def solve(self):
        super().solve()
        fd.solve(self.F == 0, self.solution, bcs=self.bcs,
                 solver_parameters=self.params)
        # self.solver.solve()

Note

To solve the discretized variational problem, we use CG with a multigrid preconditioner (see self.params in Lines 22-30). For 2D problems, one can also use direct solvers.

Implementing the shape functional

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

import firedrake as fd
from fireshape import ShapeObjective
from L2tracking_PDEconstraint import PoissonSolver


class L2trackingObjective(ShapeObjective):
    """L2 tracking functional for Poisson problem."""

    def __init__(self, pde_solver: PoissonSolver, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.u = pde_solver.solution

        # target function, exact soln is disc of radius 0.6 centered at
        # (0.5,0.5)
        (x, y) = fd.SpatialCoordinate(self.Q.mesh_m)
        self.u_target = 0.36 - (x-0.5)*(x-0.5) - (y-0.5)*(y-0.5)

    def value_form(self):
        """Evaluate misfit functional."""
        u = self.u
        return (u - self.u_target)**2 * fd.dx

Setting up and solving the problem

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

To set up the problem, we

• construct the mesh of the initial guess (Line 8, a unit square with lower left corner in the origin)

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

• choose the metric of the control space (Line 10, based on linear elasticity energy norm)

• initialize the PDE contraint on the physical mesh mesh_m (Line 15)

• specify to save the function \(u\) after each iteration in the file u.pvd by setting the function cb appropriately (Lines 19 and 22).

• initialize the shape functional (Line 22), and the reduce shape functional (Line 23),

• create a ROL optimization prolem (Lines 26-49), and solve it (Line 50).

import firedrake as fd
import fireshape as fs
import ROL
from L2tracking_PDEconstraint import PoissonSolver
from L2tracking_objective import L2trackingObjective

# setup problem
mesh = fd.UnitSquareMesh(100, 100)
Q = fs.FeControlSpace(mesh)
inner = fs.ElasticityInnerProduct(Q)
q = fs.ControlVector(Q, inner)

# setup PDE constraint
mesh_m = Q.mesh_m
e = PoissonSolver(mesh_m)

# save state variable evolution in file u.pvd
e.solve()
out = fd.File("u.pvd")

# create PDEconstrained objective functional
J_ = L2trackingObjective(e, Q, cb=lambda: out.write(e.solution))
J = fs.ReducedObjective(J_, e)

# ROL parameters
params_dict = {
    'Step': {
        'Type': 'Line Search',
        'Line Search': {
            'Descent Method': {
                'Type': 'Quasi-Newton Step'
            }
        },
    },
    'General': {
        'Secant': {
            'Type': 'Limited-Memory BFGS',
            'Maximum Storage': 10
        }
    },
    'Status Test': {
        'Gradient Tolerance': 1e-4,
        'Step Tolerance': 1e-5,
        'Iteration Limit': 15
    }
}
params = ROL.ParameterList(params_dict, "Parameters")
problem = ROL.OptimizationProblem(J, q)
solver = ROL.OptimizationSolver(problem, params)
solver.solve()

Result

Typing python3 L2tracking_main.py in the terminal returns:

Dogleg Trust-Region Solver with Limited-Memory BFGS Hessian Approximation
  iter  value          gnorm          snorm          delta          #fval     #grad     tr_flag
  0     3.984416e-03   4.253880e-02                  4.253880e-02
  1     2.380406e-03   3.243635e-02   4.253880e-02   1.063470e-01   3         2         0
  2     8.449408e-04   1.731817e-02   1.063470e-01   1.063470e-01   4         3         0
  3     4.944712e-04   8.530990e-03   2.919058e-02   2.658675e-01   5         4         0
  4     1.783406e-04   5.042658e-03   5.182347e-02   6.646687e-01   6         5         0
  5     2.395524e-05   1.342447e-03   5.472223e-02   1.661672e+00   7         6         0
  6     6.994156e-06   4.090116e-04   2.218345e-02   4.154180e+00   8         7         0
  7     4.011765e-06   2.961338e-04   1.009231e-02   1.038545e+01   9         8         0
  8     1.170509e-06   2.423493e-04   1.715142e-02   2.596362e+01   10        9         0
  9     1.170509e-06   2.423493e-04   1.372960e-02   3.432399e-03   11        9         2
  10    1.085373e-06   4.112702e-04   3.432399e-03   3.432399e-03   12        10        0
  11    7.428449e-07   1.090433e-04   3.432399e-03   8.580998e-03   13        11        0
  12    6.755854e-07   8.358801e-05   2.034473e-03   2.145249e-02   14        12        0
Optimization Terminated with Status: Converged

We can inspect the result by opening the file u.pvd with ParaView.