This tutorial is part of multiphenicsx.

Tutorial 06, case 1b: Poisson problem with distributed control

In this tutorial we solve the optimal control problem

$$\min J(y, u) = \frac{1}{2} \int_{\Omega} (y - y_d)^2 dx + \frac{\alpha}{2} \int_{\Omega} u^2 dx$$ s.t. $$\begin{cases} - \Delta y = f + u & \text{in } \Omega\\ y = 1 & \text{on } \partial\Omega \end{cases}$$

where $$\begin{align*} & \Omega & \text{unit square}\\ & u \in L^2(\Omega) & \text{control variable}\\ & y \in H^1(\Omega) & \text{state variable}\\ & \alpha > 0 & \text{penalization parameter}\\ & y_d & \text{a piecewise constant desired state}\\ & f & \text{forcing term} \end{align*}$$ using an adjoint formulation solved by a one shot approach.

The test case is from section 5.1 of

F. Negri, G. Rozza, A. Manzoni and A. Quarteroni. Reduced Basis Method for Parametrized Elliptic Optimal Control Problems. SIAM Journal on Scientific Computing, 35(5): A2316-A2340, 2013.

Note that this case does not really need multiphenicsx, and can be run with just dolfinx.

import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.mesh
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex

Geometrical parameters

L1 = 1.0
L2 = 3.0
H = 1.0
mesh_size = 0.05

Mesh

gmsh.initialize()
gmsh.model.add("mesh")

p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, mesh_size)
p1 = gmsh.model.geo.addPoint(L1, 0.0, 0.0, mesh_size)
p2 = gmsh.model.geo.addPoint(L1 + L2, 0.0, 0.0, mesh_size)
p3 = gmsh.model.geo.addPoint(L1 + L2, H, 0.0, mesh_size)
p4 = gmsh.model.geo.addPoint(L1, H, 0.0, mesh_size)
p5 = gmsh.model.geo.addPoint(0.0, H, 0.0, mesh_size)
l0 = gmsh.model.geo.addLine(p0, p1)
l1 = gmsh.model.geo.addLine(p1, p4)
l2 = gmsh.model.geo.addLine(p4, p5)
l3 = gmsh.model.geo.addLine(p5, p0)
l4 = gmsh.model.geo.addLine(p1, p2)
l5 = gmsh.model.geo.addLine(p2, p3)
l6 = gmsh.model.geo.addLine(p3, p4)
line_loop_rectangle_left = gmsh.model.geo.addCurveLoop([l0, l1, l2, l3])
line_loop_rectangle_right = gmsh.model.geo.addCurveLoop([l4, l5, l6, -l1])
rectangle_left = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_left])
rectangle_right = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_right])

gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [l0, l4, l5, l6, l2, l3], 1)
gmsh.model.addPhysicalGroup(2, [rectangle_left], 1)
gmsh.model.addPhysicalGroup(2, [rectangle_right], 2)
gmsh.model.mesh.generate(2)

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 20%] Meshing curve 2 (Line)
Info    : [ 30%] Meshing curve 3 (Line)
Info    : [ 50%] Meshing curve 4 (Line)
Info    : [ 60%] Meshing curve 5 (Line)
Info    : [ 80%] Meshing curve 6 (Line)
Info    : [ 90%] Meshing curve 7 (Line)
Info    : Done meshing 1D (Wall 0.000641894s, CPU 0.000517s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : [ 60%] Meshing surface 2 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.0401248s, CPU 0.039219s)
Info    : 1968 nodes 3960 elements

mesh, subdomains, boundaries, *_ = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()
assert subdomains is not None
assert boundaries is not None

# Create connectivities required by the rest of the code
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)

boundaries_1 = boundaries.indices[boundaries.values == 1]

# Define associated measures
dx = ufl.Measure("dx", subdomain_data=subdomains)

viskex.dolfinx.plot_mesh(mesh)

viskex.dolfinx.plot_mesh_tags(mesh, subdomains, "subdomains")

viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")

Function spaces

Y = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
U = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
Q = Y.clone()

Trial and test functions

(y, u, p) = (ufl.TrialFunction(Y), ufl.TrialFunction(U), ufl.TrialFunction(Q))
(z, v, q) = (ufl.TestFunction(Y), ufl.TestFunction(U), ufl.TestFunction(Q))

Problem data

alpha = 0.01
y_d_1 = 1.0
y_d_2 = 0.6
zero =  petsc4py.PETSc.ScalarType(0)  # type: ignore[attr-defined]
one =  petsc4py.PETSc.ScalarType(1)  # type: ignore[attr-defined]
ff = dolfinx.fem.Constant(mesh, zero)

Optimality conditions

a = [[ufl.inner(y, z) * dx, None, ufl.inner(ufl.grad(p), ufl.grad(z)) * dx],
     [None, alpha * ufl.inner(u, v) * dx, - ufl.inner(p, v) * dx],
     [ufl.inner(ufl.grad(y), ufl.grad(q)) * dx, - ufl.inner(u, q) * dx, None]]
f = [ufl.inner(y_d_1, z) * dx(1) + ufl.inner(y_d_2, z) * dx(2),
     None,
     ufl.inner(ff, q) * dx]
a[2][2] = dolfinx.fem.Constant(mesh, zero) * ufl.inner(p, q) * dx
f[1] = ufl.inner(dolfinx.fem.Constant(mesh, zero), v) * dx
bdofs_Y_1 = dolfinx.fem.locate_dofs_topological(Y, mesh.topology.dim - 1, boundaries_1)
bdofs_Q_1 = dolfinx.fem.locate_dofs_topological(Q, mesh.topology.dim - 1, boundaries_1)
bc = [dolfinx.fem.dirichletbc(one, bdofs_Y_1, Y),
      dolfinx.fem.dirichletbc(zero, bdofs_Q_1, Q)]

Solution

(y, u, p) = (dolfinx.fem.Function(Y), dolfinx.fem.Function(U), dolfinx.fem.Function(Q))

Cost functional

J = (0.5 * ufl.inner(y - y_d_1, y - y_d_1) * dx(1) + 0.5 * ufl.inner(y - y_d_2, y - y_d_2) * dx(2)
     + 0.5 * alpha * ufl.inner(u, u) * dx)
J_cpp = dolfinx.fem.form(J)

Uncontrolled functional value

# Extract state forms from the optimality conditions
a_state = ufl.replace(a[2][0], {q: z})
f_state = ufl.replace(f[2], {q: z})
bc_state = [bc[0]]

petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "ksp_error_if_not_converged": True,
}
problem_state = dolfinx.fem.petsc.LinearProblem(
    a_state, f_state, bcs=bc_state, u=y,
    petsc_options_prefix="tutorial_1b_poisson_state_", petsc_options=petsc_options
)
problem_state.solve()
del problem_state

J_uncontrolled = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(J_cpp), op=mpi4py.MPI.SUM)
print("Uncontrolled J =", J_uncontrolled)
assert np.isclose(J_uncontrolled, 0.24)

Uncontrolled J = 0.24000000000000263

viskex.dolfinx.plot_scalar_field(y, "uncontrolled state")

Optimal control

problem = dolfinx.fem.petsc.LinearProblem(
    a, f, bcs=bc, u=(y, u, p),
    petsc_options_prefix="tutorial_1b_poisson_", petsc_options=petsc_options,
    kind="mpi"
)
problem.solve()
del problem

J_controlled = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(J_cpp), op=mpi4py.MPI.SUM)
print("Optimal J =", J_controlled)
assert np.isclose(J_controlled, 0.158485065)

Optimal J = 0.15848451131272584

viskex.dolfinx.plot_scalar_field(y, "state")

viskex.dolfinx.plot_scalar_field(u, "control")

viskex.dolfinx.plot_scalar_field(p, "adjoint")