This tutorial is part of multiphenicsx.

Tutorial 06, case 3b: advection diffusion reaction control problem with Neumann control

In this tutorial we solve the optimal control problem

$$\min J(y, u) = \frac{1}{2} \int_{\Omega_3} (y - y_d)^2 dx + \frac{\alpha}{2} \int_{\Gamma_2} u^2 ds$$ s.t. $$\begin{cases} - \epsilon \Delta y + \beta \cdot \nabla y + \sigma y = f & \text{in } \Omega\\ y = g & \text{on } \Gamma_1\\ \epsilon \partial_n y = u & \text{on } \Gamma_2\\ \epsilon \partial_n y = 0 & \text{on } \Gamma_3\\ \end{cases}$$

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

The test case is from section 5.3 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.

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 petsc4py.PETSc
import ufl
import viskex

import multiphenicsx.fem
import multiphenicsx.fem.petsc

Geometrical parameters

mesh_size = 0.05

Mesh

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, [0.3, 0.7, 1.0])
[l4, l5, l6, _] = generate_rectangle_lines(1.0, 3.0, 0.0, 0.3)
[_, l7, l8, _] = generate_rectangle_lines(1.0, 3.0, 0.3, 0.7)
[_, l9, l10, _] = generate_rectangle_lines(1.0, 3.0, 0.7, 1.0)
line_loop_rectangle_left = gmsh.model.geo.addCurveLoop([l0, l1[0], l1[1], l1[2], l2, l3])
line_loop_rectangle_right_bottom = gmsh.model.geo.addCurveLoop([l4, l5[0], l6, -l1[0]])
line_loop_rectangle_right_middle = gmsh.model.geo.addCurveLoop([-l6, l7[0], l8, -l1[1]])
line_loop_rectangle_right_top = gmsh.model.geo.addCurveLoop([-l8, l9[0], l10, -l1[2]])
rectangle_left = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_left])
rectangle_right_bottom = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_right_bottom])
rectangle_right_middle = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_right_middle])
rectangle_right_top = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_right_top])

gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [l0, l2, l3], 1)
gmsh.model.addPhysicalGroup(1, [l4, l10], 2)
gmsh.model.addPhysicalGroup(1, [l5[0], l7[0], l9[0]], 3)
gmsh.model.addPhysicalGroup(2, [rectangle_left], 1)
gmsh.model.addPhysicalGroup(2, [rectangle_right_bottom, rectangle_right_top], 3)
gmsh.model.addPhysicalGroup(2, [rectangle_right_middle], 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    : [ 30%] Meshing curve 4 (Line)
Info    : [ 40%] Meshing curve 5 (Line)
Info    : [ 40%] Meshing curve 6 (Line)
Info    : [ 50%] Meshing curve 7 (Line)
Info    : [ 60%] Meshing curve 8 (Line)
Info    : [ 70%] Meshing curve 9 (Line)
Info    : [ 70%] Meshing curve 10 (Line)
Info    : [ 80%] Meshing curve 11 (Line)
Info    : [ 90%] Meshing curve 12 (Line)
Info    : [100%] Meshing curve 13 (Line)
Info    : Done meshing 1D (Wall 0.000947522s, CPU 0.001356s)
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.0291876s, CPU 0.029081s)
Info    : 1522 nodes 3152 elements

mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()

# 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]
boundaries_2 = boundaries.indices[boundaries.values == 2]

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

viskex.dolfinx.plot_mesh(mesh)

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

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

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

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

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Function spaces

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

Restrictions

dofs_Y = np.arange(0, Y.dofmap.index_map.size_local + Y.dofmap.index_map.num_ghosts)
dofs_U = dolfinx.fem.locate_dofs_topological(U, boundaries.dim, boundaries_2)
dofs_Q = dofs_Y
restriction_Y = multiphenicsx.fem.DofMapRestriction(Y.dofmap, dofs_Y)
restriction_U = multiphenicsx.fem.DofMapRestriction(U.dofmap, dofs_U)
restriction_Q = multiphenicsx.fem.DofMapRestriction(Q.dofmap, dofs_Q)
restriction = [restriction_Y, restriction_U, restriction_Q]

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.07
y_d = 2.5
epsilon = 1. / 12.
x = ufl.SpatialCoordinate(mesh)
beta = ufl.as_vector((x[1] * (1 - x[1]), 0))
sigma = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0))
ff = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0))
bc0 = petsc4py.PETSc.ScalarType(0)
bc1 = petsc4py.PETSc.ScalarType(1)

Optimality conditions

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(3), None, adjoint_operator],
     [None, alpha * ufl.inner(u, v) * ds(2), - ufl.inner(p, v) * ds(2)],
     [state_operator, - ufl.inner(u, q) * ds(2), None]]
f = [ufl.inner(y_d, z) * dx(3),
     None,
     ufl.inner(ff, q) * dx]
a[0][0] += dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)) * ufl.inner(y, z) * dx
a[2][2] = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)) * ufl.inner(p, q) * dx
f[1] = ufl.inner(dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)), v) * dx
a_cpp = dolfinx.fem.form(a)
f_cpp = dolfinx.fem.form(f)
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(bc1, bdofs_Y_1, Y),
      dolfinx.fem.dirichletbc(bc0, 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, y - y_d) * dx(3) + 0.5 * alpha * ufl.inner(u, u) * ds(2)
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})
a_state_cpp = dolfinx.fem.form(a_state)
f_state_cpp = dolfinx.fem.form(f_state)
bc_state = [bc[0]]

# Assemble the linear system for the state
A_state = dolfinx.fem.petsc.assemble_matrix(a_state_cpp, bcs=bc_state)
A_state.assemble()
F_state = dolfinx.fem.petsc.assemble_vector(f_state_cpp)
dolfinx.fem.petsc.apply_lifting(F_state, [a_state_cpp], [bc_state])
F_state.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(F_state, bc_state)

# Solve
ksp = petsc4py.PETSc.KSP()
ksp.create(mesh.comm)
ksp.setOperators(A_state)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(F_state, y.x.petsc_vec)
y.x.petsc_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
ksp.destroy()

<petsc4py.PETSc.KSP at 0x7fa9243400e0>

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, 1.35)

Uncontrolled J = 1.349999999999998

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

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Optimal control

# Assemble the block linear system for the optimality conditions
A = multiphenicsx.fem.petsc.assemble_matrix_block(a_cpp, bcs=bc, restriction=(restriction, restriction))
A.assemble()
F = multiphenicsx.fem.petsc.assemble_vector_block(f_cpp, a_cpp, bcs=bc, restriction=restriction)

# Solve
yup = multiphenicsx.fem.petsc.create_vector_block(f_cpp, restriction=restriction)
ksp = petsc4py.PETSc.KSP()
ksp.create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(F, yup)
yup.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
ksp.destroy()

<petsc4py.PETSc.KSP at 0x7fa91012aa20>

# Split the block solution in components
with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(yup, [Y.dofmap, U.dofmap, Q.dofmap], restriction) as yup_wrapper:
    for yup_wrapper_local, component in zip(yup_wrapper, (y, u, p)):
        with component.x.petsc_vec.localForm() as component_local:
            component_local[:] = yup_wrapper_local

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.035933676)

Optimal J = 0.03593400081778442

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

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

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

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

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

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen