Mistake in calculating sensitivity of boundary condition

General
1

I am trying to compute sensitivity of boundary condition using the adjoint method but getting the procedure wrong somewhere which I am unable to guess. I will appreciate if someone could point out my mistake.

For the MWE, I solve Laplace equation on a domain with two electrodes (holes) where the potential is applied. Due to geometrical symmetry and equal but opposite potential on the two, I expect the calculated sensitivity to also be equal but opposite in sign. But the result from the adjoint method suggests dissimilar magnitude. Finite difference calculation seems to pass the sanity check and suggests a number that is order of magnitude larger than the one suggested by adjoint method.

Thanks in advance for the help! Below is the MWE evaluated on dolfinx 0.9

import gmsh, dolfinx, ufl, mpi4py, petsc4py
from dolfinx.io import gmshio
import dolfinx.fem.petsc, dolfinx.nls.petsc
import matplotlib.pyplot as plt
import scifem
import numpy as np

gmsh.finalize()
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gdim = 2
occ = gmsh.model.occ

c1 = occ.addDisk(-0.2, 0, 0, 0.1, 0.1)
c2 = occ.addDisk(0.2, 0, 0, 0.1, 0.1)
fov = occ.addDisk(0, 0, 0, 1, 1)
cut = occ.cut([(gdim, fov)], [(gdim, c1), (gdim, c2)])
occ.synchronize()

# add_sub physical tags
all_doms = occ.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = occ.getEntities(gdim-1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(2)
gmsh.model.mesh.refine()

mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, 0, gdim)

Omega_V = dolfinx.fem.functionspace(mesh, ("CG", 1))
Omega_V0 = scifem.create_real_functionspace(mesh) # for electrodes

V_sol = dolfinx.fem.Function(Omega_V)
V = ufl.TrialFunction(Omega_V)
Vt = ufl.TestFunction(Omega_V)
n = ufl.FacetNormal(mesh)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft)
dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)

R = ufl.inner(ufl.grad(V_sol), ufl.grad(Vt)) * dx
R += ufl.inner(dolfinx.fem.Constant(mesh, 0.0), Vt) * dx

V1 = dolfinx.fem.Function(Omega_V0) # electrode 1 (physical tag 1)
V2 = dolfinx.fem.Function(Omega_V0) # electrode 2 (physical tag 2)
V1.x.array[:] = 1.0
V2.x.array[:] = -1.0

alpha = 10
hE = ufl.Circumradius(mesh)
R += - ufl.inner(ufl.dot(n, ufl.grad(V_sol)), Vt) * ds((1,2)) - ufl.inner(V_sol, ufl.dot(n, ufl.grad(Vt))) * ds((1,2)) + alpha/hE*ufl.inner(V_sol, Vt) * ds((1,2))
R += ufl.inner(n, ufl.grad(Vt)) * V1 * ds(1) - alpha/hE*ufl.inner(V1, Vt) * ds(1)
R += ufl.inner(n, ufl.grad(Vt)) * V2 * ds(2) - alpha/hE*ufl.inner(V2, Vt) * ds(2)

# forward problem
dbc = []
problem = dolfinx.fem.petsc.NonlinearProblem(R, V_sol, dbc)
solver = dolfinx.nls.petsc.NewtonSolver(mpi4py.MPI.COMM_WORLD, problem)
solver.solve(V_sol)

# calculate energy / objective
J = 1/2*ufl.dot(ufl.grad(V_sol), ufl.grad(V_sol)) * dx
J0 = dolfinx.fem.assemble_scalar(dolfinx.fem.form(J))
print("Energy: {}".format(J0))

# solve the adjoint problem
lmbda = dolfinx.fem.Function(Omega_V)
dRdu = ufl.derivative(R, V_sol, V)
dRdu_adj = ufl.adjoint(dRdu)
dJdu = ufl.derivative(J, V_sol, Vt)

V_dbc = dolfinx.fem.Function(Omega_V)
ext_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
ext_dofs = dolfinx.fem.locate_dofs_topological(Omega_V, gdim-1, ext_facets)
dbc_adj = [dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh, 0.0), ext_dofs, Omega_V)]

adj_problem = dolfinx.fem.petsc.LinearProblem(ufl.replace(dRdu_adj, {V_sol: Vt}), -dJdu, u=lmbda, bcs=dbc_adj)
adj_problem.solve()

dJdV1 = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(ufl.derivative(J, V1)))
dJdV1.assemble()

dJdV2 = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(ufl.derivative(J, V2)))
dJdV2.assemble()

dRdV1_form = dolfinx.fem.form(ufl.adjoint(ufl.derivative(R, V1)))
dRdV1 = dolfinx.fem.petsc.assemble_matrix(dRdV1_form)  # partial derivative
dRdV1.assemble()

dRdV2_form = dolfinx.fem.form(ufl.adjoint(ufl.derivative(R, V2)))
dRdV2 = dolfinx.fem.petsc.assemble_matrix(dRdV2_form)  # partial derivative
dRdV2.assemble()

dJdV1 += dRdV1*lmbda.x.petsc_vec
dJdV1.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
print("sensitivity of V1 from adjoint method:", dJdV1.array)

dJdV2 += dRdV2*lmbda.x.petsc_vec
dJdV2.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
print("sensitivity of V2 from adjoint method:", dJdV2.array)

# cross-validation by forward difference
eps = 1e-6
V1.x.array[:] = 1+eps
V2.x.array[:] = -1
solver.solve(V_sol)

J1 = dolfinx.fem.assemble_scalar(dolfinx.fem.form(J))
print("Energy after perturbing V1: {}".format(J1))
print("Sensitivity of V1 from finite difference: {}".format((J1-J0)/eps))

V1.x.array[:] = 1
V2.x.array[:] = -1+eps
solver.solve(V_sol)

J2 = dolfinx.fem.assemble_scalar(dolfinx.fem.form(J))
print("Energy after perturbing V2: {}".format(J2))
print("Sensitivity of V2 from finite difference: {}".format((J2-J0)/eps))
2

First of all, I would suggest that you use a direct solver for your Newton-method (You have not tuned the solver at all.
I would suggest the same for the adjoint problem.

That is the first thing to check (and also add error if not converged to either solver).

3

Thanks a lot for your feedback! I could not quickly figure how to provide petsc options to Nonlinear/NewtonSolver and just changed the formulation to linear. Even specifying direct solver does not change the output. Below is the code:

import gmsh, dolfinx, ufl, mpi4py, petsc4py
from dolfinx.io import gmshio
import dolfinx.fem.petsc, dolfinx.nls.petsc
import matplotlib.pyplot as plt
import scifem
import numpy as np

gmsh.finalize()
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gdim = 2
occ = gmsh.model.occ

c1 = occ.addDisk(-0.2, 0, 0, 0.1, 0.1)
c2 = occ.addDisk(0.2, 0, 0, 0.1, 0.1)
fov = occ.addDisk(0, 0, 0, 1, 1)
cut = occ.cut([(gdim, fov)], [(gdim, c1), (gdim, c2)])
occ.synchronize()

# add_sub physical tags
all_doms = occ.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = occ.getEntities(gdim-1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(2)
gmsh.model.mesh.refine()

mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, 0, gdim)

Omega_V = dolfinx.fem.functionspace(mesh, ("CG", 1))
Omega_V0 = scifem.create_real_functionspace(mesh) # for electrodes

V_sol = dolfinx.fem.Function(Omega_V)
V = ufl.TrialFunction(Omega_V)
Vt = ufl.TestFunction(Omega_V)
n = ufl.FacetNormal(mesh)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft)
dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)

R = ufl.inner(ufl.grad(V), ufl.grad(Vt)) * dx
R += ufl.inner(dolfinx.fem.Constant(mesh, 0.0), Vt) * dx

V1 = dolfinx.fem.Function(Omega_V0) # electrode 1 (physical tag 1)
V2 = dolfinx.fem.Function(Omega_V0) # electrode 2 (physical tag 2)
V1.x.array[:] = 1.0
V2.x.array[:] = -1.0

alpha = 10
hE = ufl.Circumradius(mesh)
R += - ufl.inner(ufl.dot(n, ufl.grad(V)), Vt) * ds((1,2)) - ufl.inner(V, ufl.dot(n, ufl.grad(Vt))) * ds((1,2)) + alpha/hE*ufl.inner(V, Vt) * ds((1,2))
R += ufl.inner(n, ufl.grad(Vt)) * V1 * ds(1) - alpha/hE*ufl.inner(V1, Vt) * ds(1)
R += ufl.inner(n, ufl.grad(Vt)) * V2 * ds(2) - alpha/hE*ufl.inner(V2, Vt) * ds(2)

# forward problem
dbc = []
petsc_options = {
    "-ksp_type": "preonly",
    "-pc_type": "lu",
    "-pc_factor_mat_solver_type": "mumps",
}

fwd_problem = dolfinx.fem.petsc.LinearProblem(ufl.lhs(R), ufl.rhs(R), bcs=dbc, u=V_sol, petsc_options=petsc_options)
fwd_problem.solve()

# calculate energy / objective
J = 1/2*ufl.dot(ufl.grad(V_sol), ufl.grad(V_sol)) * dx
J0 = dolfinx.fem.assemble_scalar(dolfinx.fem.form(J))
print("Energy: {}".format(J0))

# solve the adjoint problem
lmbda = dolfinx.fem.Function(Omega_V)
dRdu = ufl.derivative(ufl.action(R, V_sol), V_sol, V)
dRdu_adj = ufl.adjoint(dRdu)
dJdu = ufl.derivative(J, V_sol, Vt)

V_dbc = dolfinx.fem.Function(Omega_V)
ext_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
ext_dofs = dolfinx.fem.locate_dofs_topological(Omega_V, gdim-1, ext_facets)
dbc_adj = [dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh, 0.0), ext_dofs, Omega_V)]

adj_problem = dolfinx.fem.petsc.LinearProblem(ufl.replace(dRdu_adj, {V_sol: Vt}), -dJdu, u=lmbda, bcs=dbc_adj, petsc_options=petsc_options)
adj_problem.solve()

dJdV1 = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(ufl.derivative(J, V1)))
dJdV1.assemble()

dJdV2 = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(ufl.derivative(J, V2)))
dJdV2.assemble()

dRdV1_form = dolfinx.fem.form(ufl.adjoint(ufl.derivative(R, V1)))
dRdV1 = dolfinx.fem.petsc.assemble_matrix(dRdV1_form)  # partial derivative
dRdV1.assemble()

dRdV2_form = dolfinx.fem.form(ufl.adjoint(ufl.derivative(R, V2)))
dRdV2 = dolfinx.fem.petsc.assemble_matrix(dRdV2_form)  # partial derivative
dRdV2.assemble()

dJdV1 += dRdV1*lmbda.x.petsc_vec
dJdV1.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
print("sensitivity of V1 from adjoint method:", dJdV1.array)

dJdV2 += dRdV2*lmbda.x.petsc_vec
dJdV2.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
print("sensitivity of V2 from adjoint method:", dJdV2.array)

# cross-validation by finite difference
eps = 1e-6
V1.x.array[:] = 1+eps
V2.x.array[:] = -1
fwd_problem.solve()

J1 = dolfinx.fem.assemble_scalar(dolfinx.fem.form(J))
print("Energy after perturbing V1: {}".format(J1))
print("Sensitivity of V1 from finite difference: {}".format((J1-J0)/eps))

V1.x.array[:] = 1
V2.x.array[:] = -1+eps
fwd_problem.solve()

J2 = dolfinx.fem.assemble_scalar(dolfinx.fem.form(J))
print("Energy after perturbing V2: {}".format(J2))
print("Sensitivity of V2 from finite difference: {}".format((J2-J0)/eps))
4

While experimenting, I found that I get correct results if I specify homogeneous boundary condition for the adjoint problem only where the forward problem does not have dirichlet boundaries. If the experts don’t disagree, this problem is solved. Below is the code that does it

import gmsh, dolfinx, ufl, mpi4py, petsc4py
from dolfinx.io import gmshio
import dolfinx.fem.petsc, dolfinx.nls.petsc
import matplotlib.pyplot as plt
import scifem
import numpy as np

gmsh.finalize()
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gdim = 2
occ = gmsh.model.occ

c1 = occ.addDisk(-0.2, 0, 0, 0.1, 0.1)
c2 = occ.addDisk(0.2, 0, 0, 0.1, 0.1)
fov = occ.addDisk(0, 0, 0, 1, 1)
cut = occ.cut([(gdim, fov)], [(gdim, c1), (gdim, c2)])
occ.synchronize()

# add_sub physical tags
all_doms = occ.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = occ.getEntities(gdim-1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(2)
gmsh.model.mesh.refine()

mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, 0, gdim)

Omega_V = dolfinx.fem.functionspace(mesh, ("CG", 1))
Omega_V0 = scifem.create_real_functionspace(mesh) # for electrodes

V_sol = dolfinx.fem.Function(Omega_V)
V = ufl.TrialFunction(Omega_V)
Vt = ufl.TestFunction(Omega_V)
n = ufl.FacetNormal(mesh)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft)
dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)

R = ufl.inner(ufl.grad(V), ufl.grad(Vt)) * dx
R += ufl.inner(dolfinx.fem.Constant(mesh, 0.0), Vt) * dx

V1 = dolfinx.fem.Function(Omega_V0) # electrode 1 (physical tag 1)
V2 = dolfinx.fem.Function(Omega_V0) # electrode 2 (physical tag 2)
V1.x.array[:] = 1.0
V2.x.array[:] = -1.0

alpha = 10
hE = ufl.Circumradius(mesh)
R += - ufl.inner(ufl.dot(n, ufl.grad(V)), Vt) * ds((1,2)) - ufl.inner(V, ufl.dot(n, ufl.grad(Vt))) * ds((1,2)) + alpha/hE*ufl.inner(V, Vt) * ds((1,2))
R += ufl.inner(n, ufl.grad(Vt)) * V1 * ds(1) - alpha/hE*ufl.inner(V1, Vt) * ds(1)
R += ufl.inner(n, ufl.grad(Vt)) * V2 * ds(2) - alpha/hE*ufl.inner(V2, Vt) * ds(2)

# forward problem
dbc = []
petsc_options = {
    "-ksp_type": "preonly",
    "-pc_type": "lu",
    "-pc_factor_mat_solver_type": "mumps",
}

fwd_problem = dolfinx.fem.petsc.LinearProblem(ufl.lhs(R), ufl.rhs(R), bcs=dbc, u=V_sol, petsc_options=petsc_options)
fwd_problem.solve()

# calculate energy / objective
J = 1/2*ufl.dot(ufl.grad(V_sol), ufl.grad(V_sol)) * dx
J0 = dolfinx.fem.assemble_scalar(dolfinx.fem.form(J))
print("Energy: {}".format(J0))

# solve the adjoint problem
lmbda = dolfinx.fem.Function(Omega_V)
dRdu = ufl.derivative(ufl.action(R, V_sol), V_sol, V)
dRdu_adj = ufl.adjoint(dRdu)
dJdu = ufl.derivative(J, V_sol, Vt)

V_dbc = dolfinx.fem.Function(Omega_V)
ext_facets = np.hstack([ft.find(j) for j in (3,)]) # apply dirichlet boundary condition only to Neumann boundaries in the forward problem
ext_dofs = dolfinx.fem.locate_dofs_topological(Omega_V, gdim-1, ext_facets)
dbc_adj = [dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh, 0.0), ext_dofs, Omega_V)]

adj_problem = dolfinx.fem.petsc.LinearProblem(ufl.replace(dRdu_adj, {V_sol: Vt}), -dJdu, u=lmbda, bcs=dbc_adj, petsc_options=petsc_options)
adj_problem.solve()

dJdV1 = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(ufl.derivative(J, V1)))
dJdV1.assemble()

dJdV2 = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(ufl.derivative(J, V2)))
dJdV2.assemble()

dRdV1_form = dolfinx.fem.form(ufl.adjoint(ufl.derivative(R, V1)))
dRdV1 = dolfinx.fem.petsc.assemble_matrix(dRdV1_form)  # partial derivative
dRdV1.assemble()

dRdV2_form = dolfinx.fem.form(ufl.adjoint(ufl.derivative(R, V2)))
dRdV2 = dolfinx.fem.petsc.assemble_matrix(dRdV2_form)  # partial derivative
dRdV2.assemble()

dJdV1 += dRdV1*lmbda.x.petsc_vec
dJdV1.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
print("sensitivity of V1 from adjoint method:", dJdV1.array)

dJdV2 += dRdV2*lmbda.x.petsc_vec
dJdV2.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
print("sensitivity of V2 from adjoint method:", dJdV2.array)

# cross-validation by finite difference
eps = 1e-6
V1.x.array[:] = 1+eps
V2.x.array[:] = -1
fwd_problem.solve()

J1 = dolfinx.fem.assemble_scalar(dolfinx.fem.form(J))
print("Energy after perturbing V1: {}".format(J1))
print("Sensitivity of V1 from finite difference: {}".format((J1-J0)/eps))

V1.x.array[:] = 1
V2.x.array[:] = -1+eps
fwd_problem.solve()

J2 = dolfinx.fem.assemble_scalar(dolfinx.fem.form(J))
print("Energy after perturbing V2: {}".format(J2))
print("Sensitivity of V2 from finite difference: {}".format((J2-J0)/eps))
5

Let’s think about this from the strong point of view.
If your equation is:

\min_g J(u, u(g)) = \int_\Omega ...~\mathrm{d}x

such that

\begin{align} -\Delta u &= f && \text{in }\Omega\\ - \frac{\partial u}{\partial n} &= h &&\text{on} \partial\Omega_N\\ u &= g &&\text{on } \partial \Omega_D \end{align}

you will observe when setting up the adjoint equation that the Neumann conditions becomes homogenous (and the same for the Dirichlet condition as you differentiate the residual with respect to u).

It is not clear to me why you decided to enforce the Dirichlet conditions strongly for the adjoint when you chose to use weak enforcement for the forward problem.
Just removing the enforcement

and use dbc_adj = []
should resolve your issue.