In this tutorial we solve the optimal control problem
$$\min J(y, u) = \frac{1}{2} \int_{\Omega_1 \cup \Omega_2} (y - y_d)^2 dx + \frac{\alpha}{2} \int_{\Omega} u^2 dx$$ s.t. $$\begin{cases} - \epsilon \Delta y + \beta \cdot \nabla y + \sigma y = f + u & \text{in } \Omega\\ y = g & \text{on } \partial\Omega \end{cases}$$
where $$\begin{align*} & \Omega & \text{domain}\\ & u \in L^2(\Omega) & \text{control variable}\\ & y \in H^1_0(\Omega) & \text{state variable}\\ & \alpha > 0 & \text{penalization parameter}\\ & y_d & \text{desired state}\\ & \epsilon > 0 & \text{diffusion coefficient}\\ & \beta \in \mathbb{R}^2 & \text{advection field}\\ & \sigma > 0 & \text{reaction coefficient}\\ & f & \text{forcing term}\\ & g & \text{non homogeneous piecewise constant Dirichlet BC}\\ \end{align*}$$ using an adjoint formulation solved by a one shot approach.
The test case is from section 5.2 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 typing
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.mesh
import gmsh
import mpi4py.MPI
import numpy as np
import numpy.typing as npt
import petsc4py.PETSc
import ufl
import viskex
mesh_size = 0.05
class GenerateRectangleLines:
"""Generate a rectangle."""
def __init__(self) -> None:
self.points: dict[tuple[float, float], int] = dict()
self.lines: dict[tuple[int, int], int] = dict()
def add_point(self, x: float, y: float) -> int:
"""Add a point to gmsh, if not present already."""
key = (x, y)
try:
return self.points[key]
except KeyError:
p = gmsh.model.geo.addPoint(x, y, 0.0, mesh_size)
self.points[key] = p
return p # type: ignore[no-any-return]
def add_line(self, p0: int, p1: int) -> int:
"""Add a line to gmsh, if not present already."""
try:
return self.lines[(p0, p1)]
except KeyError:
l01 = gmsh.model.geo.addLine(p0, p1)
self.lines[(p0, p1)] = l01
self.lines[(p1, p0)] = -l01
return l01 # type: ignore[no-any-return]
def __call__(
self, x_min: float, x_max: float, y_min: float, y_max: typing.Union[float, list[float]]
) -> tuple[int, list[int], int, int]:
"""Add points and lines that define a rectangle with the provided coordinates."""
p0 = self.add_point(x_min, y_min)
p1 = self.add_point(x_max, y_min)
if isinstance(y_max, list):
p2 = [self.add_point(x_max, y) for y in y_max]
p3 = self.add_point(x_min, y_max[-1])
else:
p2 = [self.add_point(x_max, y_max)]
p3 = self.add_point(x_min, y_max)
l0 = self.add_line(p0, p1)
p1_p2 = [p1, *p2]
l1 = [self.add_line(p1_p2[i], p1_p2[i + 1]) for i in range(len(p2))]
l2 = self.add_line(p2[-1], p3)
l3 = self.add_line(p3, p0)
return (l0, l1, l2, l3)
gmsh.initialize()
gmsh.model.add("mesh")
generate_rectangle_lines = GenerateRectangleLines()
[l0, l1, l2, l3] = generate_rectangle_lines(0.0, 1.0, 0.0, 1.0)
[l4, l5, l6, _] = generate_rectangle_lines(1.0, 2.5, 0.0, [0.3, 0.7, 1.0])
[l7, l8, l9, l10] = generate_rectangle_lines(0.2, 0.8, 0.3, 0.7)
[l11, l12, l13, l14] = generate_rectangle_lines(1.2, 2.5, 0.3, 0.7)
line_loop_rectangle_outer_left = gmsh.model.geo.addCurveLoop([l0, l1[0], l2, l3])
line_loop_rectangle_outer_right = gmsh.model.geo.addCurveLoop([l4, l5[0], l5[1], l5[2], l6, -l1[0]])
line_loop_rectangle_inner_left = gmsh.model.geo.addCurveLoop([l7, l8[0], l9, l10])
line_loop_rectangle_inner_right = gmsh.model.geo.addCurveLoop([l11, l12[0], l13, l14])
rectangle_outer_left = gmsh.model.geo.addPlaneSurface(
[line_loop_rectangle_outer_left, line_loop_rectangle_inner_left])
rectangle_outer_right = gmsh.model.geo.addPlaneSurface(
[line_loop_rectangle_outer_right, line_loop_rectangle_inner_right])
rectangle_inner_left = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_inner_left])
rectangle_inner_right = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_inner_right])
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [l0, l2, l3], 1)
gmsh.model.addPhysicalGroup(1, [l4, l6], 2)
gmsh.model.addPhysicalGroup(1, [l5[0], l5[1], l5[2]], 3)
gmsh.model.addPhysicalGroup(2, [rectangle_outer_left], 3)
gmsh.model.addPhysicalGroup(2, [rectangle_outer_right], 4)
gmsh.model.addPhysicalGroup(2, [rectangle_inner_left], 1)
gmsh.model.addPhysicalGroup(2, [rectangle_inner_right], 2)
gmsh.model.mesh.generate(2)
Info : Meshing 1D... Info : [ 0%] Meshing curve 1 (Line) Info : [ 10%] Meshing curve 2 (Line) Info : [ 20%] Meshing curve 3 (Line) Info : [ 20%] Meshing curve 4 (Line) Info : [ 30%] Meshing curve 5 (Line) Info : [ 40%] Meshing curve 6 (Line) Info : [ 40%] Meshing curve 7 (Line) Info : [ 50%] Meshing curve 8 (Line) Info : [ 60%] Meshing curve 9 (Line) Info : [ 60%] Meshing curve 10 (Line) Info : [ 70%] Meshing curve 11 (Line) Info : [ 70%] Meshing curve 12 (Line) Info : [ 80%] Meshing curve 13 (Line) Info : [ 90%] Meshing curve 14 (Line) Info : [ 90%] Meshing curve 15 (Line) Info : [100%] Meshing curve 16 (Line) Info : Done meshing 1D (Wall 0.00142417s, CPU 0.002076s) Info : Meshing 2D... Info : [ 0%] Meshing surface 1 (Plane, Frontal-Delaunay) Info : [ 30%] Meshing surface 2 (Plane, Frontal-Delaunay) Info : [ 60%] Meshing surface 3 (Plane, Frontal-Delaunay) Info : [ 80%] Meshing surface 4 (Plane, Frontal-Delaunay) Info : Done meshing 2D (Wall 0.0258163s, CPU 0.026103s) Info : 1301 nodes 2734 elements
mesh, subdomains, boundaries, *other_tags = 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)
# 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")
Y = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
U = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
Q = Y.clone()
(y, u, p) = (ufl.TrialFunction(Y), ufl.TrialFunction(U), ufl.TrialFunction(Q))
(z, v, q) = (ufl.TestFunction(Y), ufl.TestFunction(U), ufl.TestFunction(Q))
alpha = 0.01
y_d_1 = 0.6
y_d_2 = 1.8
epsilon = 1. / 15.
x = ufl.SpatialCoordinate(mesh)
beta = ufl.as_vector((x[1] * (1 - x[1]), 0))
zero = petsc4py.PETSc.ScalarType(0) # type: ignore[attr-defined]
sigma = dolfinx.fem.Constant(mesh, zero)
ff = dolfinx.fem.Constant(mesh, zero)
state_operator = (epsilon * ufl.inner(ufl.grad(y), ufl.grad(q)) * dx
+ ufl.inner(ufl.dot(beta, ufl.grad(y)), q) * dx + sigma * ufl.inner(y, q) * dx)
adjoint_operator = (epsilon * ufl.inner(ufl.grad(p), ufl.grad(z)) * dx
- ufl.inner(ufl.dot(beta, ufl.grad(p)), z) * dx + sigma * ufl.inner(p, z) * dx)
a = [[ufl.inner(y, z) * (dx(1) + dx(2)), None, adjoint_operator],
[None, alpha * ufl.inner(u, v) * dx, - ufl.inner(p, v) * dx],
[state_operator, - 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[0][0] += dolfinx.fem.Constant(mesh, zero) * ufl.inner(y, z) * 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
def bdofs_Y(idx: int) -> npt.NDArray[np.int32]:
"""Return DOFs of the space Y located on the boundary `idx`."""
assert boundaries is not None
return dolfinx.fem.locate_dofs_topological(
Y, mesh.topology.dim - 1, boundaries.indices[boundaries.values == idx])
def bdofs_Q(idx: int) -> npt.NDArray[np.int32]:
"""Return DOFs of the space Q located on the boundary `idx`."""
assert boundaries is not None
return dolfinx.fem.locate_dofs_topological(
Q, mesh.topology.dim - 1, boundaries.indices[boundaries.values == idx])
bc_state = [
dolfinx.fem.dirichletbc(petsc4py.PETSc.ScalarType(idx), bdofs_Y(idx), Y) # type: ignore[attr-defined]
for idx in (1, 2)
]
bc_adjoint = [dolfinx.fem.dirichletbc(zero, bdofs_Q(idx), Q) for idx in (1, 2)]
(y, u, p) = (dolfinx.fem.Function(Y), dolfinx.fem.Function(U), dolfinx.fem.Function(Q))
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)
# Extract state forms from the optimality conditions
a_state = ufl.replace(a[2][0], {q: z})
f_state = ufl.replace(f[2], {q: z})
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_2b_advection_diffusion_reaction_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.028096831)
Uncontrolled J = 0.028096639120195353
viskex.dolfinx.plot_scalar_field(y, "uncontrolled state")
problem = dolfinx.fem.petsc.LinearProblem(
a, f, bcs=[*bc_state, *bc_adjoint], u=(y, u, p),
petsc_options_prefix="tutorial_2b_advection_diffusion_reaction_", 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.001775304)
Optimal J = 0.0017753044656044726
viskex.dolfinx.plot_scalar_field(y, "state")
viskex.dolfinx.plot_scalar_field(u, "control")
viskex.dolfinx.plot_scalar_field(p, "adjoint")