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}\).

import firedrake as fd
import fireshape as fs


class LevelsetFunctional(fs.ShapeObjective):
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

        # physical mesh
        self.mesh_m = self.Q.mesh_m

        # global function defined in terms of physical coordinates
        x, y = fd.SpatialCoordinate(self.mesh_m)
        self.f = (x - 0.5)**2 + (y - 0.5)**2 - 0.5

    def value_form(self):
        # volume integral
        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).

import firedrake as fd
import fireshape as fs
import ROL
from levelsetfunctional import LevelsetFunctional

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

# save shape evolution in file domain.pvd
out = fd.File("domain.pvd")

# create objective functional
J = LevelsetFunctional(Q, cb=lambda: out.write(Q.mesh_m.coordinates))

# ROL parameters
params_dict = {
    'Step': {
        'Type': 'Line Search',
        'Line Search': {
            'Descent Method': {
                'Type': 'Quasi-Newton Step'
            }
        }
    },
    'General': {
        'Secant': {
            'Type': 'Limited-Memory BFGS',
            'Maximum Storage': 25
        }
    },
    'Status Test': {
        'Gradient Tolerance': 1e-4,
        'Step Tolerance': 1e-10,
        'Iteration Limit': 30
    }
}

params = ROL.ParameterList(params_dict, "Parameters")
problem = ROL.OptimizationProblem(J, q)
solver = ROL.OptimizationSolver(problem, params)
solver.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.