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

Legacy DOLFIN
1

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 *
from fenics_adjoint 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)
parameters["form_compiler"]["quadrature_degree"] = 2

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 = (
    rho_fluid*dot(dot(u,nabla_grad(u)),v)*dx
    + mu_fluid*inner(grad(u),grad(v))*dx
    - 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.

2

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
3
4

Hi @dokken,

can you check on this?

5

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?

6

Hi Dokken,

thank you for your reply. Just so I can properly try to address your recommendations, could you please elaborate a little bit more on what you call “properly format your code”?

7

Use 3x ` encapsulation ie

```python
#add code here
```
1 Like
8

Hey @dokken,

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

from fenics import *
from fenics_adjoint 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)
parameters["form_compiler"]["quadrature_degree"] = 2

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 = (
    rho_fluid*dot(dot(u,nabla_grad(u)),v)*dx
    + mu_fluid*inner(grad(u),grad(v))*dx
    - 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,
    adjoint, VectorElement, FiniteElement, MixedElement, TestFunctions,
    CellDiameter, nabla_grad, div, split
)
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
radius = L
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)
dx = Measure("dx", domain=domain, metadata={"quadrature_degree": 2})
F_navier_stokes = (
    rho_fluid*dot(dot(u,nabla_grad(u)),v)*dx
    + mu_fluid*inner(grad(u),grad(v))*dx
    - 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)

bcs_adjoint = []
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))
bcs_adjoint     = [bcu_no_slip_right,  bcu_no_slip_left,\
                   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_matrix(dresidual_dU_form, bcs_adjoint)
dresidual_dU_matrix.assemble()
dpressure_drop_dp_vector = create_vector(dpressure_drop_dp_form)
assemble_vector(dpressure_drop_dp_vector, dpressure_drop_dp_form)
apply_lifting(dpressure_drop_dp_vector, [dresidual_dU_form], [bcs_adjoint])
dpressure_drop_dp_vector.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(dpressure_drop_dp_vector, bcs_adjoint)

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.

9

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.