This tutorial is part of multiphenicsx.

Tutorial 06, case 6: Navier-Stokes problem with distributed control

In this tutorial we solve the optimal control problem

$$\min J(y, u) = \frac{1}{2} \int_{\Omega} |v - v_d|^2 dx + \frac{\alpha}{2} \int_{\Omega} |u|^2 dx$$ s.t. $$\begin{cases} - \nu \Delta v + v \cdot \nabla v + \nabla p = f + u & \text{in } \Omega\\ \text{div} v = 0 & \text{in } \Omega\\ v = 0 & \text{on } \partial\Omega \end{cases}$$

where $$\begin{align*} & \Omega & \text{unit square}\\ & u \in [L^2(\Omega)]^2 & \text{control variable}\\ & v \in [H^1_0(\Omega)]^2 & \text{state velocity variable}\\ & p \in L^2(\Omega) & \text{state pressure variable}\\ & \alpha > 0 & \text{penalization parameter}\\ & v_d & \text{desired state}\\ & \nu & \text{kinematic viscosity}\\ & 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 sympy
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_velocity = dolfinx.fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim, )))
Y_pressure = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
U = dolfinx.fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim, )))
Q_velocity = Y_velocity.clone()
Q_pressure = Y_pressure.clone()

Trial and test functions

(dv, dp) = (ufl.TrialFunction(Y_velocity), ufl.TrialFunction(Y_pressure))
(w, q) = (ufl.TestFunction(Y_velocity), ufl.TestFunction(Y_pressure))
du = ufl.TrialFunction(U)
r = ufl.TestFunction(U)
(dz, db) = (ufl.TrialFunction(Q_velocity), ufl.TrialFunction(Q_pressure))
(s, d) = (ufl.TestFunction(Q_velocity), ufl.TestFunction(Q_pressure))

Solution

(v, p) = (dolfinx.fem.Function(Y_velocity), dolfinx.fem.Function(Y_pressure))
u = dolfinx.fem.Function(U)
(z, b) = (dolfinx.fem.Function(Q_velocity), dolfinx.fem.Function(Q_pressure))

Problem data

alpha = 1.e-5
epsilon = 1.e-5
x, y = sympy.symbols("x[0], x[1]")
psi_d = 10 * (1 - sympy.cos(0.8 * np.pi * x)) * (1 - sympy.cos(0.8 * np.pi * y)) * (1 - x)**2 * (1 - y)**2
v_d_x = sympy.lambdify([x, y], psi_d.diff(y, 1))
v_d_y = sympy.lambdify([x, y], - psi_d.diff(x, 1))
v_d = dolfinx.fem.Function(Y_velocity)
v_d.interpolate(lambda x: np.stack((v_d_x(x[0], x[1]), v_d_y(x[0], x[1])), axis=0))
nu = 0.01
zero_scalar =  petsc4py.PETSc.ScalarType(0)  # type: ignore[attr-defined]
zero_vector = np.zeros((2, ), dtype=petsc4py.PETSc.ScalarType)  # type: ignore[attr-defined]
ff = dolfinx.fem.Constant(mesh, zero_vector)

Optimality conditions

F = [nu * ufl.inner(ufl.grad(z), ufl.grad(w)) * ufl.dx
     + ufl.inner(z, ufl.grad(w) * v) * ufl.dx + ufl.inner(z, ufl.grad(v) * w) * ufl.dx
     - ufl.inner(b, ufl.div(w)) * ufl.dx + ufl.inner(v - v_d, w) * ufl.dx,
     - ufl.inner(ufl.div(z), q) * ufl.dx + epsilon * ufl.inner(b, q) * ufl.dx,
     alpha * ufl.inner(u, r) * ufl.dx - ufl.inner(z, r) * ufl.dx,
     nu * ufl.inner(ufl.grad(v), ufl.grad(s)) * ufl.dx + ufl.inner(ufl.grad(v) * v, s) * ufl.dx
     - ufl.inner(p, ufl.div(s)) * ufl.dx - ufl.inner(u + ff, s) * ufl.dx,
     - ufl.inner(ufl.div(v), d) * ufl.dx + epsilon * ufl.inner(p, d) * ufl.dx]
dF = [[ufl.derivative(F_i, u_j, du_j) for (u_j, du_j) in zip((v, p, u, z, b), (dv, dp, du, dz, db))] for F_i in F]
dF[3][3] = dolfinx.fem.Constant(mesh, zero_scalar) * ufl.inner(dz, s) * ufl.dx
bdofs_Y_velocity_1234 = dolfinx.fem.locate_dofs_topological(
    Y_velocity, mesh.topology.dim - 1, boundaries_1234)
bdofs_Q_velocity_1234 = dolfinx.fem.locate_dofs_topological(
    Q_velocity, mesh.topology.dim - 1, boundaries_1234)
bc = [dolfinx.fem.dirichletbc(zero_vector, bdofs_Y_velocity_1234, Y_velocity),
      dolfinx.fem.dirichletbc(zero_vector, bdofs_Q_velocity_1234, Q_velocity)]

Cost functional

J = 0.5 * ufl.inner(v - v_d, v - v_d) * ufl.dx + 0.5 * alpha * ufl.inner(u, u) * ufl.dx
J_cpp = dolfinx.fem.form(J)

Uncontrolled functional value

# Create problem by extracting state forms from the optimality conditions
F_state = [ufl.replace(F[i], {s: w, d: q}) for i in (3, 4)]
dF_state = [[ufl.replace(dF[i][j], {s: w, d: q}) for j in (0, 1)] for i in (3, 4)]
bc_state = [bc[0]]

petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "ksp_error_if_not_converged": True,
    "snes_monitor": None,
    "snes_error_if_not_converged": True
}
problem_state = dolfinx.fem.petsc.NonlinearProblem(
    F_state, (v, p), bcs=bc_state, J=dF_state,
    petsc_options_prefix="tutorial_6_navier_stokes_state_", petsc_options=petsc_options,
    kind="mpi"
)
problem_state.solve()
del problem_state

  0 SNES Function norm 0.000000000000e+00

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

Uncontrolled J = 0.17845361493420545

viskex.dolfinx.plot_vector_field(v, "uncontrolled state velocity", glyph_factor=1e-1)

viskex.dolfinx.plot_scalar_field(p, "uncontrolled state pressure")

Optimal control

# PYTEST_XFAIL_AND_SKIP_NEXT: ufl.derivative(adjoint of the trilinear term) introduces spurious conj(trial)
"""Expect this cell to fail.

assert not np.issubdtype(petsc4py.PETSc.ScalarType, np.complexfloating)  # type: ignore[attr-defined]
"""

'Expect this cell to fail.\n\nassert not np.issubdtype(petsc4py.PETSc.ScalarType, np.complexfloating)  # type: ignore[attr-defined]\n'

"""Skip cell due to a previously xfailed cell.

problem = dolfinx.fem.petsc.NonlinearProblem(
    F, (v, p, u, z, b), bcs=bc, J=dF,
    petsc_options_prefix="tutorial_6_navier_stokes_", petsc_options=petsc_options,
    kind="mpi"
)
problem.solve()
del problem
"""

'Skip cell due to a previously xfailed cell.\n\nproblem = dolfinx.fem.petsc.NonlinearProblem(\n    F, (v, p, u, z, b), bcs=bc, J=dF,\n    petsc_options_prefix="tutorial_6_navier_stokes_", petsc_options=petsc_options,\n    kind="mpi"\n)\nproblem.solve()\ndel problem\n'

"""Skip cell due to a previously xfailed cell.

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, 9.9127635e-7)
"""

'Skip cell due to a previously xfailed cell.\n\nJ_controlled = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(J_cpp), op=mpi4py.MPI.SUM)\nprint("Optimal J =", J_controlled)\nassert np.isclose(J_controlled, 9.9127635e-7)\n'

"""Skip cell due to a previously xfailed cell.

viskex.dolfinx.plot_vector_field(v, "state velocity", glyph_factor=1e-1)
"""

'Skip cell due to a previously xfailed cell.\n\nviskex.dolfinx.plot_vector_field(v, "state velocity", glyph_factor=1e-1)\n'

"""Skip cell due to a previously xfailed cell.

viskex.dolfinx.plot_scalar_field(p, "state pressure")
"""

'Skip cell due to a previously xfailed cell.\n\nviskex.dolfinx.plot_scalar_field(p, "state pressure")\n'

"""Skip cell due to a previously xfailed cell.

viskex.dolfinx.plot_vector_field(u, "control", glyph_factor=1e-1)
"""

'Skip cell due to a previously xfailed cell.\n\nviskex.dolfinx.plot_vector_field(u, "control", glyph_factor=1e-1)\n'

"""Skip cell due to a previously xfailed cell.

viskex.dolfinx.plot_vector_field(z, "adjoint velocity", glyph_factor=1e4)
"""

'Skip cell due to a previously xfailed cell.\n\nviskex.dolfinx.plot_vector_field(z, "adjoint velocity", glyph_factor=1e4)\n'

"""Skip cell due to a previously xfailed cell.

viskex.dolfinx.plot_scalar_field(b, "adjoint pressure")
"""

'Skip cell due to a previously xfailed cell.\n\nviskex.dolfinx.plot_scalar_field(b, "adjoint pressure")\n'