Example 1: Level Set

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

\[\mathcal{J}(\Omega) = \int_\Omega f(\mathbf{x}) \,\mathrm{d}\mathbf{x}\,.\]

where \(f:\mathbb{R}^2\to\mathbb{R}\) is a scalar function. In particular, we consider

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

The domain that minimizes this shape functional is a disc of radius \(1/\sqrt{2}\) centered at \((0.5,0.5)\).

Implementing the shape functional

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

 1import firedrake as fd
 2import fireshape as fs
 3
 4
 5class LevelsetFunctional(fs.ShapeObjective):
 6    def __init__(self, *args, **kwargs):
 7        super().__init__(*args, **kwargs)
 8
 9        # physical mesh
10        self.mesh_m = self.Q.mesh_m
11
12        # global function defined in terms of physical coordinates
13        x, y = fd.SpatialCoordinate(self.mesh_m)
14        self.f = (x - 0.5)**2 + (y - 0.5)**2 - 0.5
15
16    def value_form(self):
17        # volume integral
18        return self.f * fd.dx

Setting up and solving the problem

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

To set up the problem, we need to:

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

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

  • choose the metric of the control space (Line 9, \(H^1\)-seminorm)

Then, we specify to save the mesh after each iteration in the file domain.pvd by setting the function cb appropriately (Lines 13 and 16).

Finally, we initialize the shape functional (Line 16), create a ROL optimization prolem (Lines 19-43), and solve it (Line 44).

 1import firedrake as fd
 2import fireshape as fs
 3import ROL
 4from levelsetfunctional import LevelsetFunctional
 5
 6# setup problem
 7mesh = fd.UnitSquareMesh(30, 30)
 8Q = fs.FeControlSpace(mesh)
 9inner = fs.LaplaceInnerProduct(Q)
10q = fs.ControlVector(Q, inner)
11
12# save shape evolution in file domain.pvd
13out = fd.File("domain.pvd")
14
15# create objective functional
16J = LevelsetFunctional(Q, cb=lambda: out.write(Q.mesh_m.coordinates))
17
18# ROL parameters
19params_dict = {
20    'Step': {
21        'Type': 'Line Search',
22        'Line Search': {
23            'Descent Method': {
24                'Type': 'Quasi-Newton Step'
25            }
26        }
27    },
28    'General': {
29        'Secant': {
30            'Type': 'Limited-Memory BFGS',
31            'Maximum Storage': 25
32        }
33    },
34    'Status Test': {
35        'Gradient Tolerance': 1e-4,
36        'Step Tolerance': 1e-10,
37        'Iteration Limit': 30
38    }
39}
40
41params = ROL.ParameterList(params_dict, "Parameters")
42problem = ROL.OptimizationProblem(J, q)
43solver = ROL.OptimizationSolver(problem, params)
44solver.solve()

Result

Typing python3 levelset.py in the terminal returns:

Quasi-Newton Method with Limited-Memory BFGS
Line Search: Cubic Interpolation satisfying Strong Wolfe Conditions
  iter  value          gnorm          snorm          #fval     #grad     ls_#fval  ls_#grad
  0     -3.333333e-01  2.426719e-01
  1     -3.765250e-01  1.121533e-01   2.426719e-01   2         2         1         0
  2     -3.886767e-01  8.040797e-02   2.525978e-01   3         3         1         0
  3     -3.919733e-01  2.331695e-02   6.622577e-02   4         4         1         0
  4     -3.926208e-01  4.331993e-03   5.628606e-02   5         5         1         0
  5     -3.926603e-01  2.906313e-03   1.239420e-02   6         6         1         0
  6     -3.926945e-01  9.456530e-04   2.100085e-02   7         7         1         0
  7     -3.926980e-01  3.102278e-04   6.952015e-03   8         8         1         0
  8     -3.926987e-01  1.778454e-04   3.840828e-03   9         9         1         0
  9     -3.926989e-01  9.001788e-05   2.672387e-03   10        10        1         0
Optimization Terminated with Status: Converged

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