Shape Optimization with a Poisson Problem

Problem Formulation

In this demo, we investigate the basics of CASHOCS for shape optimization problems. As a model problem, we investigate the following one from Etling, Herzog, Loayza, Wachsmuth, First and Second Order Shape Optimization Based on Restricted Mesh Deformations

\[\begin{split}&\min_\Omega J(u, \Omega) = \int_\Omega u \text{ d}x \\ &\text{subject to} \quad \left\lbrace \quad \begin{alignedat}{2} -\Delta u &= f \quad &&\text{ in } \Omega,\\ u &= 0 \quad &&\text{ on } \Gamma. \end{alignedat} \right.\end{split}\]

For the initial domain, we use the unit disc \(\Omega = \{ x \in \mathbb{R}^2 \,\mid\, \lvert\lvert x \rvert\rvert_2 < 1 \}\), and the right-hand side \(f\) is given by

\[f(x) = 2.5 \left( x_1 + 0.4 - x_2^2 \right)^2 + x_1^2 + x_2^2 - 1.\]

Implementation

The complete python code can be found in the file demo_shape_poisson.py, and the corresponding config can be found in config.ini.

Initialization

We start the problem by using a wildcard import for FEniCS, and by importing CASHOCS

from fenics import *
import cashocs

As for the case of optimal control problems, we can specify the verbosity of CASHOCS with the line

cashocs.set_log_level(cashocs.LogLevel.INFO)

which is documented at cashocs.set_log_level() (cf. Distributed Control of a Poisson Problem).

Similarly to the optimal control case, we also require config files for shape optimization problems in CASHOCS. A detailed discussion of the config files for shape optimization is given in Documentation of the Config Files for Shape Optimization Problems. We read the config file with the load_config command

config = cashocs.load_config('./config.ini')

Next, we have to define the mesh. As the above problem is posed on the unit disc initially, we define this via FEniCS commands (CASHOCS only has rectangular meshes built in). This is done via the following code

meshlevel = 10
degree = 1
dim = 2
mesh = UnitDiscMesh.create(MPI.comm_world, meshlevel, degree, dim)

Next up, we define the fenics.Measure objects, which we need to define the problem. For the volume measure, we can simply invoke

dx = Measure('dx', mesh)

However, for the surface measure, we need to mark the boundary. This is required since CASHOCS distinguishes between three types of boundaries: The deformable boundary, the fixed boundary, and boundaries that can only be deformed perpendicular to a certain coordinate axis (see the relevant documentation of the config files). Here, we investigate the case of a completely deformable boundary, which makes things rather easy. We mark this boundary with the marker 1 with the following piece of code

boundary = CompiledSubDomain('on_boundary')
boundaries = MeshFunction('size_t', mesh, dim=1)
boundary.mark(boundaries, 1)
ds = Measure('ds', mesh, subdomain_data=boundaries)

Note

In config.ini, in the section ShapeGradient, there is the line

shape_bdry_def = [1]

which specifies that the boundary marked with 1 is deformable. For our example this is exactly what we want, as this means that the entire boundary is variable, due to the previous commands. For a detailed documentation we refer to the corresponding documentation of the ShapeGradient section.

Note, that all of the alternative ways of marking subdomains or boundaries with numbers, as explained in Langtangen and Logg, Solving PDEs in Python also work here. If it is valid for FEniCS, it is also for CASHOCS.

After having defined the initial geometry, we define a fenics.FunctionSpace consisting of piecewise linear Lagrange elements via

V = FunctionSpace(mesh, 'CG', 1)
u = Function(V)
p = Function(V)

This also defines our state variable \(u\) as u, and the adjoint state \(p\) is given by p.

Note

As remarked in Distributed Control of a Poisson Problem, in classical FEniCS syntax we would use a fenics.TrialFunction for u and a fenics.TestFunction for p. However, for CASHOCS this must not be the case. Instead, the state and adjoint variables have to be fenics.Function objects.

The right-hand side of the PDE constraint is then defined as

x = SpatialCoordinate(mesh)
f = 2.5*pow(x[0] + 0.4 - pow(x[1], 2), 2) + pow(x[0], 2) + pow(x[1], 2) - 1

which allows us to define the weak form of the state equation via

e = inner(grad(u), grad(p))*dx - f*p*dx
bcs = DirichletBC(V, Constant(0), boundaries, 1)

The optimization problem and its solution

We are now almost done, the only thing left to do is to define the cost functional

J = u*dx

and the shape optimization problem

sop = cashocs.ShapeOptimizationProblem(e, bcs, J, u, p, boundaries, config)

This can then be solved in complete analogy to Distributed Control of a Poisson Problem with the sop.solve() command

sop.solve()

The result of the optimization looks like this

../../_images/img_shape_poisson.png

Note

As in Distributed Control of a Poisson Problem we can specify some keyword arguments for the solve command. If none are given, then the settings from the config file are used, but if some are given, they override the parameters specified in the config file. In particular, these arguments are

  • algorithm : Specifies which solution algorithm shall be used.

  • rtol : The relative tolerance for the optimization algorithm.

  • atol : The absolute tolerance for the optimization algorithm.

  • max_iter : The maximum amount of iterations that can be carried out.

The possible choices for these parameters are discussed in detail in Section OptimizationRoutine and the documentation of the solve method.

As before, it is not strictly necessary to supply config files to CASHOCS, but it is very strongly recommended to do so. In case one does not supply a config file, one has to at least specify the solution algorithm in the call to the solve method.