This tutorial is part of multiphenicsx.

Tutorial 06, case 1a: 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 = 0 & \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_0(\Omega) & \text{state variable}\\ & \alpha > 0 & \text{penalization parameter}\\ & y_d & \text{desired state}\\ & f & \text{forcing term} \end{align*}$$ using an adjoint formulation solved by a one shot approach.

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 mpi4py.MPI
import numpy as np
import numpy.typing as npt
import petsc4py.PETSc
import ufl
import viskex

Mesh

mesh = dolfinx.mesh.create_unit_square(mpi4py.MPI.COMM_WORLD, 32, 32)

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

def bottom(x: npt.NDArray[np.float64]) -> npt.NDArray[np.bool_]:
    """Condition that defines the bottom boundary."""
    return abs(x[1] - 0.) < np.finfo(float).eps  # type: ignore[no-any-return]


def left(x: npt.NDArray[np.float64]) -> npt.NDArray[np.bool_]:
    """Condition that defines the left boundary."""
    return abs(x[0] - 0.) < np.finfo(float).eps  # type: ignore[no-any-return]


def top(x: npt.NDArray[np.float64]) -> npt.NDArray[np.bool_]:
    """Condition that defines the top boundary."""
    return abs(x[1] - 1.) < np.finfo(float).eps  # type: ignore[no-any-return]


def right(x: npt.NDArray[np.float64]) -> npt.NDArray[np.bool_]:
    """Condition that defines the right boundary."""
    return abs(x[0] - 1.) < np.finfo(float).eps  # type: ignore[no-any-return]


boundaries_entities = dict()
boundaries_values = dict()
for (boundary, boundary_id) in zip((bottom, left, top, right), (1, 2, 3, 4)):
    boundaries_entities[boundary_id] = dolfinx.mesh.locate_entities_boundary(
        mesh, mesh.topology.dim - 1, boundary)
    boundaries_values[boundary_id] = np.full(
        boundaries_entities[boundary_id].shape, boundary_id, dtype=np.int32)
boundaries_entities_unsorted = np.hstack(list(boundaries_entities.values()))
boundaries_values_unsorted = np.hstack(list(boundaries_values.values()))
boundaries_entities_argsort = np.argsort(boundaries_entities_unsorted)
boundaries_entities_sorted = boundaries_entities_unsorted[boundaries_entities_argsort]
boundaries_values_sorted = boundaries_values_unsorted[boundaries_entities_argsort]
boundaries = dolfinx.mesh.meshtags(
    mesh, mesh.topology.dim - 1,
    boundaries_entities_sorted, boundaries_values_sorted)
boundaries.name = "boundaries"

boundaries_1234 = boundaries.indices[np.isin(boundaries.values, (1, 2, 3, 4))]

viskex.dolfinx.plot_mesh(mesh)

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 = 1.e-5
x = ufl.SpatialCoordinate(mesh)
y_d = 10 * x[0] * (1 - x[0]) * x[1] * (1 - x[1])
zero =  petsc4py.PETSc.ScalarType(0)  # type: ignore[attr-defined]
ff = dolfinx.fem.Constant(mesh, zero)

Optimality conditions

a = [[ufl.inner(y, z) * ufl.dx, None, ufl.inner(ufl.grad(p), ufl.grad(z)) * ufl.dx],
     [None, alpha * ufl.inner(u, v) * ufl.dx, - ufl.inner(p, v) * ufl.dx],
     [ufl.inner(ufl.grad(y), ufl.grad(q)) * ufl.dx, - ufl.inner(u, q) * ufl.dx, None]]
f = [ufl.inner(y_d, z) * ufl.dx,
     None,
     ufl.inner(ff, q) * ufl.dx]
a[2][2] = dolfinx.fem.Constant(mesh, zero) * ufl.inner(p, q) * ufl.dx
f[1] = ufl.inner(dolfinx.fem.Constant(mesh, zero), v) * ufl.dx
bdofs_Y_1234 = dolfinx.fem.locate_dofs_topological(Y, mesh.topology.dim - 1, boundaries_1234)
bdofs_Q_1234 = dolfinx.fem.locate_dofs_topological(Q, mesh.topology.dim - 1, boundaries_1234)
bc = [dolfinx.fem.dirichletbc(zero, bdofs_Y_1234, Y),
      dolfinx.fem.dirichletbc(zero, bdofs_Q_1234, 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) * ufl.dx + 0.5 * alpha * ufl.inner(u, u) * ufl.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_1a_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.055555555)

Uncontrolled J = 0.055555555555555344

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_1a_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.0002337096)

Optimal J = 0.00023370965741544301

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

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

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