Dirichlet BC control of the Stokes equations#
This demo is based on the example stokes-bc-control.py from the dolfin-adjoint source tree,
authored by Simon W. Funke and André Massing.
The updated version has been made by Jørgen S. Dokken dokken@simula.no.
This example demonstrates how to compute the sensitivity with respect to the Dirichlet boundary conditions in pyadjoint.
Problem definition#
Consider the problem of minimising the compliance
subject to the Stokes equations
with Dirichlet boundary conditions
where :math:\Omega is the domain of interest,
:math:u:\Omega \to \mathbb R^2 is the unknown velocity,
:math:p:\Omega \to \mathbb R is the unknown pressure, :math:\nu
is the viscosity, :math:\alpha is the regularisation parameter,
:math:f denotes the value for the Dirichlet inflow boundary
condition, and :math:g is the control variable that specifies the
Dirichlet boundary condition on the circle.
Physically, this setup corresponds to minimising the loss of flow energy into heat by actively controlling the in/outflow at the circle boundary. To avoid excessive control solutions, non-zero control values are penalised via the regularisation term.
Implementation#
First, we implement the various modules we will use in this demo
import itertools
from dolfinx_adjoint import Constant, Function, LinearProblem, assemble_scalar, dirichletbc
try:
from dolfinx.io import gmsh as gmshio
except ImportError:
from dolfinx.io import gmshio # type: ignore[attr-defined, no-redef]
The geometry is the one the original demo meshes with mshr in its own make-mesh.py:
a 30x10 rectangle with a circle of radius 2.5 centred at (10, 5) removed. Keeping those
dimensions matters, because the inflow profile below is written for them.
The gmsh commands themselves follow DOLFINx-tutorial Navier-Stokes benchmark, which explains the mesh generation process in more detail. The mesh is graded towards the circle, since that is where the control acts and the flow varies most.
gmsh.initialize()
L = 30.0
H = 10.0
c_x, c_y = 10.0, 5.0
r = 2.5
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
fluid_marker = 1
res_min = r / 3
if mesh_comm.rank == model_rank:
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
_obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)
fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, _obstacle)])
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=gdim)
assert len(volumes) == 1
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(center_of_mass, [0, H / 2, 0]):
inflow.append(boundary[1])
elif np.allclose(center_of_mass, [L, H / 2, 0]):
outflow.append(boundary[1])
elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]):
walls.append(boundary[1])
else:
obstacle.append(boundary[1])
gmsh.model.addPhysicalGroup(1, walls, wall_marker)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")
distance_field = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
threshold_field = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.25 * H)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * H)
min_field = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")
mesh_data = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
mesh = mesh_data.mesh
assert mesh_data.facet_tags is not None
ft = mesh_data.facet_tags
ft.name = "Facet markers"
Info : [ 0%] Difference
Info : [ 10%] Difference
Info : [ 20%] Difference
Info : [ 30%] Difference - Performing Face-Face intersection
Info : [ 70%] Difference - Performing intersection of shapes
Info : [ 80%] Difference - Making faces
Info : [ 90%] Difference - Adding holes
Info : Meshing 1D...
Info : [ 0%] Meshing curve 5 (Ellipse)
Info : [ 30%] Meshing curve 6 (Line)
Info : [ 50%] Meshing curve 7 (Line)
Info : [ 70%] Meshing curve 8 (Line)
Info : [ 90%] Meshing curve 9 (Line)
Info : Done meshing 1D (Wall 0.00604822s, CPU 0.006191s)
Info : Meshing 2D...
Info : Meshing surface 1 (Plane, Frontal-Delaunay for Quads)
Info : Simple recombination completed (Wall 0.00111728s, CPU 0.001118s): 61 quads, 9 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = 0.799046, min Q = 0.429862
Info : Simple recombination completed (Wall 0.00251103s, CPU 0.002511s): 271 quads, 0 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = 0.833, min Q = 0.348885
Info : Done meshing 2D (Wall 0.00562271s, CPU 0.005673s)
Info : Refining mesh...
Info : Meshing order 2 (curvilinear on)...
Info : [ 0%] Meshing curve 5 order 2
Info : [ 20%] Meshing curve 6 order 2
Info : [ 40%] Meshing curve 7 order 2
Info : [ 60%] Meshing curve 8 order 2
Info : [ 70%] Meshing curve 9 order 2
Info : [ 90%] Meshing surface 1 order 2
Info : Done meshing order 2 (Wall 0.00174981s, CPU 0.001785s)
Info : Done refining mesh (Wall 0.00197892s, CPU 0.002092s)
Info : 1170 nodes 1261 elements
Info : Meshing order 2 (curvilinear on)...
Info : [ 0%] Meshing curve 5 order 2
Info : [ 20%] Meshing curve 6 order 2
Info : [ 40%] Meshing curve 7 order 2
Info : [ 60%] Meshing curve 8 order 2
Info : [ 70%] Meshing curve 9 order 2
Info : [ 90%] Meshing surface 1 order 2
Info : Done meshing order 2 (Wall 0.0058518s, CPU 0.005961s)
Info : Optimizing mesh (Netgen)...
Info : Done optimizing mesh (Wall 1.333e-06s, CPU 2e-06s)
Visualising the mesh and the boundary markers#
Before solving anything, we look at what we have meshed. The first panel shows the graded
mesh; the second draws only the tagged facets, coloured by the marker each carries, over a
faint wireframe of the domain. The control acts on the Obstacle facets alone, so it is
worth confirming those are the ones that were tagged.
Defining the function spaces and boundary conditions#
Then, we define the discrete function spaces. A Taylor-Hood finite-element pair is a suitable choice for the Stokes equations. The control function is the Dirichlet boundary value on the velocity field and is hence be a function on the velocity space
el_u = basix.ufl.element("Lagrange", mesh.basix_cell(), 2, shape=(gdim,))
el_p = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)
V = dolfinx.fem.functionspace(mesh, el_u)
Q = dolfinx.fem.functionspace(mesh, el_p)
W = ufl.MixedFunctionSpace(V, Q)
u, p = ufl.TrialFunctions(W)
v, q = ufl.TestFunctions(W)
Our functional requires the computation of a boundary integral
over :math:\partial \Omega_{\mathrm{circle}}. Therefore, we need
to create a measure for this integral, which will be accessible as
:py:data:ds(2) in the definition of the functional. In addition, we
define our strong Dirichlet boundary conditions.
ds = ufl.ds(subdomain_data=ft)
# Define boundary conditions
x = ufl.SpatialCoordinate(mesh)
u_inflow = ufl.as_vector((x[1] * (10 - x[1]) / 25, 0))
noslip = Constant(mesh, (0, 0))
Locate the degrees of freedom for the various boundary facets.
noslip_dofs = dolfinx.fem.locate_dofs_topological(V, ft.dim, ft.find(wall_marker))
noslip = dirichletbc(noslip, noslip_dofs, V)
inflow_dofs = dolfinx.fem.locate_dofs_topological(V, ft.dim, ft.find(inlet_marker))
inflow = dirichletbc(u_inflow, inflow_dofs, V)
circle_dofs = dolfinx.fem.locate_dofs_topological(V, ft.dim, ft.find(obstacle_marker))
circle = dirichletbc(g, circle_dofs, V)
bcs = [inflow, noslip, circle]
We derive the standard weak formulation of
the Stokes problem: Find :math:u, p such that for all test
functions :math:v, q
with
In code, this becomes:
a = (
nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
- ufl.inner(p, ufl.div(v)) * ufl.dx
- ufl.inner(q, ufl.div(u)) * ufl.dx
)
# Named `L_forms`, not `L`: the channel length above is already called `L`, and rebinding
# it here would leave the plotting code below unable to refer to the geometry.
L_forms = [ufl.ZeroBaseForm((v,)), ufl.ZeroBaseForm((q,))]
Next we assemble and solve the system once to record it with
:py:mod:dolin-adjoint.
direct_solver_options = {
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"ksp_error_if_not_converged": True,
}
problem = LinearProblem(
ufl.extract_blocks(a),
L_forms,
u=[uh, ph],
bcs=bcs,
petsc_options=direct_solver_options,
adjoint_petsc_options=direct_solver_options,
tlm_petsc_options=direct_solver_options,
)
problem.solve()
[Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, quadrilateral, 2, equispaced, unset, False, float64, []), (2,)), 0), blocked element (Basix element (P, quadrilateral, 2, gll_warped, unset, False, float64, []), (2,))), 0),
Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, quadrilateral, 2, equispaced, unset, False, float64, []), (2,)), 0), Basix element (P, quadrilateral, 1, gll_warped, unset, False, float64, [])), 1)]
Next we define the functional of interest :math:J, the
optimisation parameter :math:g, and create the reduced
functional.
# The iteration monitor below wants the gradient norm as well as the functional value, but
# SciPy's callback is handed neither. `derivative_cb_post` fires every time the reduced
# functional computes a gradient -- which L-BFGS-B does once per iterate anyway -- so
# stashing the norm here costs no extra adjoint solve. Note pyadjoint assigns this callback's
# return value back over the derivatives, so it has to hand the list straight back.
gradient_norm = [float("nan")]
def record_gradient(_functional_value, derivatives, _control_values):
(gradient,) = derivatives
owned = gradient.x.index_map.size_local * gradient.x.block_size
local_max = np.abs(gradient.x.array[:owned]).max(initial=0.0)
gradient_norm[0] = mesh.comm.allreduce(local_max, op=MPI.MAX)
return derivatives
Jhat = pyadjoint.ReducedFunctional(J, m, derivative_cb_post=record_gradient)
Now, everything is set up to run the optimisation and to plot the
results. By default, :py:func:minimize uses the L-BFGS-B
algorithm.
# SciPy deprecated L-BFGS-B's own `disp`/`iprint` options in 1.15 and they print nothing as
# of 1.17 (removal is scheduled for 1.18), so `options={"disp": True}` is silently ignored.
# The supported replacement is a callback, which SciPy calls once per accepted iteration;
# {py:func}`pyadjoint.minimize` forwards any extra keyword straight to
# {py:func}`scipy.optimize.minimize`.
# Naming the parameter `intermediate_result` is what selects SciPy's newer callback API --
# the one handed the functional value, rather than only the control vector.
iteration = itertools.count()
def log_iteration(intermediate_result):
"""Log the iteration number, functional value, and gradient norm.
The gradient reported is the last one computed, which is the one at this iterate.
`max|dJ/dg|` is the unconstrained analogue of the `|proj g|` column the original
dolfin-adjoint demo printed.
"""
print(
f" iteration {next(iteration):3d} J = {intermediate_result.fun:.8e} max|dJ/dg| = {gradient_norm[0]:.3e}",
flush=True,
)
g_opt = pyadjoint.minimize(Jhat, method="L-BFGS-B", callback=log_iteration)
iteration 0 J = 3.25884382e+01 max|dJ/dg| = 2.487e+00
iteration 1 J = 2.04757012e+01 max|dJ/dg| = 3.806e-01
iteration 2 J = 2.01558026e+01 max|dJ/dg| = 2.138e-01
iteration 3 J = 1.99831092e+01 max|dJ/dg| = 3.914e-02
iteration 4 J = 1.99818393e+01 max|dJ/dg| = 2.817e-02
iteration 5 J = 1.99813771e+01 max|dJ/dg| = 2.801e-03
iteration 6 J = 1.99813637e+01 max|dJ/dg| = 1.480e-03
iteration 7 J = 1.99813603e+01 max|dJ/dg| = 8.727e-04
iteration 8 J = 1.99813600e+01 max|dJ/dg| = 5.614e-04
iteration 9 J = 1.99813598e+01 max|dJ/dg| = 1.778e-05
iteration 10 J = 1.99813598e+01 max|dJ/dg| = 1.154e-05
Results#
minimize returns the optimal control, but it leaves the tape at whatever point the line
search last probed. Re-evaluating the reduced functional at g_opt replays the forward
problem there, which both reports the optimal functional value and leaves uh/ph holding
the corresponding state – the velocity and pressure we want to look at.
J_opt = Jhat(g_opt)
print(f"J(g_opt) = {J_opt:.6e}")
J(g_opt) = 1.998136e+01
We visualise the three fields together: the optimised boundary control on the circle, and the velocity and pressure it induces.
The control and the velocity are vector fields, so we draw them as glyphs, following
the tutorial’s Navier-Stokes demo.
dolfinx.plot.vtk_mesh(V) places one point per degree of freedom of the (P2) velocity
space, so the glyph grid is built on the space rather than on the mesh.
The optimal control draws fluid into the circle on the upstream side and expels it downstream: rather than forcing the flow around the obstacle, it lets the obstacle pass flow through, which is what reduces the energy dissipated into heat. The strongest control sits on the upstream face, where the oncoming flow would otherwise stagnate.
The regularisation term alpha/2 * inner(g, g) * ds(obstacle_marker) is what keeps that
control finite; raising alpha shrinks it towards zero.
Verifying the implementation#
Finally we check the derivative that drove the optimisation, with a Taylor test: perturb the
control by eps * h and watch how fast the remainder of a Taylor expansion of the reduced
functional decays as eps shrinks.
Dropping the derivative term (dJdm=0) leaves a remainder of size \(\mathcal{O}(\epsilon)\),
so it should converge at rate 1 – that only confirms the functional is being re-evaluated
at the perturbed control. Including the adjoint gradient leaves
\(\mathcal{O}(\epsilon^2)\) and so rate 2, and that is the test of the adjoint: rate 2 is
reached only if Jhat.derivative() really is the gradient of Jhat.
The expansion point is the zero control, not g_opt. At the optimum the gradient vanishes,
which would take the first-order remainder with it and collapse the rate-1 check into the
rate-2 one, testing nothing. The direction is g_opt itself, which is supported exactly on
the circle – a direction supported elsewhere would leave the functional unchanged, since
the control enters the problem only through the boundary condition on that circle.
with pyadjoint.stop_annotating():
expansion_point = Function(V, name="TaylorExpansionPoint")
rate_0 = pyadjoint.taylor_test(Jhat, expansion_point, g_opt, dJdm=0)
rate_1 = pyadjoint.taylor_test(Jhat, expansion_point, g_opt)
print(f"Taylor convergence rates: zeroth order {rate_0:.3f} (expect 1), first order {rate_1:.3f} (expect 2)")
assert np.isclose(rate_0, 1.0, atol=0.1), f"zeroth-order Taylor rate {rate_0} is not 1"
assert np.isclose(rate_1, 2.0, atol=0.1), f"first-order Taylor rate {rate_1} is not 2"
Running Taylor test
Computed residuals: [0.545462411674194, 0.2734164601107949, 0.1368795436239978, 0.0684826002039216]
Computed convergence rates: [0.9963796843863362, 0.9981932433249101, 0.9990974694676765]
Running Taylor test
Computed residuals: [0.0027410170933611644, 0.0006852542729826605, 0.0001713135678909905, 4.282839202278399e-05]
Computed convergence rates: [2.000000000752935, 2.0000000029868463, 1.9999999983145007]
Taylor convergence rates: zeroth order 0.996 (expect 1), first order 2.000 (expect 2)