I believe that the issue I am facing is related to an improper use of scatter_forward as discussed here.
I am attempting to differentiate the functional J(u) = 1/2 integral (k1 + k2 u^2) grad(u)^T grad(u) dx. Do a finite difference check for consistency that is |(J(u0 + eps uhat) - J(u0)) / eps - gradJ(u_0)^T uhat| = O(eps). This seems to be done correctly in serial and parallel as indicated by my finite difference plots. I try to do the same finite difference check but for the Hessian of J and it work in serial but not in parallel. That is I check that || (gradJ(u0 + eps * uhat) - gradJ(u0)) / eps - HessianJ(u0) * uhat || = eps.
Any suggestions would be appreciated and if someone can point me to more details on what is occuring with scatter_forward and scatter_backward operations that would be greatly appreciated.
# Import FEniCSx
import dolfinx as dl
import ufl
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import matplotlib.pyplot as plt
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
# Define the finite element mesh. The mesh size h is 1/nx
nx = 16
ny = nx
mesh = dl.mesh.create_unit_square(comm, nx, ny, dl.mesh.CellType.triangle)
# Define the finite element space V_h as the space of piecewise linear functions on the elements of the mesh.
degree = 1
Vh = dl.fem.FunctionSpace(mesh, ("CG", degree))
# define the homogeneous Dirichlet boundary condition at x = 0
facet_dim = mesh.topology.dim-1
facets_D = dl.mesh.locate_entities_boundary(mesh, dim=facet_dim, \
marker=lambda x: np.isclose(x[0], 0.0))
dofs_D = dl.fem.locate_dofs_topological(V=Vh, entity_dim=facet_dim, entities=facets_D)
u_bc = dl.fem.Constant(mesh, PETSc.ScalarType(0.0))
bcs = [dl.fem.dirichletbc(u_bc, dofs_D, Vh)]
# ----------- define the functional ------------------
# J(u) = \int_{\Omega} (0.5 * (k1 + k2 u^2) grad(u)^T grad(u) dx
uh = dl.fem.Function(Vh)
k1 = dl.fem.Constant(mesh, PETSc.ScalarType(0.05))
k2 = dl.fem.Constant(mesh, PETSc.ScalarType(1.0))
Jform = PETSc.ScalarType(0.5)*(k1 + k2*uh*uh)*ufl.inner(ufl.grad(uh), ufl.grad(uh))*ufl.dx
J = dl.fem.form(Jform)
gradform = ufl.derivative(Jform, uh)
grad = dl.fem.form(gradform)
u0 = dl.fem.Function(Vh)
u0.interpolate(lambda x: x[0]*(x[0]-1)*x[1]*(x[1]-1))
# --- evaluate J at u = u0
uh.vector.aypx(0.0, u0.vector) # uh = u0 + 0.0 * uh
J0 = mesh.comm.allreduce(dl.fem.assemble_scalar(J), op = MPI.SUM)
# evaluate the gradient of J at u = u0
grad0 = dl.fem.petsc.assemble_vector(grad)
dl.fem.petsc.set_bc(grad0, bcs)
grad0.assemble()
# set up a random vector that obeys the Dirichlet conditions
uhat = dl.la.create_petsc_vector(Vh.dofmap.index_map, Vh.dofmap.index_map_bs)
uhat.setRandom()
dl.fem.petsc.set_bc(uhat, bcs)
# grad(J)^T uhat
grad0_uhat = grad0.dot(uhat)
n_eps = 32
epss = 1e-2*np.power(2., -np.arange(n_eps))
fdgrad_residual = np.zeros(n_eps)
for i, eps in enumerate(epss):
uh.vector.aypx(0.0, u0.vector)
uh.vector.axpy(eps, uhat) # uh = uh + eps[i]*dir
Jplus = mesh.comm.allreduce(dl.fem.assemble_scalar(J), op = MPI.SUM)
fdgrad_residual[i] = abs( (Jplus - J0)/eps - grad0_uhat)
if rank == 0:
plt.loglog(epss, fdgrad_residual, "-ob", label="Error Grad")
plt.loglog(epss, (.5*fdgrad_residual[0]/epss[0])*epss, "-.k", label=r"First Order, $\propto\varepsilon$")
plt.title("Finite difference check of the first variation")
plt.xlabel(r"$\varepsilon$")
plt.ylabel(r"$r_{1}(\varepsilon)$, finite-difference error")
plt.legend(loc = "upper left")
plt.show()
Hform = ufl.derivative(gradform, uh)
H = dl.fem.form(Hform)
# ---- evaluate the Hessian of J at u0
uh.vector.aypx(0.0, u0.vector)
H_0 = dl.fem.petsc.assemble_matrix(H, bcs)
H_0.assemble()
# ----- compute H(u0) * udir
H_0udir = H_0.createVecLeft()
H_0.mult(uhat, H_0udir)
H_0udir.assemble()
# -----
fdH_residuals = np.zeros(n_eps)
diff_grad = H_0.createVecLeft()
for i, eps in enumerate(epss):
# grad(J) evaluated at u = u0 + eps * uhat
uh.vector.aypx(0.0, u0.vector)
uh.vector.axpy(eps, uhat)
grad_plus = dl.fem.petsc.assemble_vector(grad)
dl.fem.petsc.set_bc(grad_plus, bcs)
grad_plus.assemble()
# ------- diff_grad = [ dJ(u0 + eps * uhat) - dJ(u0) ] / eps - H(u0) * uhat
diff_grad.zeroEntries()
diff_grad.axpy(1., grad_plus)
diff_grad.axpy(-1., grad0)
diff_grad.scale(1./eps)
diff_grad.axpy(-1, H_0udir)
diff_grad.assemble()
# ------ || diff_grad ||_2 = O(eps)
fdH_residuals[i] = diff_grad.norm(2)
if rank == 0:
plt.figure()
plt.loglog(epss, fdH_residuals, "-ob", label="Error Hessian")
plt.loglog(epss, (.5*fdH_residuals[0]/epss[0])*epss, "-.k", label="First Order")
plt.title("Finite difference check of the second variation")
plt.xlabel(r"$\varepsilon$")
plt.ylabel(r"$r_{2}(\varepsilon)$, finite-difference error")
plt.legend(loc = "upper left")
plt.show()