Different results when solving linear elastic problem using fem.petsc.LinearProblem vs building PETSc KSP directly

I’m using the Arnold-Falk-Winther elements in dolfinx 0.9.0 to solve the linear elastic problem and am getting a different solution when I solve using fem.petsc.LinearProblem than when I create the matrix and load vector (using fem.petsc.assemble_matrix_nest and fem.petsc.assemble_vector_nest respectively) and solve with PETSc directly. fem.petsc.LinearProblem is giving me correct results, while solving the problem in PETSc gives incorrect results. To illustrate the problem I solve a cook membrane like problem on the unit square (zero displacement on left of square, zero traction on top/bottom, and applied traction on right side) for a MWE and evaluate the displacement at the top right corner of the unit square. The “exact” solution (using a high order primal formulation on a refined mesh) should be [-0.25185959 0.59833939]. Here is the code using fem.petsc.LinearProblem:

#Import packages
from dolfinx import mesh as dmesh, fem, geometry
import dolfinx.fem.petsc
import basix
import ufl
from mpi4py import MPI
import numpy as np

#Define material parameters
E = 70                        #Young modulus
nu = 0.4                      #Poisson's ratio
lam_ = E*nu/((1+nu)*(1-2*nu)) #Lame's first parameter Pa
mu_ = E/(2*(1+nu))            #Shear modulus Pa
trac_val = 6.25               #Traction on right edge
deg_fe = 1

def in_right(x):
    return np.isclose(x[0], 1.)

def in_top_bottom(x):
    return np.isclose(np.abs(x[1]-0.5), 0.5)

def traction(x):
    values = np.zeros((2, x.shape[1]))
    values[0] = np.full(x.shape[1], trac_val)
    return values

#Create mesh
mesh = dmesh.create_unit_square(MPI.COMM_WORLD, 5, 5, cell_type=dmesh.CellType.triangle, diagonal=dmesh.DiagonalType.crossed)

#Set up AFW element
Se = basix.ufl.element('N2F', mesh.basix_cell(), deg_fe)
Ve = basix.ufl.element('DG', mesh.basix_cell(), deg_fe-1, shape=(2,))
Ke = basix.ufl.element('DG', mesh.basix_cell(), deg_fe-1)
AFW = fem.functionspace(mesh, basix.ufl.mixed_element([Se, Se, Ve, Ke]))

#Define zero traction (essential) boundary conditions on top and bottom
Sigma_h = [AFW.sub(i).collapse()[0] for i in range(2)]
top_bottom_facets = dmesh.locate_entities_boundary(mesh, 1, in_top_bottom)
top_bottom_dofs = [fem.locate_dofs_topological((AFW.sub(i), Sigma_h[i]), 1, top_bottom_facets) for i in range(2)]
sig_zeros = [fem.Function(Sigma_h[i]) for i in range(2)]
for sig in sig_zeros:
    sig.x.array[:] = 0.
bcs_top_bottom = [fem.dirichletbc(sig_zeros[i], top_bottom_dofs[i], AFW.sub(i)) for i in range(2)]

#Define applied traction (essential) boundary on right
right_facets = dmesh.locate_entities_boundary(mesh, 1, in_right)
right_dofs = [fem.locate_dofs_topological((AFW.sub(i), Sigma_h[i]), 1, right_facets) for i in range(2)]
tbar = fem.Function(Sigma_h[1])
tbar.interpolate(traction)
bcs_right = [fem.dirichletbc(sig_zeros[0], right_dofs[0], AFW.sub(0)), fem.dirichletbc(tbar, right_dofs[1], AFW.sub(1))]
bcs = bcs_top_bottom + bcs_right

#Build forms
(sig0, sig1, u, zeta) = ufl.TrialFunctions(AFW)
(tau0, tau1, v, eta) = ufl.TestFunctions(AFW)
a_form = 1./(2*mu_)*(ufl.inner(sig0, tau0) + ufl.inner(sig1, tau1))*ufl.dx - lam_/((4*mu_)*(mu_+lam_))*(sig0[0]+sig1[1])*(tau0[0]+tau1[1])*ufl.dx
b_form = (ufl.div(sig0)*v[0] + ufl.div(sig1)*v[1])*ufl.dx
bT_form = (ufl.div(tau0)*u[0] + ufl.div(tau1)*u[1])*ufl.dx
c_form = (sig0[1] - sig1[0])*eta*ufl.dx
cT_form = (tau0[1] - tau1[0])*zeta*ufl.dx
lhs = a_form + b_form + bT_form + c_form + cT_form
rhs = ufl.inner(fem.Constant(mesh, (0., 0.)), v)*ufl.dx

#Solve
problem = fem.petsc.LinearProblem(lhs, rhs, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
w = problem.solve()
w.x.scatter_forward()
(sig0h, sig1h, uh, _) = w.split()

#Evaluate displacement at top right corner
p = np.array([[1., 1., 0.]], dtype=np.float64)
bb_tree = geometry.bb_tree(mesh, 2)
cell_candidates = geometry.compute_collisions_points(bb_tree, p)
cells = geometry.compute_colliding_cells(mesh, cell_candidates, p).array
print(f'Displacement at (1., 1.) is {uh.eval(p, cells[0])}')

The output of the above code is “Displacement at (1., 1.) is [-0.21290887 0.5055914]” which is close to the correct answer. The following is the exact same problem, the only difference is that I solve using PETSc directly:

#Import packages
from dolfinx import mesh as dmesh, fem, la, geometry
import dolfinx.fem.petsc
from petsc4py import PETSc
import basix
import ufl
from mpi4py import MPI
import numpy as np

#Define material parameters
E = 70                        #Young modulus
nu = 0.4                      #Poisson's ratio
lam_ = E*nu/((1+nu)*(1-2*nu)) #Lame's first parameter Pa
mu_ = E/(2*(1+nu))            #Shear modulus Pa
trac_val = 6.25               #Traction on right edge
deg_fe = 1

def in_right(x):
    return np.isclose(x[0], 1.)

def in_top_bottom(x):
    return np.isclose(np.abs(x[1]-0.5), 0.5)

def traction(x):
    values = np.zeros((2, x.shape[1]))
    values[0] = np.full(x.shape[1], trac_val)
    return values

#Create mesh
mesh = dmesh.create_unit_square(MPI.COMM_WORLD, 5, 5, cell_type=dmesh.CellType.triangle, diagonal=dmesh.DiagonalType.crossed)

#Set up AFW element
Se = basix.ufl.element('N2F', mesh.basix_cell(), deg_fe)
Ve = basix.ufl.element('DG', mesh.basix_cell(), deg_fe-1, shape=(2,))
Ke = basix.ufl.element('DG', mesh.basix_cell(), deg_fe-1)
S = fem.functionspace(mesh, basix.ufl.mixed_element([Se, Se]))
V = fem.functionspace(mesh, Ve)
K = fem.functionspace(mesh, Ke)

#Define zero traction boundary conditions on top and bottom
Sigma_h = [S.sub(i).collapse()[0] for i in range(2)]
top_bottom_facets = dmesh.locate_entities_boundary(mesh, 1, in_top_bottom)
top_bottom_dofs = [fem.locate_dofs_topological((S.sub(i), Sigma_h[i]), 1, top_bottom_facets) for i in range(2)]
sig_zeros = [fem.Function(Sigma_h[i]) for i in range(2)]
for sig in sig_zeros:
    sig.x.array[:] = 0.
bcs_top_bottom = [fem.dirichletbc(sig_zeros[i], top_bottom_dofs[i], S.sub(i)) for i in range(2)]

#Define applied traction boundary on right
right_facets = dmesh.locate_entities_boundary(mesh, 1, in_right)
right_dofs = [fem.locate_dofs_topological((S.sub(i), Sigma_h[i]), 1, right_facets) for i in range(2)]
tbar = fem.Function(Sigma_h[1])
tbar.interpolate(traction)
bcs_right = [fem.dirichletbc(sig_zeros[0], right_dofs[0], S.sub(0)), fem.dirichletbc(tbar, right_dofs[1], S.sub(1))]
bcs = bcs_top_bottom + bcs_right
bcs_by_block = [bcs, [], []]

#Build forms
(sig0, sig1), (tau0, tau1) = ufl.TrialFunctions(S), ufl.TestFunctions(S)
(u, v) = ufl.TrialFunction(V), ufl.TestFunction(V)
(zeta, eta) = ufl.TrialFunction(K), ufl.TestFunction(K)
a_form = 1./(2*mu_)*(ufl.inner(sig0, tau0) + ufl.inner(sig1, tau1))*ufl.dx - lam_/((4*mu_)*(mu_+lam_))*(sig0[0]+sig1[1])*(tau0[0]+tau1[1])*ufl.dx
b_form = (ufl.div(sig0)*v[0] + ufl.div(sig1)*v[1])*ufl.dx
bT_form = (ufl.div(tau0)*u[0] + ufl.div(tau1)*u[1])*ufl.dx
c_form = (sig0[1] - sig1[0])*eta*ufl.dx
cT_form = (tau0[1] - tau1[0])*zeta*ufl.dx
tau_zero = ufl.inner(fem.Constant(mesh, (0., 0.)), tau0)*ufl.dx
v_zero = ufl.inner(fem.Constant(mesh, (0., 0.)), v)*ufl.dx
eta_zero = fem.Constant(mesh, 0.)*eta*ufl.dx

#Create blocked forms
lhs_form = fem.form([[a_form,  bT_form, cT_form], \
                     [b_form, None, None], \
                     [c_form, None, None]])
rhs_form = fem.form([tau_zero, v_zero, eta_zero])

#Assemble matrices
A = fem.petsc.assemble_matrix_nest(lhs_form, bcs=bcs)
A.assemble()
b = fem.petsc.assemble_vector_nest(rhs_form)
fem.petsc.apply_lifting_nest(b, lhs_form, bcs=bcs)
for b_sub in b.getNestSubVecs():
    b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc_nest(b, bcs=bcs_by_block)

#Solve
ksp = PETSc.KSP()
ksp.create(comm=mesh.comm)
ksp.setOperators(A)
pc = ksp.getPC()
ksp.setType(PETSc.KSP.Type.PREONLY)
pc.setType(PETSc.PC.Type.LU)
pc.setFactorSolverType(PETSc.Mat.SolverType.MUMPS)
sigh, uh, zetah = fem.Function(S), fem.Function(V), fem.Function(K)
x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(sigh.x), la.create_petsc_vector_wrap(uh.x), la.create_petsc_vector_wrap(zetah.x)])
ksp.solve(b, x)
sigh.x.scatter_forward(), uh.x.scatter_forward(), zetah.x.scatter_forward()

#Evaluate
p = np.array([[1., 1., 0.]], dtype=np.float64)
bb_tree = geometry.bb_tree(mesh, 2)
cell_candidates = geometry.compute_collisions_points(bb_tree, p)
cells = geometry.compute_colliding_cells(mesh, cell_candidates, p).array
print(f'Displacement at (1., 1.) is {uh.eval(p, cells[0])}')

The output of the above is “Displacement at (1., 1.) is [-0.10308539 0.19198263]” which is not the same. My initial though was that the essential boundary conditions were being applied incorrectly and so I integrated the stress tensor along the top/bottom boundary and right boundary and get the correct results. Does anyone see what I’m doing incorrect? I’d appreciate any help.

Here is the code to integrate the rows of the stress tensor:

#Create meshtags for boundary integrals
nfacet = mesh.topology.index_map(1).size_local + mesh.topology.index_map(1).num_ghosts
local_facets = np.arange(nfacet, dtype=np.int32)
right_facets = dmesh.locate_entities_boundary(mesh, 1, in_right)
top_bottom_facets = dmesh.locate_entities_boundary(mesh, 1, in_top_bottom)
values = np.full(nfacet, -1, dtype=np.int32)
values[top_bottom_facets] = 1
values[right_facets] = 2
mt = dmesh.meshtags(mesh, 1, local_facets, values)

#Check if solution satisfies essential BC
(sig0h, sig1h) = sigh.split()
n = ufl.FacetNormal(mesh)
ds = ufl.Measure('ds', domain=mesh, subdomain_data=mt)
print('sig0 dot n on top/bottom', fem.assemble_scalar(fem.form(ufl.inner(sig0h, n)*ds(1)))/fem.assemble_scalar(fem.form(fem.Constant(mesh, 1.)*ds(1))))
print('sig1 dot n on top/bottom', fem.assemble_scalar(fem.form(ufl.inner(sig1h, n)*ds(1)))/fem.assemble_scalar(fem.form(fem.Constant(mesh, 1.)*ds(1))))
print('sig0 dot n on rhs', fem.assemble_scalar(fem.form(ufl.inner(sig0h, n)*ds(2)))/fem.assemble_scalar(fem.form(fem.Constant(mesh, 1.)*ds(2))))
print('sig1 dot n on rhs', fem.assemble_scalar(fem.form(ufl.inner(sig1h, n)*ds(2)))/fem.assemble_scalar(fem.form(fem.Constant(mesh, 1.)*ds(2))))