Adjoint method using nonlinear residual and mixed elements

Hi,

I’m trying to implement the adjoint method using a nonlinear residual and mixed elements. I based my code in the following post: Compute the sensitivity in parallel with DOLFINx
However, when I compare my implementation with finite difference I get almost an order of magnitude in relative error. Can someone help me figuring out why such an unexpected difference? Below follow my code. I’m currently running dolfinx 0.6.0.

Thank you.

from mpi4py import MPI
from dolfinx.fem import (
    Function, FunctionSpace, form, locate_dofs_topological, dirichletbc,
    assemble_scalar
)
from dolfinx.fem.petsc import (
    assemble_matrix, assemble_vector, set_bc, LinearProblem, 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,
    nabla_grad, div, split
)
from petsc4py.PETSc import ScalarType
import numpy as np
from petsc4py import PETSc
import warnings

comm = MPI.COMM_WORLD
number_processors = comm.Get_size()
warnings.filterwarnings("ignore")

mu_fluid = 0.0089
rho_fluid = 1.0           
L = 1.0             
u_inlet = 0.1             
p_outlet = 0.0

(lx, ly, lz) = (2.0, 1.0, 1.0)
(nelx, nely, nelz) = (30, 15, 15)
domain = create_box(comm=MPI.COMM_WORLD, points=[(0., 0., 0.), (lx, ly, lz)], n=[nelx, nely, nelz], cell_type=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))
design_variable = 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

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, _      = 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, _              = 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]

design_variable.vector.array = np.ones(int(nelx*nely*nelz))
design_variable.x.scatter_forward()
w      = Function(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)
dw     = TrialFunction(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
    + design_variable*dot(v,u)*dx
    + q*div(u)*dx
)
jacobian_navier_stokes = derivative(F_navier_stokes, w, dw)
# jacobian_navier_stokes = derivative(F_navier_stokes, w)
problem_navier_stokes = NonlinearProblem(F_navier_stokes, w, boundary_conditions_navier_stokes, jacobian_navier_stokes)
solver_navier_stokes  = NewtonSolver(domain.comm, problem_navier_stokes)
solver_navier_stokes.convergence_criteria = "incremental"
solver_navier_stokes.report = True
solver_navier_stokes.rtol = 1E-5
solver_navier_stokes.max_it = 400
krylov_subspace = solver_navier_stokes.krylov_solver
options = PETSc.Options()
option_prefix = krylov_subspace.getOptionsPrefix()
options[f"{option_prefix}ksp_type"] = "preonly"
krylov_subspace.setFromOptions()
newton_solver_iterations, converged = solver_navier_stokes.solve(w)
assert (converged)
(u_solution, p_solution) = w.split()

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)
# Calculate pressure drop
pressure_drop = (1/inlet_area_global)*w.sub(1)*dGamma(0)
J = pressure_drop

adjoint_lambdas   = Function(W)
dpressure_drop_dU = -derivative(pressure_drop, w)
dResidual_dU      = adjoint(derivative(F_navier_stokes, w)) 
problem = LinearProblem(dResidual_dU, dpressure_drop_dU, boundary_conditions_navier_stokes)
adjoint_lambdas = problem.solve()

dJ_dvol        = assemble_vector(form(derivative(J, design_variable)))
dJ_dvol.assemble()
dResidual_dvol = assemble_matrix(form(adjoint(derivative(F_navier_stokes, design_variable))))
dResidual_dvol.assemble()

dJ_dvol = dJ_dvol + dResidual_dvol*adjoint_lambdas.vector
dJ_dvol.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
print("dolfinx sensitivity element 41: ", dJ_dvol.array[41], flush=True)

# Finite difference method starts here
J_local   = assemble_scalar(form(J))
J_global  = comm.allreduce(J_local, op=MPI.SUM)
objective = [J_global]
print('Pressure drop inlet     = {:.16e}'.format(objective[0]))
if number_processors > 1:
    raise RuntimeError("Sensitivity check must be performed on a single process.")

artificial_perturbation = [1e-4]
element_id = [41]
initial_design_variable = design_variable.vector.array.copy()
DJ_dolfinx           = [dJ_dvol.array[element_id]]
DJ_finite_difference = [np.zeros((len(element_id), len(artificial_perturbation))) for _ in range(len(DJ_dolfinx))]
error_absolute       = [np.zeros((len(element_id), len(artificial_perturbation))) for _ in range(len(DJ_dolfinx))]
error_relative       = [np.zeros((len(element_id), len(artificial_perturbation))) for _ in range(len(DJ_dolfinx))]
for i in range(len(element_id)):
    for j in range(len(artificial_perturbation)):
        design_variable.vector.array[:] = initial_design_variable.copy()
        design_variable.vector.array[element_id[i]] -= artificial_perturbation[j]
        design_variable.x.scatter_forward()

        w      = Function(W)
        (u, p) = split(w)
        (v, q) = TestFunctions(W)
        dw     = TrialFunction(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
            + design_variable*dot(v,u)*dx
            + q*div(u)*dx
        )
        jacobian_navier_stokes = derivative(F_navier_stokes, w, dw)
        # jacobian_navier_stokes = derivative(F_navier_stokes, w)
        problem_navier_stokes = NonlinearProblem(F_navier_stokes, w, boundary_conditions_navier_stokes, jacobian_navier_stokes)
        solver_navier_stokes  = NewtonSolver(domain.comm, problem_navier_stokes)
        solver_navier_stokes.convergence_criteria = "incremental"
        solver_navier_stokes.report = True
        solver_navier_stokes.rtol = 1E-5
        solver_navier_stokes.max_it = 400
        krylov_subspace = solver_navier_stokes.krylov_solver
        options = PETSc.Options()
        option_prefix = krylov_subspace.getOptionsPrefix()
        options[f"{option_prefix}ksp_type"] = "preonly"
        krylov_subspace.setFromOptions()
        newton_solver_iterations, converged = solver_navier_stokes.solve(w)
        assert (converged)
        (u_solution, p_solution) = w.split()

        pressure_drop_new = (1/inlet_area_global)*w.sub(1)*dGamma(0)
        J_new = pressure_drop_new
        J_local_new = assemble_scalar(form(J_new))
        J_global_new = comm.allreduce(J_local_new, op=MPI.SUM)
        objective_perturbed = [J_global_new]
        print('Pressure drop inlet new = {:.16e}'.format(objective_perturbed[0]))

        print("Artificial Perturbation - Element id - Dolfinx - Finite Difference - Absolute Difference - Relative Difference", flush=True)
        for k in range(len(objective)):
            DJ_finite_difference[k][i, j] = (objective[k] - objective_perturbed[k])/artificial_perturbation[j]
            error_absolute[k][i, j]       = np.abs(DJ_dolfinx[k][i] - DJ_finite_difference[k][i, j])
            error_relative[k][i, j]       = error_absolute[k][i, j]/np.abs(DJ_dolfinx[k][i])
            print("%10.2e - %d - %10.16e - %10.16e - %10.8e - %10.8e" % (artificial_perturbation[j],\
                      element_id[k], DJ_dolfinx[k][i], DJ_finite_difference[k][i,j],\
                      error_absolute[k][i,j], error_relative[k][i,j]), flush=True)

Hi there,

I have an update regarding this question. Using mixed elements, I create my cost function like below:

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)*w.sub(1)*dGamma(0)

but when I take the derivative of that cost function with respect to the state variables, like below:

dpressure_drop_dU = -derivative(pressure_drop, w)

I get a vector of zero values:

Vector of derivatives:  [0. 0. 0. ... 0. 0. 0.]
Norm of vector of derivatives 0.0

Can someone help me figuring out why using this mixed formulation my cost function seems to have no relationship with respect to my state variables? My full code is:

from mpi4py import MPI
from dolfinx.fem import (
    Function, FunctionSpace, form, locate_dofs_topological, dirichletbc,
    assemble_scalar
)
from dolfinx.fem.petsc import (
    assemble_vector, 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,
    nabla_grad, div, split
)
from petsc4py.PETSc import ScalarType
import numpy as np
from petsc4py import PETSc
import warnings

comm = MPI.COMM_WORLD
warnings.filterwarnings("ignore")

mu_fluid = 0.0089
rho_fluid = 1.0           
L = 1.0             
u_inlet = 0.1             
p_outlet = 0.0

(lx, ly, lz) = (2.0, 1.0, 1.0)
(nelx, nely, nelz) = (30, 15, 15)
domain = create_box(comm, points=[(0., 0., 0.), (lx, ly, lz)], n=[nelx, nely, nelz], cell_type=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))
design_variable = 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

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, _      = 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, _              = 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]

design_variable.vector.array = np.ones(int(nelx*nely*nelz))
design_variable.x.scatter_forward()
w      = Function(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)
dw     = TrialFunction(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
    + design_variable*dot(v,u)*dx
    + q*div(u)*dx
)
jacobian_navier_stokes = derivative(F_navier_stokes, w, dw)
problem_navier_stokes = NonlinearProblem(F_navier_stokes, w, boundary_conditions_navier_stokes, jacobian_navier_stokes)
solver_navier_stokes  = NewtonSolver(domain.comm, problem_navier_stokes)
solver_navier_stokes.convergence_criteria = "incremental"
solver_navier_stokes.report = True
solver_navier_stokes.rtol = 1E-5
solver_navier_stokes.max_it = 400
krylov_subspace = solver_navier_stokes.krylov_solver
options = PETSc.Options()
option_prefix = krylov_subspace.getOptionsPrefix()
options[f"{option_prefix}ksp_type"] = "preonly"
krylov_subspace.setFromOptions()
newton_solver_iterations, converged = solver_navier_stokes.solve(w)
assert (converged)
(u_solution, p_solution) = w.split()

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)*w.sub(1)*dGamma(0)
dpressure_drop_dU = -derivative(pressure_drop, w)
dpressure_drop_dU_vector = assemble_vector(form(dpressure_drop_dU))
dpressure_drop_dU_vector.assemble()
print("Vector of derivatives: ", dpressure_drop_dU_vector.array, flush=True)
print("Norm of vector of derivatives", np.linalg.norm(dpressure_drop_dU_vector.array), flush=True)