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)