This tutorial is part of multiphenicsx.

Tutorial 06, case 7a: Stokes problem with Dirichlet control

In this tutorial we solve the optimal control problem

$$\min J(y, u) = \frac{1}{2} \int_{\Omega_{obs}} |\text{curl} v|^2 dx + \frac{\alpha}{2} \int_{\Gamma_C} |\nabla_{\mathbf{t}} u|^2 ds$$ s.t. $$\begin{cases} - \nu \Delta v + \nabla p = f & \text{in } \Omega\\ \text{div} v = 0 & \text{in } \Omega\\ v = g & \text{on } \Gamma_{in}\\ v = 0 & \text{on } \Gamma_{w}\\ v \cdot \mathbf{n} = u & \text{on } \Gamma_{C}\\ v \cdot \mathbf{t} = 0 & \text{on } \Gamma_{C}\\ v \cdot \mathbf{n} = 0 & \text{on } \Gamma_{s}\\ \nu \partial_n v \cdot \mathbf{t} = 0 & \text{on } \Gamma_{s}\\ p n - \nu \partial_n v = 0 & \text{on } \Gamma_{N} \end{cases}$$

where $$\begin{align*} & \Omega & \text{unit square}\\ & \Gamma_{in} & \text{has boundary id 1}\\ & \Gamma_{s} & \text{has boundary id 2}\\ & \Gamma_{N} & \text{has boundary id 3}\\ & \Gamma_{C} & \text{has boundary id 4}\\ & \Gamma_{w} & \text{has boundary id 5}\\ & u \in L^2(\Gamma_C) & \text{control variable}\\ & v \in [H^1(\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}\\ & f & \text{forcing term}\\ & g & \text{inlet profile}\\ \end{align*}$$ using an adjoint formulation solved by a one shot approach.

The test case is from section 5 of

F. Negri, A. Manzoni and G. Rozza. Reduced basis approximation of parametrized optimal flow control problems for the Stokes equations. Computer and Mathematics with Applications, 69(4):319-336, 2015.

import dolfinx.fem
import dolfinx.io
import dolfinx.mesh
import gmsh
import mpi4py.MPI
import numpy as np
import numpy.typing
import petsc4py.PETSc
import ufl
import viskex

import multiphenicsx.fem
import multiphenicsx.fem.petsc

Geometrical parameters

L1 = 0.9
L2 = 0.35
L3 = 0.55
L4 = 0.2
H = 1.0
r = 0.1
mesh_size = 0.025

Mesh

gmsh.initialize()
gmsh.model.add("mesh")

p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, mesh_size)
p1 = gmsh.model.geo.addPoint(L1, 0.0, 0.0, mesh_size)
p2 = gmsh.model.geo.addPoint(L1 + L2, 0.0, 0.0, mesh_size)
p3 = gmsh.model.geo.addPoint(L1 + L2 + L3, 0.0, 0.0, mesh_size)
p4 = gmsh.model.geo.addPoint(L1 + L2 + L3 + L4, 0.0, 0.0, mesh_size)
p5 = gmsh.model.geo.addPoint(L1 + L2 + L3 + L4, H, 0.0, mesh_size)
p6 = gmsh.model.geo.addPoint(L1 + L2 + L3, H, 0.0, mesh_size)
p7 = gmsh.model.geo.addPoint(L1 + L2, H, 0.0, mesh_size)
p8 = gmsh.model.geo.addPoint(L1, H, 0.0, mesh_size)
p9 = gmsh.model.geo.addPoint(0.0, H, 0.0, mesh_size)
p10 = gmsh.model.geo.addPoint(L1, H / 2, 0.0, mesh_size)
p11 = gmsh.model.geo.addPoint(L1, H / 2 + r, 0.0, mesh_size)
p12 = gmsh.model.geo.addPoint(L1, H / 2 - r, 0.0, mesh_size)
p13 = gmsh.model.geo.addPoint(L1 + L2, H / 2 - r, 0.0, mesh_size)
p14 = gmsh.model.geo.addPoint(L1 + L2 + L3, H / 2 - 3 * r, 0.0, mesh_size)
p15 = gmsh.model.geo.addPoint(L1 + L2 + L3, H / 2 + 3 * r, 0.0, mesh_size)
p16 = gmsh.model.geo.addPoint(L1 + L2, H / 2 + r, 0.0, mesh_size)
l0 = gmsh.model.geo.addLine(p0, p1)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p5)
l5 = gmsh.model.geo.addLine(p5, p6)
l6 = gmsh.model.geo.addLine(p6, p7)
l7 = gmsh.model.geo.addLine(p7, p8)
l8 = gmsh.model.geo.addLine(p8, p9)
l9 = gmsh.model.geo.addLine(p9, p0)
l10 = gmsh.model.geo.addLine(p12, p13)
l11 = gmsh.model.geo.addLine(p13, p14)
l12 = gmsh.model.geo.addLine(p14, p15)
l13 = gmsh.model.geo.addLine(p15, p16)
l14 = gmsh.model.geo.addLine(p16, p11)
l15 = gmsh.model.geo.addLine(p13, p16)
l16 = gmsh.model.geo.addLine(p1, p12)
l17 = gmsh.model.geo.addLine(p11, p8)
l18 = gmsh.model.geo.addLine(p2, p13)
l19 = gmsh.model.geo.addLine(p16, p7)
l20 = gmsh.model.geo.addLine(p3, p14)
l21 = gmsh.model.geo.addLine(p15, p6)
c0 = gmsh.model.geo.addCircleArc(p11, p10, p12)
line_loop_subdomain1 = gmsh.model.geo.addCurveLoop([l0, l16, -c0, l17, l8, l9])
line_loop_subdomain2a = gmsh.model.geo.addCurveLoop([l1, l18, -l10, -l16])
line_loop_subdomain2b = gmsh.model.geo.addCurveLoop([l7, -l17, -l14, l19])
line_loop_subdomain3a = gmsh.model.geo.addCurveLoop([l2, l20, -l11, -l18])
line_loop_subdomain3b = gmsh.model.geo.addCurveLoop([l6, -l19, -l13, l21])
line_loop_subdomain3c = gmsh.model.geo.addCurveLoop([l3, l4, l5, -l21, -l12, -l20])
line_loop_subdomain4 = gmsh.model.geo.addCurveLoop([l11, l12, l13, -l15])
subdomain1 = gmsh.model.geo.addPlaneSurface([line_loop_subdomain1])
subdomain2a = gmsh.model.geo.addPlaneSurface([line_loop_subdomain2a])
subdomain2b = gmsh.model.geo.addPlaneSurface([line_loop_subdomain2b])
subdomain3a = gmsh.model.geo.addPlaneSurface([line_loop_subdomain3a])
subdomain3b = gmsh.model.geo.addPlaneSurface([line_loop_subdomain3b])
subdomain3c = gmsh.model.geo.addPlaneSurface([line_loop_subdomain3c])
subdomain4 = gmsh.model.geo.addPlaneSurface([line_loop_subdomain4])

gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [l9], 1)
gmsh.model.addPhysicalGroup(1, [l0, l1, l2, l3, l5, l6, l7, l8], 2)
gmsh.model.addPhysicalGroup(1, [l4], 3)
gmsh.model.addPhysicalGroup(1, [l10, l14], 4)
gmsh.model.addPhysicalGroup(1, [c0, l15], 5)
gmsh.model.addPhysicalGroup(2, [subdomain1], 1)
gmsh.model.addPhysicalGroup(2, [subdomain2a, subdomain2b], 2)
gmsh.model.addPhysicalGroup(2, [subdomain3a, subdomain3b, subdomain3c], 3)
gmsh.model.addPhysicalGroup(2, [subdomain4], 4)
gmsh.model.mesh.generate(2)

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 10%] Meshing curve 2 (Line)
Info    : [ 10%] Meshing curve 3 (Line)
Info    : [ 20%] Meshing curve 4 (Line)
Info    : [ 20%] Meshing curve 5 (Line)
Info    : [ 30%] Meshing curve 6 (Line)
Info    : [ 30%] Meshing curve 7 (Line)
Info    : [ 40%] Meshing curve 8 (Line)
Info    : [ 40%] Meshing curve 9 (Line)
Info    : [ 40%] Meshing curve 10 (Line)
Info    : [ 50%] Meshing curve 11 (Line)
Info    : [ 50%] Meshing curve 12 (Line)
Info    : [ 60%] Meshing curve 13 (Line)
Info    : [ 60%] Meshing curve 14 (Line)
Info    : [ 70%] Meshing curve 15 (Line)
Info    : [ 70%] Meshing curve 16 (Line)
Info    : [ 70%] Meshing curve 17 (Line)
Info    : [ 80%] Meshing curve 18 (Line)
Info    : [ 80%] Meshing curve 19 (Line)
Info    : [ 90%] Meshing curve 20 (Line)
Info    : [ 90%] Meshing curve 21 (Line)
Info    : [100%] Meshing curve 22 (Line)
Info    : [100%] Meshing curve 23 (Circle)
Info    : Done meshing 1D (Wall 0.00266516s, CPU 0.003343s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : [ 20%] Meshing surface 2 (Plane, Frontal-Delaunay)
Info    : [ 30%] Meshing surface 3 (Plane, Frontal-Delaunay)
Info    : [ 50%] Meshing surface 4 (Plane, Frontal-Delaunay)
Info    : [ 60%] Meshing surface 5 (Plane, Frontal-Delaunay)
Info    : [ 80%] Meshing surface 6 (Plane, Frontal-Delaunay)
Info    : [ 90%] Meshing surface 7 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.0777169s, CPU 0.075036s)
Info    : 3824 nodes 7815 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)

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

# Normal and tangent
n = ufl.FacetNormal(mesh)
t = ufl.as_vector([n[1], -n[0]])

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_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))
L = U.clone()
Q_velocity = Y_velocity.clone()
Q_pressure = Y_pressure.clone()

Y_velocity_0 = Y_velocity.sub(0)
Y_velocity_1 = Y_velocity.sub(1)
Q_velocity_0 = Q_velocity.sub(0)
Q_velocity_1 = Q_velocity.sub(1)

Restrictions

dofs_Y_velocity = np.arange(0, Y_velocity.dofmap.index_map.size_local + Y_velocity.dofmap.index_map.num_ghosts)
dofs_Y_pressure = np.arange(0, Y_pressure.dofmap.index_map.size_local + Y_pressure.dofmap.index_map.num_ghosts)
dofs_U = dolfinx.fem.locate_dofs_topological(U, boundaries.dim, boundaries.indices[boundaries.values == 4])
dofs_L = dofs_U
dofs_Q_velocity = dofs_Y_velocity
dofs_Q_pressure = dofs_Y_pressure
restriction_Y_velocity = multiphenicsx.fem.DofMapRestriction(Y_velocity.dofmap, dofs_Y_velocity)
restriction_Y_pressure = multiphenicsx.fem.DofMapRestriction(Y_pressure.dofmap, dofs_Y_pressure)
restriction_U = multiphenicsx.fem.DofMapRestriction(U.dofmap, dofs_U)
restriction_L = multiphenicsx.fem.DofMapRestriction(L.dofmap, dofs_L)
restriction_Q_velocity = multiphenicsx.fem.DofMapRestriction(Q_velocity.dofmap, dofs_Q_velocity)
restriction_Q_pressure = multiphenicsx.fem.DofMapRestriction(Q_pressure.dofmap, dofs_Q_pressure)
restriction = [
    restriction_Y_velocity, restriction_Y_pressure, restriction_U, restriction_L,
    restriction_Q_velocity, restriction_Q_pressure]

Trial and test functions

(v, p) = (ufl.TrialFunction(Y_velocity), ufl.TrialFunction(Y_pressure))
(w, q) = (ufl.TestFunction(Y_velocity), ufl.TestFunction(Y_pressure))
(u, l) = (ufl.TrialFunction(U), ufl.TrialFunction(L))
(r, m) = (ufl.TestFunction(U), ufl.TestFunction(L))
(z, b) = (ufl.TrialFunction(Q_velocity), ufl.TrialFunction(Q_pressure))
(s, d) = (ufl.TestFunction(Q_velocity), ufl.TestFunction(Q_pressure))

Problem data

def non_zero_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
        petsc4py.PETSc.ScalarType]:
    """Return the flat velocity profile at the inlet."""
    values = np.zeros((2, x.shape[1]))
    values[0, :] = 2.5
    return values


nu = 1.
alpha = 1.e-2
ff = dolfinx.fem.Constant(mesh, tuple(petsc4py.PETSc.ScalarType(0) for _ in range(2)))
bc0 = dolfinx.fem.Function(Y_velocity)
bc0_component = dolfinx.fem.Function(Y_velocity_0.collapse()[0])
bc1 = dolfinx.fem.Function(Y_velocity)
bc1.interpolate(non_zero_eval)

Optimality conditions

def vorticity(v: ufl.Argument, w: ufl.Argument) -> ufl.core.expr.Expr:  # type: ignore[no-any-unimported]
    """Return the UFL expression corresponding to the inner(curl, curl) operator."""
    return ufl.inner(ufl.curl(v), ufl.curl(w))


def penalty(u: ufl.Argument, r: ufl.Argument) -> ufl.core.expr.Expr:  # type: ignore[no-any-unimported]
    """Return the UFL expression corresponding to the penalty term."""
    return alpha * ufl.inner(ufl.dot(ufl.grad(u), t), ufl.dot(ufl.grad(r), t))


a = [[vorticity(v, w) * dx(4), None, None, ufl.inner(l, ufl.dot(w, n)) * ds(4),
      nu * ufl.inner(ufl.grad(z), ufl.grad(w)) * dx, - ufl.inner(b, ufl.div(w)) * dx],
     [None, None, None, None, - ufl.inner(ufl.div(z), q) * dx, None],
     [None, None, penalty(u, r) * ds(4), - ufl.inner(l, r) * ds(4), None, None],
     [ufl.inner(ufl.dot(v, n), m) * ds(4), None, - ufl.inner(u, m) * ds(4), None, None, None],
     [nu * ufl.inner(ufl.grad(v), ufl.grad(s)) * dx, - ufl.inner(p, ufl.div(s)) * dx, None, None, None, None],
     [- ufl.inner(ufl.div(v), d) * dx, None, None, None, None, None]]
f = [None,
     None,
     None,
     None,
     ufl.inner(ff, s) * dx,
     None]
a[0][0] += dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)) * ufl.inner(v, w) * dx
a[4][4] = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)) * ufl.inner(z, s) * dx
f[0] = ufl.inner(dolfinx.fem.Constant(mesh, tuple(petsc4py.PETSc.ScalarType(0) for _ in range(2))), w) * dx
f[1] = ufl.inner(dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)), q) * dx
f[2] = ufl.inner(dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)), r) * dx
f[3] = ufl.inner(dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)), m) * dx
f[5] = ufl.inner(dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)), d) * dx
a_cpp = dolfinx.fem.form(a)
f_cpp = dolfinx.fem.form(f)


def bdofs(
    space_from: dolfinx.fem.FunctionSpace, space_to: dolfinx.fem.FunctionSpace, idx: int
) -> np.typing.NDArray[np.int32]:
    """Locate DOFs on the boundary `idx`."""
    return dolfinx.fem.locate_dofs_topological(
        (space_from, space_to), mesh.topology.dim - 1, boundaries.indices[boundaries.values == idx])


bc = [
    dolfinx.fem.dirichletbc(
        bc1, bdofs(Y_velocity, bc1.function_space, 1), Y_velocity),
    dolfinx.fem.dirichletbc(
        bc0_component, bdofs(Y_velocity_1, bc0_component.function_space, 2), Y_velocity_1),
    dolfinx.fem.dirichletbc(
        bc0_component, bdofs(Y_velocity_0, bc0_component.function_space, 4), Y_velocity_0),
    dolfinx.fem.dirichletbc(
        bc0, bdofs(Y_velocity, bc0.function_space, 5), Y_velocity),
    dolfinx.fem.dirichletbc(
        bc0, bdofs(Q_velocity, bc0.function_space, 1), Q_velocity),
    dolfinx.fem.dirichletbc(
        bc0_component, bdofs(Q_velocity_1, bc0_component.function_space, 2), Q_velocity_1),
    dolfinx.fem.dirichletbc(
        bc0, bdofs(Q_velocity, bc0.function_space, 4), Q_velocity),
    dolfinx.fem.dirichletbc(
        bc0, bdofs(Q_velocity, bc0.function_space, 5), Q_velocity)
]

Solution

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

Cost functional

J = 0.5 * vorticity(v, v) * dx(4) + 0.5 * penalty(u, u) * ds(4)
J_cpp = dolfinx.fem.form(J)

Uncontrolled functional value

# Extract state forms from the optimality conditions
a_state = [[ufl.replace(a[i][j], {s: w, d: q}) if a[i][j] is not None else None for j in (0, 1)] for i in (4, 5)]
f_state = [ufl.replace(f[i], {s: w, d: q}) for i in (4, 5)]
a_state_cpp = dolfinx.fem.form(a_state)
f_state_cpp = dolfinx.fem.form(f_state)
bc_state = [
    dolfinx.fem.dirichletbc(
        bc1, bdofs(Y_velocity, bc1.function_space, 1), Y_velocity),
    dolfinx.fem.dirichletbc(
        bc0_component, bdofs(Y_velocity_1, bc0_component.function_space, 2), Y_velocity_1),
    dolfinx.fem.dirichletbc(
        bc0, bdofs(Y_velocity, bc0.function_space, 4), Y_velocity),
    dolfinx.fem.dirichletbc(
        bc0, bdofs(Y_velocity, bc0.function_space, 5), Y_velocity)
]

# Assemble the block linear system for the state
A_state = multiphenicsx.fem.petsc.assemble_matrix_block(
    a_state_cpp, bcs=bc_state, restriction=([restriction[i] for i in (4, 5)], [restriction[j] for j in (0, 1)]))
A_state.assemble()
F_state = multiphenicsx.fem.petsc.assemble_vector_block(
    f_state_cpp, a_state_cpp, bcs=bc_state, restriction=[restriction[i] for i in (4, 5)])

# Solve
vp = multiphenicsx.fem.petsc.create_vector_block(
    [f_cpp[j] for j in (0, 1)], restriction=[restriction[j] for j in (0, 1)])
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, vp)
vp.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
ksp.destroy()

<petsc4py.PETSc.KSP at 0x7f2dc20d0130>

# Split the block solution in components
with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(
        vp, [c.function_space.dofmap for c in (v, p)], [restriction[j] for j in (0, 1)]) as vp_wrapper:
    for vp_wrapper_local, component in zip(vp_wrapper, (v, p)):
        with component.x.petsc_vec.localForm() as component_local:
            component_local[:] = vp_wrapper_local
vp.destroy()

<petsc4py.PETSc.Vec at 0x7f2d501f9fd0>

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

Uncontrolled J = 2.923619390750833

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

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, "uncontrolled state pressure")

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
vpulzb = 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, vpulzb)
vpulzb.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
ksp.destroy()

<petsc4py.PETSc.KSP at 0x7f2dcd6d2700>

# Split the block solution in components
with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(
        vpulzb, [c.function_space.dofmap for c in (v, p, u, l, z, b)], restriction) as vpulzb_wrapper:
    for vpulzb_wrapper_local, component in zip(vpulzb_wrapper, (v, p, u, l, z, b)):
        with component.x.petsc_vec.localForm() as component_local:
            component_local[:] = vpulzb_wrapper_local
vpulzb.destroy()

<petsc4py.PETSc.Vec at 0x7f2d5026c9f0>

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

Optimal J = 1.710273971369316

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

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, "state pressure")

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(l, "lambda")

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_vector_field(z, "adjoint velocity", glyph_factor=1)

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(b, "adjoint pressure")

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