# Unable to assemble a system created after the derivative of mixed elements

Hi there,

the mwe below fails when I try to assemble a system composed by partial derivatives calculated using a state variable and a residual, both created with mixed elements:

from fenics import *
import numpy as np
import ufl as ufl
from mpi4py import MPI

mu_fluid = 0.0089
rho_fluid = 1.0
L = 1.0
u_inlet = 0.1
p_outlet = 0.0
Re = rho_fluid*u_inlet*L/mu_fluid
lx, ly, lz = 2.0, 1.0, 1.0
nelx, nely, nelz = 20, 10, 10
mesh = BoxMesh.create(MPI.COMM_WORLD, [Point(0.0, 0.0, 0.0), Point(lx, ly, lz)],
[nelx, nely, nelz],
CellType.Type.hexahedron)

V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
solution_space = MixedElement([V,Q])
W = FunctionSpace(mesh, solution_space)

bcs = []
fluid_boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
fluid_boundary_parts.set_all(3)
no_slip_boundary_right  = CompiledSubDomain("near(x[1],0.0)")
no_slip_boundary_left   = CompiledSubDomain("near(x[1],ly)", ly = ly)
no_slip_boundary_bottom = CompiledSubDomain("near(x[2],0.0)")
no_slip_boundary_top    = CompiledSubDomain("near(x[2],lz)", lz=lz)
inlet_front = CompiledSubDomain("on_boundary && near(x[0],0.0)")
outlet_back = CompiledSubDomain("on_boundary && near(x[0],lx)", lx=lx)
no_slip_boundary_right.mark(fluid_boundary_parts, 2)
no_slip_boundary_left.mark(fluid_boundary_parts, 2)
no_slip_boundary_bottom.mark(fluid_boundary_parts, 2)
no_slip_boundary_top.mark(fluid_boundary_parts, 2)
inlet_front.mark(fluid_boundary_parts, 0)
outlet_back.mark(fluid_boundary_parts, 1)
bcu_no_slip_right  = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_right)
bcu_no_slip_left   = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_left)
bcu_no_slip_bottom = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_bottom)
bcu_no_slip_top    = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_top)
inflow_profile    = Expression(('u_inlet','0','0'), u_inlet=u_inlet, degree=2)
bcu_inlet         = DirichletBC(W.sub(0), inflow_profile, inlet_front)
bcu_back_pressure = DirichletBC(W.sub(1), Constant(0), outlet_back)
bcs = [bcu_no_slip_right,  bcu_no_slip_left,\
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet,          bcu_back_pressure]

w      = Function(W)
dw     = TrialFunction(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)
F_navier_stokes = (
- inner(p,div(v))*dx
+ q*div(u)*dx
)
J_navier_stokes = derivative(F_navier_stokes, w, dw)
problem_navier_stokes = NonlinearVariationalProblem(F_navier_stokes, w, bcs, J_navier_stokes)
solver_navier_stokes = NonlinearVariationalSolver(problem_navier_stokes)
solver_navier_stokes.parameters["newton_solver"]["absolute_tolerance"] = 1E-4
solver_navier_stokes.parameters["newton_solver"]["relative_tolerance"] = 1E-3
solver_navier_stokes.parameters["newton_solver"]["maximum_iterations"] = 1000
solver_navier_stokes.parameters["newton_solver"]["linear_solver"] = "mumps"
solver_navier_stokes.solve()
(u_dimensionless, p_dimensionless) = w.split()

dGamma = Measure("ds")(subdomain_data=fluid_boundary_parts)
inlet_area_local = L**2
pressure_drop= (1/inlet_area_local)*p_dimensionless*dGamma(0)
dpressure_drop_dU        = derivative(pressure_drop, p_dimensionless)
dpressure_drop_dU_vector = assemble(dpressure_drop_dU)
dpressure_drop_dU_vector_array = dpressure_drop_dU_vector.get_local()

bcs_new = []
inflow_profile    = Expression(('0','0','0'), u_inlet=u_inlet, degree=2)
bcu_inlet_new         = DirichletBC(W.sub(0), inflow_profile, inlet_front)
bcs_new = [bcu_no_slip_right,  bcu_no_slip_left,\
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet_new,      bcu_back_pressure]

dResidual_dU = derivative(F_navier_stokes, w, dw)
dResidual_dU_matrix = assemble(dResidual_dU)
A, b = assemble_system(dResidual_dU, -dpressure_drop_dU, bcs_new)

The error I get is:

*** Error:   Unable to assemble system.
*** Reason:  expected forms (a, L) to share a FunctionSpace.
*** Where:   This error was encountered inside SystemAssembler.cpp.

Can someone help me fixing this?
Thank you.

Be careful about the use of split(w) and w.split(), see Function.split() vs Split(Function) - #2 by dokken
I would remove the w.split() and simply use p instead of p_dimensionless.
Also derivative should take 3 arguments: a form, a function and a direction, not sure what is happening with only 2 arguments in derivative(pressure_drop, p_dimensionless).
Finally, J_navier_stokes is the same as dResidual_dU, so be consistent in your definitions, avoid redefining everything

2 Likes

Hi @dokken,

can you check on this?

As you havenâ€™t properly formatted your code, I do not have bandwidth to go through it. Please properly format your code.

I would also start by looking at the solution in Paraview. do they look similar? Are boundary conditions properly applied?

Hi Dokken,

Use 3x ` encapsulation ie

```python
```
1 Like

Hey @dokken,

The solution for legacy FEniCS after @bleyerj recommendations worked well and the updated code looks like below:

from fenics import *
import numpy as np
import ufl as ufl
from mpi4py import MPI

mu_fluid = 0.0089
rho_fluid = 1.0
L = 1.0
u_inlet = 0.1
p_outlet = 0.0
Re = rho_fluid*u_inlet*L/mu_fluid
lx, ly, lz = 2.0, 1.0, 1.0
nelx, nely, nelz = 20, 10, 10
mesh = BoxMesh.create(MPI.COMM_WORLD, [Point(0.0, 0.0, 0.0), Point(lx, ly, lz)],
[nelx, nely, nelz],
CellType.Type.hexahedron)

V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
solution_space = MixedElement([V,Q])
W = FunctionSpace(mesh, solution_space)

bcs = []
fluid_boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
fluid_boundary_parts.set_all(3)
no_slip_boundary_right  = CompiledSubDomain("near(x[1],0.0)")
no_slip_boundary_left   = CompiledSubDomain("near(x[1],ly)", ly = ly)
no_slip_boundary_bottom = CompiledSubDomain("near(x[2],0.0)")
no_slip_boundary_top    = CompiledSubDomain("near(x[2],lz)", lz=lz)
inlet_front = CompiledSubDomain("on_boundary && near(x[0],0.0)")
outlet_back = CompiledSubDomain("on_boundary && near(x[0],lx)", lx=lx)
no_slip_boundary_right.mark(fluid_boundary_parts, 2)
no_slip_boundary_left.mark(fluid_boundary_parts, 2)
no_slip_boundary_bottom.mark(fluid_boundary_parts, 2)
no_slip_boundary_top.mark(fluid_boundary_parts, 2)
inlet_front.mark(fluid_boundary_parts, 0)
outlet_back.mark(fluid_boundary_parts, 1)
bcu_no_slip_right  = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_right)
bcu_no_slip_left   = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_left)
bcu_no_slip_bottom = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_bottom)
bcu_no_slip_top    = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_top)
inflow_profile    = Expression(('u_inlet','0','0'), u_inlet=u_inlet, degree=2)
bcu_inlet         = DirichletBC(W.sub(0), inflow_profile, inlet_front)
bcu_back_pressure = DirichletBC(W.sub(1), Constant(0), outlet_back)
bcs = [bcu_no_slip_right,  bcu_no_slip_left,\
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet,          bcu_back_pressure]

w        = Function(W)
dw       = TrialFunction(W)
(u, p)   = split(w)
(du, dp) = TrialFunctions(W)
(v, q)   = TestFunctions(W)
F_navier_stokes = (
- inner(p,div(v))*dx
+ q*div(u)*dx
)
J_navier_stokes = derivative(F_navier_stokes, w, dw)
problem_navier_stokes = NonlinearVariationalProblem(F_navier_stokes, w, bcs, J_navier_stokes)
solver_navier_stokes = NonlinearVariationalSolver(problem_navier_stokes)
solver_navier_stokes.parameters["newton_solver"]["absolute_tolerance"] = 1E-4
solver_navier_stokes.parameters["newton_solver"]["relative_tolerance"] = 1E-3
solver_navier_stokes.parameters["newton_solver"]["maximum_iterations"] = 1000
solver_navier_stokes.parameters["newton_solver"]["linear_solver"] = "mumps"
solver_navier_stokes.solve()

dGamma = Measure("ds")(subdomain_data=fluid_boundary_parts)
inlet_area_local = L**2
pressure_drop = (1/inlet_area_local)*p*dGamma(0)
dpressure_drop_dU        = derivative(pressure_drop, p, dp)
dpressure_drop_dU_vector = assemble(dpressure_drop_dU)

bcs_new = []
inflow_profile    = Expression(('0','0','0'), u_inlet=u_inlet, degree=2)
bcu_inlet_new         = DirichletBC(W.sub(0), inflow_profile, inlet_front)
bcs_new = [bcu_no_slip_right,  bcu_no_slip_left,\
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet_new,      bcu_back_pressure]

dResidual_dU = J_navier_stokes
A, b = assemble_system(dResidual_dU, -dpressure_drop_dU, bcs_new)

solution      = Function(W)
linear_solver = KrylovSolver('default', 'ilu')
linear_solver.parameters["relative_tolerance"] = 1E-5
linear_solver.solve(A, solution.vector(), b)
print("Norm of solution: ", norm(solution.vector(), 'l2'), flush=True)
print("Solution[23]:     ", solution.vector()[23], flush=True)
where the norm of solution and value at element 23 of dolfin mesh are:

Norm of solution:  2.55711615777551e-35
Solution[23]:      -9.427552233559459e-38

In the translated code to dolfinx syntax, the solution of the element of dolfinx mesh located at the same position of element 23 (numbering is different between dolfin and dolfinx) is completely different. The code in dolfinx looks like below. Iâ€™m using dolfinx 0.6.0:

from mpi4py import MPI
from dolfinx.fem import (
Function, FunctionSpace, form, locate_dofs_topological, dirichletbc,
Constant, assemble_scalar
)
from dolfinx.fem.petsc import (
apply_lifting, assemble_matrix, assemble_vector, create_vector,
set_bc, NonlinearProblem
)
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.mesh import (
CellType, create_box, locate_entities_boundary, meshtags
)
from ufl import (
dx, grad, inner, derivative, TrialFunction, TrialFunctions, Measure, dot,
)
from petsc4py.PETSc import ScalarType
from dolfinx.io import XDMFFile
import numpy as np
from dolfinx import log
from petsc4py import PETSc

mu_fluid = 0.0089
rho_fluid = 1.0
L = 1.0
u_inlet = 0.1
p_outlet = 0.0
alpha_solid = 1.0e6
alpha_fluid = 0.0
lx, ly, lz = 2.0, 1.0, 1.0
nelx, nely, nelz = 20, 10, 10
domain = create_box(MPI.COMM_WORLD, [(0., 0., 0.), (lx, ly, lz)], [nelx, nely, nelz], CellType.hexahedron)

V = VectorElement("CG", domain.ufl_cell(), 2)
Q = FiniteElement("CG", domain.ufl_cell(), 1)
solution_space = MixedElement([V,Q])
W = FunctionSpace(domain, solution_space)
D = FunctionSpace(domain, ("DG", 0))
volume_fractions = Function(D)

tolerance = 1e-6
boundary_conditions_navier_stokes = []
def no_slip_boundary_right(x):
return (x[1] < tolerance)
def no_slip_boundary_left(x):
return (x[1] > (ly - tolerance))
def no_slip_boundary_bottom(x):
return (x[2] < tolerance)
def no_slip_boundary_top(x):
return (x[2] > (lz - tolerance))
def inlet_front(x):
return (x[0] < tolerance)
def outlet_back(x):
return (x[0] > (lx - tolerance))
def inlet_velocity(x):
values = np.zeros((3, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = u_inlet
return values
def inlet_velocity_homogeneous(x):
values = np.zeros((3, x.shape[1]), dtype=PETSc.ScalarType)
return values

facets_no_slip_boundary_right  = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_right)
facets_no_slip_boundary_left   = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_left)
facets_no_slip_boundary_bottom = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_bottom)
facets_no_slip_boundary_top    = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_top)
facets_inlet_front             = locate_entities_boundary(domain, domain.topology.dim - 1, inlet_front)
facets_outlet_back             = locate_entities_boundary(domain, domain.topology.dim - 1, outlet_back)
facet_indices = [facets_no_slip_boundary_right,  facets_no_slip_boundary_left,\
facets_no_slip_boundary_bottom, facets_no_slip_boundary_top,\
facets_inlet_front,             facets_outlet_back]
facet_markers = []
facet_markers.append(np.full_like(facet_indices[0], 2))
facet_markers.append(np.full_like(facet_indices[1], 2))
facet_markers.append(np.full_like(facet_indices[2], 2))
facet_markers.append(np.full_like(facet_indices[3], 2))
facet_markers.append(np.full_like(facet_indices[4], 0))
facet_markers.append(np.full_like(facet_indices[5], 1))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tags    = meshtags(domain, domain.topology.dim - 1, facet_indices[sorted_facets], facet_markers[sorted_facets])
B, B_to_W = W.sub(0).collapse()
u_no_slip = Function(B)
u_no_slip.x.set(0)
bcu_no_slip_right  = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_right), W.sub(0))
bcu_no_slip_left   = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_left), W.sub(0))
bcu_no_slip_bottom = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_bottom), W.sub(0))
bcu_no_slip_top    = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_top), W.sub(0))
u_inlet_profile    = Function(B)
u_inlet_profile.interpolate(inlet_velocity)
bcu_inlet          = dirichletbc(u_inlet_profile, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_inlet_front), W.sub(0))
G, G_to_W          = W.sub(1).collapse()
zero_pressure      = Function(G)
zero_pressure.x.set(0)
bcu_back_pressure = dirichletbc(zero_pressure, locate_dofs_topological((W.sub(1), G), domain.topology.dim - 1, facets_outlet_back), W.sub(1))
boundary_conditions_navier_stokes = [bcu_no_slip_right,  bcu_no_slip_left,\
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet,          bcu_back_pressure]
def alpha(volume_fractions):
return alpha_solid + (alpha_fluid - alpha_solid)*volume_fractions
volume_fractions.x.set(1.0)
w        = Function(W)
dw       = TrialFunction(W)
(u, p)   = split(w)
(du, dp) = TrialFunctions(W)
(v, q)   = TestFunctions(W)
F_navier_stokes = (
- inner(p,div(v))*dx
+ alpha(volume_fractions)*dot(v,u)*dx
+ q*div(u)*dx
)
J_navier_stokes = derivative(F_navier_stokes, w, dw)

problem_navier_stokes = NonlinearProblem(F_navier_stokes, w,
boundary_conditions_navier_stokes,
J_navier_stokes)
solver_navier_stokes = NewtonSolver(MPI.COMM_WORLD, problem_navier_stokes)
solver_navier_stokes.convergence_criteria = "incremental"
solver_navier_stokes.report = True
solver_navier_stokes.atol   = 1E-4
solver_navier_stokes.rtol   = 1E-3
solver_navier_stokes.max_it = 1000
krylov_subspace = solver_navier_stokes.krylov_solver
options = PETSc.Options()
option_prefix = krylov_subspace.getOptionsPrefix()
options[f"{option_prefix}kps_type"] = "preonly"
krylov_subspace.setFromOptions()
newton_solver_iterations, converged = solver_navier_stokes.solve(w)
assert (converged)
w.x.scatter_forward()

dGamma = Measure("ds", domain=domain, subdomain_data=facet_tags, metadata=None)
inlet_area             = L**2
pressure_drop          = (1/inlet_area)*p*dGamma(0)
dpressure_drop_dp      = derivative(pressure_drop, p, dp)
dpressure_drop_dp_form = form(dpressure_drop_dp)

u_inlet_profile = Function(B)
u_inlet_profile.interpolate(inlet_velocity_homogeneous)
bcu_inlet_new   = dirichletbc(u_inlet_profile, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_inlet_front), W.sub(0))
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet_new,      bcu_back_pressure]

solution            = Function(W)
dresidual_dU        = J_navier_stokes
dresidual_dU_form   = form(dresidual_dU)
dresidual_dU_matrix.assemble()
dpressure_drop_dp_vector = create_vector(dpressure_drop_dp_form)
assemble_vector(dpressure_drop_dp_vector, dpressure_drop_dp_form)

options = PETSc.Options()
linear_solver = PETSc.KSP().create(MPI.COMM_WORLD)
linear_solver.setType("preonly")
linear_solver.setTolerances(rtol=1e-8)
linear_solver.getPC().setType("ilu")
options = PETSc.Options()
linear_solver.setFromOptions()
linear_solver.setInitialGuessNonzero(False)
linear_solver.setOperators(dresidual_dU_matrix)
linear_solver.solve(dpressure_drop_dp_vector, solution.vector)
solution.x.scatter_forward()
print("Norm of solution:   ", np.linalg.norm(solution.x.array), flush=True)
print("solution[21]:       ", solution.x.array[21])
where the norm of solution and value obtained at element 21 of dolfinx mesh are:

Norm of solution:    237.3371158413562
solution[21]:        -9.345168135762837

Iâ€™ve already checked the solution in paraview, they are the same. I also checked the boundary conditions application and couldnâ€™t find any difference. Actually, I printed the locations where I apply the boundary conditions to confirm they were being properly selected. For me it seems a matter of syntax, i.e., on how to work with mixed elements in dolfinx. What do you think? Thank you.

Please note that the solution w seems similar for both problems (comparing the norm).
This means that you are likely not solving the same problem in the adjoint step.

If you change the solver for the adjoint solution to â€śluâ€ť, instead of â€śiluâ€ť you will see a dramatic change in the solution (and you will get consistent norms for the solution.vector().

Next you will note that the local numbering [21] and [23] is not consistent, i.e. you are not evaluating at the same point in both problems.

Opening either of the original ILU solutions in paraview seems to yield unphysical solutions.

Hi @dokken,

Iâ€™m finally back at this issue with an update after your comments. Previously I was checking the norm, now Iâ€™m checking total derivative values at certain elements. The code I have written after your comments is:

from mpi4py import MPI
from dolfinx.fem import (
Function, FunctionSpace, form, locate_dofs_topological, dirichletbc,
assemble_scalar
)
from dolfinx.fem.petsc import (
apply_lifting, assemble_matrix, assemble_vector, create_vector,
set_bc, NonlinearProblem
)
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.mesh import (
CellType, create_box, locate_entities_boundary, meshtags
)
from ufl import (
dx, grad, inner, derivative, TrialFunction, TrialFunctions, Measure, dot,
)
from petsc4py.PETSc import ScalarType
import numpy as np
from petsc4py import PETSc

comm = MPI.COMM_WORLD

mu_fluid = 0.01
rho_fluid = 1110
L = 1.0
u_inlet = 0.000022523
p_outlet = 0.0
alpha_solid = 1.0e6
alpha_fluid = 0.0
print("Reynolds: ", rho_fluid*u_inlet*L/mu_fluid, flush=True)
lx, ly, lz = 2.0, 1.0, 1.0
nelx, nely, nelz = 40, 20, 20
domain = create_box(MPI.COMM_WORLD, [(0., 0., 0.), (lx, ly, lz)], [nelx, nely, nelz], CellType.hexahedron)

V = VectorElement("CG", domain.ufl_cell(), 2)
Q = FiniteElement("CG", domain.ufl_cell(), 1)
solution_space = MixedElement([V,Q])
W = FunctionSpace(domain, solution_space)
D = FunctionSpace(domain, ("DG", 0))
volume_fractions = Function(D)

tolerance = 1e-6
boundary_conditions_navier_stokes = []
def no_slip_boundary_right(x):
return (x[1] < tolerance)
def no_slip_boundary_left(x):
return (x[1] > (ly - tolerance))
def no_slip_boundary_bottom(x):
return (x[2] < tolerance)
def no_slip_boundary_top(x):
return (x[2] > (lz - tolerance))
def inlet_front(x):
return (x[0] < tolerance)
def outlet_back(x):
return (x[0] > (lx - tolerance))
def inlet_velocity(x):
values = np.zeros((3, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = u_inlet
return values
def inlet_velocity_homogeneous(x):
values = np.zeros((3, x.shape[1]), dtype=PETSc.ScalarType)
return values

facets_no_slip_boundary_right  = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_right)
facets_no_slip_boundary_left   = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_left)
facets_no_slip_boundary_bottom = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_bottom)
facets_no_slip_boundary_top    = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_top)
facets_inlet_front             = locate_entities_boundary(domain, domain.topology.dim - 1, inlet_front)
facets_outlet_back             = locate_entities_boundary(domain, domain.topology.dim - 1, outlet_back)
facet_indices = [facets_no_slip_boundary_right,  facets_no_slip_boundary_left,\
facets_no_slip_boundary_bottom, facets_no_slip_boundary_top,\
facets_inlet_front,             facets_outlet_back]
facet_markers = []
facet_markers.append(np.full_like(facet_indices[0], 2))
facet_markers.append(np.full_like(facet_indices[1], 2))
facet_markers.append(np.full_like(facet_indices[2], 2))
facet_markers.append(np.full_like(facet_indices[3], 2))
facet_markers.append(np.full_like(facet_indices[4], 0))
facet_markers.append(np.full_like(facet_indices[5], 1))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tags    = meshtags(domain, domain.topology.dim - 1, facet_indices[sorted_facets], facet_markers[sorted_facets])
B, B_to_W = W.sub(0).collapse()
u_no_slip = Function(B)
u_no_slip.x.set(0)
bcu_no_slip_right  = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_right), W.sub(0))
bcu_no_slip_left   = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_left), W.sub(0))
bcu_no_slip_bottom = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_bottom), W.sub(0))
bcu_no_slip_top    = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_top), W.sub(0))
u_inlet_profile    = Function(B)
u_inlet_profile.interpolate(inlet_velocity)
bcu_inlet          = dirichletbc(u_inlet_profile, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_inlet_front), W.sub(0))
G, G_to_W          = W.sub(1).collapse()
zero_pressure      = Function(G)
zero_pressure.x.set(0)
bcu_back_pressure = dirichletbc(zero_pressure, locate_dofs_topological((W.sub(1), G), domain.topology.dim - 1, facets_outlet_back), W.sub(1))
boundary_conditions_navier_stokes = [bcu_no_slip_right,  bcu_no_slip_left,\
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet,          bcu_back_pressure]
def alpha(volume_fractions):
return alpha_solid + (alpha_fluid - alpha_solid)*volume_fractions
volume_fractions.x.set(1.0)
w        = Function(W)
dw       = TrialFunction(W)
(u, p)   = split(w)
(du, dp) = TrialFunctions(W)
(v, q)   = TestFunctions(W)
F_navier_stokes = (
- inner(p,div(v))*dx
+ alpha(volume_fractions)*dot(v,u)*dx
+ q*div(u)*dx
)
J_navier_stokes = derivative(F_navier_stokes, w, dw)

problem_navier_stokes = NonlinearProblem(F_navier_stokes, w,
boundary_conditions_navier_stokes,
J_navier_stokes)
solver_navier_stokes = NewtonSolver(MPI.COMM_WORLD, problem_navier_stokes)
solver_navier_stokes.convergence_criteria = "incremental"
solver_navier_stokes.report = True
solver_navier_stokes.atol   = 1E-50
solver_navier_stokes.rtol   = 1E-7
solver_navier_stokes.max_it = 1000
krylov_subspace = solver_navier_stokes.krylov_solver
options = PETSc.Options()
option_prefix = krylov_subspace.getOptionsPrefix()
options[f"{option_prefix}kps_type"] = "preonly"
options[f"{option_prefix}pc_type"]  = "lu"
krylov_subspace.setFromOptions()
newton_solver_iterations, converged = solver_navier_stokes.solve(w)
assert (converged)
w.x.scatter_forward()

dGamma = Measure("ds", domain=domain, subdomain_data=facet_tags, metadata=None)
inlet_area_local       = assemble_scalar(form(1*dGamma(0)))
inlet_area_global      = comm.allreduce(inlet_area_local, op=MPI.SUM)
pressure_drop          = (1/inlet_area_global)*p*dGamma(0)
dpressure_drop_dp      = derivative(pressure_drop, p, dp)
dpressure_drop_dp_form = form(dpressure_drop_dp)

u_inlet_profile = Function(B)
u_inlet_profile.interpolate(inlet_velocity_homogeneous)
bcu_inlet_new   = dirichletbc(u_inlet_profile, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_inlet_front), W.sub(0))
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet_new,      bcu_back_pressure]
solution            = Function(W)
dresidual_dU        = J_navier_stokes
dresidual_dU_form   = form(dresidual_dU)
dresidual_dU_matrix.assemble()
dpressure_drop_dp_vector = create_vector(dpressure_drop_dp_form)
assemble_vector(dpressure_drop_dp_vector, dpressure_drop_dp_form)
options = PETSc.Options()
linear_solver = PETSc.KSP().create(MPI.COMM_WORLD)
linear_solver.setType("preonly")
linear_solver.getPC().setType("ilu")
linear_solver.setTolerances(rtol=1e-10)
options = PETSc.Options()
linear_solver.setFromOptions()
linear_solver.setInitialGuessNonzero(False)
linear_solver.setOperators(dresidual_dU_matrix)
linear_solver.solve(dpressure_drop_dp_vector, solution.vector)
solution.x.scatter_forward()

dresidual_dvol_matrix.assemble()
DJ_Dvol = Function(D)
DJ_Dvol.vector.array = dresidual_dvol_matrix*solution.vector

print("Sensitivity of element 3189:  ", DJ_Dvol.x.array[3189], flush=True)
print("Sensitivity of element 2272:  ", DJ_Dvol.x.array[2272], flush=True)
print("Sensitivity of element 9380:  ", DJ_Dvol.x.array[9380], flush=True)

This code outputs the following results:
Sensitivity of element 3189: 7.148662250159949e-05
Sensitivity of element 2272: -0.0028377509913678333
Sensitivity of element 9380: -4.144243370744519e-06

A comparison with finite difference is established with finite difference at the same element IDs and the results seem off to me:
Finite difference of pressure drop function in element 3189: -0.0016184140569784253
Finite difference of pressure drop function in element 2272: -0.0029172783317461728
Finite difference of pressure drop function in element 9380: -6.039570156924495e-06
Do you have any idea what is going? I mean, I expected a much closer result, and element 3189 seems completely off.

Thank you.

You are still using ilu for solving the adjoint

and havenâ€™t commented on my feedback from last time around

Without any discussion as to why you have disregarded this information, I do not have the bandwidth to debug your rather lengthy code.