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} - \epsilon \Delta y + \beta \cdot \nabla y + \sigma 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}\\ & \epsilon > 0 & \text{diffusion coefficient}\\ & \beta \in \mathbb{R}^2 & \text{advection field}\\ & \sigma > 0 & \text{reaction coefficient}\\ & 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 = 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)
error: XDG_RUNTIME_DIR is invalid or not set in the environment. error: XDG_RUNTIME_DIR is invalid or not set in the environment. error: XDG_RUNTIME_DIR is invalid or not set in the environment.
viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")
error: XDG_RUNTIME_DIR is invalid or not set in the environment. error: XDG_RUNTIME_DIR is invalid or not set in the environment.
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 = 1.e-5
y_d = 1.
epsilon = 1.e-1
beta = ufl.as_vector((-1., -2.))
sigma = 1.
ff = 1.
bc0 = petsc4py.PETSc.ScalarType(0)
state_operator = (epsilon * ufl.inner(ufl.grad(y), ufl.grad(q)) * ufl.dx
+ ufl.inner(ufl.dot(beta, ufl.grad(y)), q) * ufl. dx + sigma * ufl.inner(y, q) * ufl.dx)
adjoint_operator = (epsilon * ufl.inner(ufl.grad(p), ufl.grad(z)) * ufl.dx
- ufl.inner(ufl.dot(beta, ufl.grad(p)), z) * ufl.dx + sigma * ufl.inner(p, z) * ufl.dx)
a = [[ufl.inner(y, z) * ufl.dx, None, adjoint_operator],
[None, alpha * ufl.inner(u, v) * ufl.dx, - ufl.inner(p, v) * ufl.dx],
[state_operator, - 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, petsc4py.PETSc.ScalarType(0)) * ufl.inner(p, q) * ufl.dx
f[1] = ufl.inner(dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)), 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(bc0, bdofs_Y_1234, Y),
dolfinx.fem.dirichletbc(bc0, bdofs_Q_1234, Q)]
(y, u, p) = (dolfinx.fem.Function(Y), dolfinx.fem.Function(U), dolfinx.fem.Function(Q))
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)
# 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_2a_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.37987224)
Uncontrolled J = 0.37987224632571426
viskex.dolfinx.plot_scalar_field(y, "uncontrolled state")
error: XDG_RUNTIME_DIR is invalid or not set in the environment. error: XDG_RUNTIME_DIR is invalid or not set in the environment.
problem = dolfinx.fem.petsc.LinearProblem(
a, f, bcs=bc, u=(y, u, p),
petsc_options_prefix="tutorial_2a_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.02899012)
Optimal J = 0.028990128663930536
viskex.dolfinx.plot_scalar_field(y, "state")
error: XDG_RUNTIME_DIR is invalid or not set in the environment. error: XDG_RUNTIME_DIR is invalid or not set in the environment.
viskex.dolfinx.plot_scalar_field(u, "control")
error: XDG_RUNTIME_DIR is invalid or not set in the environment. error: XDG_RUNTIME_DIR is invalid or not set in the environment.
viskex.dolfinx.plot_scalar_field(p, "adjoint")
error: XDG_RUNTIME_DIR is invalid or not set in the environment. error: XDG_RUNTIME_DIR is invalid or not set in the environment.