Hi everyone,
I’m fairly new to dolfinx and currently working on an adjoint-based optimization problem where I try to estimate an unknown inflow boundary condition.
I use a mixed function space on a 2D domain:
elem_u = element("Lagrange", domain.basix_cell(), 2, shape=(2,))
elem_p = element("Lagrange", domain.basix_cell(), 1)
W = fem.functionspace(domain, mixed_element([elem_u, elem_p]))
W0, W1 = W.sub(0), W.sub(1)
V, V_to_W = W0.collapse()
For the initial inflow condition, I set the x-velocity to 1.0:
inflow_facets = facet_tags.find(inflow_id)
inflow_dofs_W0 = fem.locate_dofs_topological((W0, V), domain.topology.dim - 1, inflow_facets)
u_inflow = fem.Function(V)
u_inflow.x.array[::2] = 1.0
bc_inflow = fem.dirichletbc(u_inflow, inflow_dofs_W0, W0)
This works fine for the forward solve.
However, inside my optimization loop, I want to update the inflow boundary based on a gradient computed in the same space V:
gradJ = -lam.x.array + gamma * u_inflow.x.array
u_inflow.x.array[inflow_dofs_V] -= alpha * gradJ[inflow_dofs_V]
where
inflow_dofs_V = fem.locate_dofs_topological(V, domain.topology.dim - 1, inflow_facets)
The problem is that inflow_dofs_V does not correspond to the DOFs actually constrained by the Dirichlet BC defined via (W0, V).
So my BC updates have no effect during optimization.
How can I correctly identify the velocity DOFs on the inflow boundary (in the collapsed space V) that correspond to the BC DOFs used in the mixed space W0, so that I can safely update them during optimization?
To illustrate the issue, I provide the MWE below that shows that accessing the boundary with ...[inflow_dofs_V] is not working.
Thanks in advance for your help!
import numpy as np
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from basix.ufl import element, mixed_element
import pyvista as pv
# --- Helper to visualize velocity BCs ---
def visualize_bc(w, scale=0.1):
u = w.sub(0).collapse()
grid = pv.UnstructuredGrid(*plot.vtk_mesh(u.function_space))
vals = np.zeros((len(u.x.array)//2, 3))
vals[:, :2] = u.x.array.reshape(-1, 2)
grid["u"] = vals
pl = pv.Plotter()
pl.add_mesh(grid, color="lightgray", opacity=0.4)
pl.add_mesh(grid.glyph(orient="u", factor=scale))
pl.view_xy()
pl.show()
# --- Mesh and spaces ---
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 5)
elem_u = element("Lagrange", domain.basix_cell(), 2, shape=(2,))
elem_p = element("Lagrange", domain.basix_cell(), 1)
W = fem.functionspace(domain, mixed_element([elem_u, elem_p]))
W0 = W.sub(0)
V, V_to_W = W0.collapse()
# --- Inflow boundary (x = 0) ---
facets = mesh.locate_entities_boundary(domain, 1, lambda x: np.isclose(x[0], 0))
inflow_dofs_W0 = fem.locate_dofs_topological((W0, V), 1, facets)
inflow_dofs_V = fem.locate_dofs_topological(V, 1, facets)
# --- Works: BC defined via (W0, V) ---
u1 = fem.Function(V)
u1.x.array[::2] = 1.0
bc1 = fem.dirichletbc(u1, inflow_dofs_W0, W0)
tmp = fem.Function(W)
bc1.set(tmp.x.array)
tmp.x.scatter_forward()
visualize_bc(tmp) # inflow visible on left boundary
# --- Fails: using inflow_dofs_V ---
u2 = fem.Function(V)
u2.x.array[inflow_dofs_V] = 1.0
bc2 = fem.dirichletbc(u2, inflow_dofs_W0, W0)
tmp2 = fem.Function(W)
bc2.set(tmp2.x.array)
tmp2.x.scatter_forward()
visualize_bc(tmp2) # no inflow visible
(dolfinx 0.9.0)