Example 2: Dido’s problem

In this example, we solve Dido’s problem: find the two-dimensional geometry with biggest area for a given perimeter. Mathematically, this means minimizing the negative area

\[J(\Omega) = \int_\Omega (-1) \,\mathrm{d}\mathbf{x}\]

while keeping the perimeter

\[P(\Omega) = \int_\Omega (-1) \,\mathrm{d}s\]

constant. The solution to the problem is a disc with radius \(r = P(\Omega)/(2\pi)\). In this example, the initial domain is a unit square, which implies that \(P(\Omega)= 4\). We thus aim to reconstruct a disc of radius \(r = 2/\pi \approx 0.636619772367581\) and area \(r^2\pi \approx 1.273239544735163\).

In the following, we describe how to solve this problem in Fireshape. The entire script is contained in the Python file dido.py, which is saved in the Fireshape repository and can be found at the following link: link-to-dido-example.

Import modules

We begin by importing Firedrake, Fireshape, and ROL. We also import fireshop.zoo, which contains utilities to set up the perimeter constraint.

from firedrake import *
from fireshape import *
import fireshape.zoo as fsz
import ROL

Implement the shape function

To implement the shape function \(J\), we use Fireshape’s class PDEconstrainedObjective. This requires specifying how to evaluate \(J\) in the method PDEconstrainedObjective.objective_value.

class NegativeArea(PDEconstrainedObjective):
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

    def objective_value(self):
        return assemble(Constant(-1) * dx(self.Q.mesh_m))

Note

Although \(J\) is not technically constrained to a boundary value problem, it is convenient to use the class PDEconstrainedObjective as this automatically returns NaN on poor quality meshes.

Select initial guess, control space, and inner product

We select a unit square centered at \((0.5,0.5)\) as initial domain. To modify the domain, we create a control space of geometric transformations discretized using finite elements of degree 2. (Note the additional input add_to_degree_r=1 in the definition of Q.) To compute descent directions, we employ Riesz representatives of shape derivatives with respect to a full \(H^1\)-inner product. With these, we create a control variable q that will be updated by the optimization algorithm.

# Select initial guess, control space, and inner product
mesh = UnitSquareMesh(5, 5)
Q = FeControlSpace(mesh, add_to_degree_r=1)
IP = H1InnerProduct(Q)
q = ControlVector(Q, IP)

Instantiate objective function J

We instantiate \(J(\Omega)\) using the class NegativeArea we have created. During instantiation, we also pass a call back function cb that stores the shape iterates whenever J is evaluated.

# Instantiate objective function J
out = VTKFile("domain.pvd")
J = NegativeArea(Q, cb=lambda: out.write(Q.mesh_m.coordinates))

Instantiate equality constraint P

We instantiate \(P(\Omega)\) using the class SurfaceAreaFunctional from fireshop.zoo.

# Set up perimeter constraint
perimeter = fsz.SurfaceAreaFunctional(Q)
initial_perimeter = perimeter.value(q, None)
econ = EqualityConstraint([perimeter], target_value=[initial_perimeter])
emul = ROL.StdVector(1)

Select the optimization algorithm and solve the problem

Finally, we select an augmented Lagrangian optimization algorithm that uses a trust-region optimization algorithm (with l-BFGS Hessian) to solve the augmented Lagrangian subproblems. We carefully select and set the optimization parameters in the dictionary pd. This, together with J, q, econ, and emul are passed to ROL, which solves the problem.

# Select the optimization algorithm and solve the problem
pd = {'General': {'Print Verbosity': 0,
                  'Secant': {'Type': 'Limited-Memory BFGS',
                             'Maximum Storage': 10}},
      'Step': {'Type': 'Augmented Lagrangian',
               'Augmented Lagrangian':
               {'Use Default Problem Scaling': False,
                'Constraint Scaling': 1.5,
                'Subproblem Step Type': 'Trust Region',
                'Print Intermediate Optimization History': False,
                'Subproblem Iteration Limit': 10}},
      'Status Test': {'Gradient Tolerance': 1e-2,
                      'Step Tolerance': 1e-3,
                      'Constraint Tolerance': 1e-1,
                      'Iteration Limit': 10}}
params = ROL.ParameterList(pd, "Parameters")
problem = ROL.OptimizationProblem(J, q, econ=econ, emul=emul)
solver = ROL.OptimizationSolver(problem, params)
solver.solve()

Result

Typing python3 dido.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.000000e+00  0.000000e+00   1.359583e+00                  1.00e+01  1.26e-01  1.36e-02
  1     -1.477728e+00  3.093114e-01   3.244788e-02   4.698565e-01   1.00e+02  1.26e-01  1.00e-01  15      13      24      10
  2     -1.274740e+00  2.943785e-03   1.042650e-01   5.678527e-01   1.00e+02  7.94e-02  1.00e-01  28      23      45      10
  3     -1.272798e+00  1.229090e-04   7.903207e-02   1.968245e-02   1.00e+02  5.01e-02  1.00e-03  38      28      58      7
  4     -1.273378e+00  3.158611e-04   1.664273e-01   5.434476e-02   1.00e+02  3.16e-02  1.00e-03  51      39      80      10
  5     -1.273019e+00  3.230874e-04   1.760829e-02   5.272820e-02   1.00e+02  2.00e-02  1.00e-03  64      49      101     10
  6     -1.273236e+00  4.254024e-06   7.637262e-03   1.757102e-03   1.00e+02  1.26e-02  1.00e-03  77      59      122     10
Optimization Terminated with Status: Converged

This output implies that the retrieved optimal negative area is approximately -1.273236e+00 and that the (perimeter) equality constraint violation is approximately (only) 4.254024e-06.

We can inspect the result by opening the file levelset_domain.pvd with ParaView. In the GIF below, we see that the domain (black grid) converges to the right shape (red circle). Note that the mesh presents bent edges due to using finite elements of degree 2 to discretize domain updates.