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.