I was able to generate more information about the issue (and this may be a useful step for others who get the dreaded Error Code 76) by implementing the custom Newton solver from the FEniCSx tutorial but with the KSP setting suggested above.
NewtonSolverCustom.py:
## setting up Custom Newton's solvers in FEniCSx
import numpy as np
import gmsh
import ufl
import basix
from mpi4py import MPI
from dolfinx import mesh, fem, io
import dolfinx
import dolfinx.fem.petsc as petsc
from dolfinx.nls.petsc import NewtonSolver
from petsc4py import PETSc
class NewtonSolverCustom():
def __init__(self, domain, F, u, bcs, params = {}):
# set up the class and create all of the things you need to do the solving
self.domain = domain
self.F = F
self.u = u
self.bcs = bcs
self.V = self.u.function_space
self.residual = dolfinx.fem.form(self.F) # construct the residual form for the problem
self.J = ufl.derivative(self.F, self.u) # calculate the jacobian of the problem
self.jacobian = dolfinx.fem.form(self.J)
# set up elements of the linear equations to solve
self.A = dolfinx.fem.petsc.create_matrix(self.jacobian)
self.L = dolfinx.fem.petsc.create_vector(self.residual)
# create the linear algebra backend to solve the system
self.solver = PETSc.KSP().create(self.domain.comm)
self.solver.setOperators(self.A)
self.solver.setType("preonly")
self.solver.setErrorIfNotConverged(True)
self.solver.getPC().setType("lu")
self.solver.getPC().setFactorSolverType("mumps")
self.du = dolfinx.fem.Function(self.V) # step to take in the Newton Solver
# parameters of the solver
self.max_iterations = 1000
self.atol = 1e-8
self.verbose = False
def solve(self,u):
i = 0
while i < self.max_iterations:
# Assemble Jacobian and residual
with self.L.localForm() as loc_L:
loc_L.set(0)
self.A.zeroEntries()
dolfinx.fem.petsc.assemble_matrix(self.A, self.jacobian, bcs=self.bcs)
self.A.assemble()
dolfinx.fem.petsc.assemble_vector(self.L, self.residual)
self.L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
self.L.scale(-1)
# Compute b - J(u_D-u_(i-1))
dolfinx.fem.petsc.apply_lifting(self.L, [self.jacobian], [self.bcs], x0=[self.u.x.petsc_vec], alpha=1)
# Set du|_bc = u_{i-1}-u_D
dolfinx.fem.petsc.set_bc(self.L, self.bcs, self.u.x.petsc_vec, 1.0)
self.L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
# Solve linear problem
self.solver.solve(self.L, self.du.x.petsc_vec)
self.du.x.scatter_forward()
# Update u_{i+1} = u_i + delta u_i
u.x.array[:] += self.du.x.array
i += 1
# Compute norm of update
correction_norm = self.du.x.petsc_vec.norm(0)
if self.verbose:
print(f"Iteration {i}: Correction norm {correction_norm}")
if correction_norm < self.atol:
converged = 1
break
if i>=self.max_iterations:
converged = 0
print("Max iterations hit. System did not converge.")
return converged
With this, the code looks like the following:
import numpy as np
import gmsh
import ufl
import basix
from mpi4py import MPI
from dolfinx import mesh, fem, io
import dolfinx
import dolfinx.fem.petsc as petsc
from dolfinx.nls.petsc import NewtonSolver
from NewtonSolverCustom import *
#################################################################
# CONSTRUCT BEAM MESH #
#################################################################
# Initialize gmsh
gmsh.initialize()
# Parameters
degree = 3 # polynomial order of the functions
N = 20*(2**2) # number of elements in the beam
L = 300e-3
char_len = L/(N-1)
d = 1
# Creating points and lines
p0 = gmsh.model.occ.addPoint(0, 0, 0, d)
p1 = gmsh.model.occ.addPoint(L, 0, 0, d)
line_ = gmsh.model.occ.addLine(p0, p1)
gmsh.model.occ.synchronize()
lines = gmsh.model.get_entities(1)
gmsh.model.add_physical_group(1, [l[1] for l in lines], 1)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",char_len)
gmsh.model.mesh.generate(dim=1)
mesh_data = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0)
try:
domain, _, _ = mesh_data
except ValueError:
domain = mesh_data.mesh
gmsh.finalize()
gdim = domain.geometry.dim
tdim = domain.topology.dim
print(f"Geometrical dimension = {gdim}")
print(f"Topological dimension = {tdim}")
#################################################################
# CONSTRUCT REFERENCE CONFIGURATION #
#################################################################
dx_dX = ufl.Jacobian(domain)[:, 0]
t = dx_dX
print(type(t))
t /= ufl.sqrt(ufl.dot(t,t)) # uncomment to normalize tangent vector
ez = ufl.as_vector([0, 0, 1])
a1 = ufl.cross(ez, t)
a1 /= ufl.sqrt(ufl.dot(a1, a1))
a2 = ufl.cross(t, a1)
a2 /= ufl.sqrt(ufl.dot(a2, a2))
def dds(u):
# use this to get derivatives along the rest configuration of the beam (d/ds)
return ufl.dot(ufl.grad(u), t)
#################################################################
# DEFINE PARAMETERS OF THE BEAM #
#################################################################
thick = fem.Constant(domain, 0.78e-3)
width = fem.Constant(domain, 30.4e-3)
E = fem.Constant(domain, 200e9)
nu = fem.Constant(domain, 0.3)
G = E/2/(1+nu)
rho = fem.Constant(domain, 2.7e-3)*1.
g = fem.Constant(domain, 9.81)
S = thick*width
ES = E*S
EI1 = E*width*thick**3/12
EI2 = E*width**3*thick/12
GJ = G*0.26*thick*width**3
kappa = fem.Constant(domain, 5./6.)
GS1 = kappa*G*S
GS2 = kappa*G*S
#################################################################
# SET UP ELEMENTS AND FUNCTION SPACES #
#################################################################
Ue = basix.ufl.element("P", domain.basix_cell(), degree, shape=(gdim,)) # here "P" means Lagrange polynomial elements (here cubic)
Us = basix.ufl.element("P", domain.basix_cell(), 1)
W = fem.functionspace(domain, basix.ufl.mixed_element([Ue, Ue]))
dstate = ufl.TestFunction(W)
dr, dd = ufl.split(dstate)
state = fem.Function(W)
r,d1 = ufl.split(state) # r: position vector, d: normal director
d2 = ufl.cross(dds(r),d1) # third director of the beam (bi-normal)
# setting up initial condition
def r0_func(x):
return x
def d0_func(x):
values = np.zeros((3,x.shape[1]))
values[1,:] = 1.0
return values
state.sub(0).interpolate(r0_func)
state.sub(1).interpolate(d0_func)
state.x.scatter_forward()
#################################################################
# SET BOUNDARY CONDITIONS AND INTEGRATION MEASURES #
#################################################################
def left(x):
return np.isclose(x[0], 0)
def right(x):
return np.isclose(x[0], L)
fdim = domain.topology.dim -1
right_facets = mesh.locate_entities_boundary(domain,fdim,right)
marked_facets = np.hstack([right_facets])
marked_values = np.hstack([np.full(len(right_facets),1,dtype=np.int32)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain,fdim,marked_facets[sorted_facets],marked_values[sorted_facets])
metadata={"quadrature_scheme":"default", "quadrature_degree": 4*degree}
ds = ufl.Measure('ds',domain=domain,subdomain_data=facet_tag,metadata=metadata)
dx = ufl.Measure('dx',domain=domain,subdomain_data=facet_tag,metadata=metadata)
Vu, Vu_to_W = W.sub(0).collapse()
Vt, Vt_to_W = W.sub(1).collapse()
Veps = fem.functionspace(domain,Us)
u_dofs = fem.locate_dofs_geometrical((W.sub(0), Vu), left)
theta_dofs = fem.locate_dofs_geometrical((W.sub(1), Vt), left)
r0 = fem.Function(Vu)
r0.interpolate(r0_func)
d0 = fem.Function(Vt)
d0.interpolate(d0_func)
# Boundary conditions: no displacement, no rotation
bcs = [fem.dirichletbc(r0, u_dofs, W.sub(0)), fem.dirichletbc(d0, theta_dofs, W.sub(1))]
#################################################################
# DEFINING THE STRAIN METRICS #
#################################################################
def norm_stretch(r,d1,d2):
eps = ufl.dot(dds(r), dds(r))
return eps
def curvature(r,d1,d2):
kap1 = ufl.dot(dds(r), dds(d1))
kap2 = ufl.dot(dds(r), dds(d2))
return kap1,kap2
def shear_strain(r,d1,d2):
gam1 = ufl.dot(dds(r), d1)
gam2 = ufl.dot(dds(r), d2)
return gam1,gam2
def shear_couples(r,d1,d2):
delt1 = ufl.dot(dds(d1),d1)
delt2 = ufl.dot(dds(d2),d2)
return delt1, delt2
def thickness_stretch(r,d1,d2):
lam1 = ufl.dot(d1,d1)
lam2 = ufl.dot(d2,d2)
return lam1, lam2
#################################################################
# DEFINING THE STRAIN ENERGY #
#################################################################
def psi(r,d1,d2):
psi = 0
eps = norm_stretch(r,d1,d2)
kap1,kap2 = curvature(r,d1,d2)
gam1,gam2 = shear_strain(r,d1,d2)
delt1,delt2 = shear_couples(r,d1,d2)
lam1,lam2 = thickness_stretch(r,d1,d2)
# axial stress contributions
psi += (ES / (2*(1-nu**2)))*( (nu/4)*(eps - 1)**2 )
# shear stress constributions
psi += (GS1/2)*gam1**2 + (GS2/2)*gam2**2
# primary bending moments contributions
psi += (EI1/2)*kap1**2 + (EI2/2)*kap2**2
# shear couple contributions
psi += 0
# CSA change contributions
psi += (ES/2)*( (lam1-1)**2 + (lam2-1)**2 )
return psi
# total potential energy
PSI = psi(r,d1,d2) * ufl.dx
# defining the work
T = fem.Constant(domain, dolfinx.default_scalar_type((0, 0, 0)))
M = fem.Constant(domain, dolfinx.default_scalar_type((0, 0, 0)))
grav = fem.Constant(domain, dolfinx.default_scalar_type((0, 1.54, 0)))
Work = ufl.dot(T,dr)*ds(1) + ufl.dot(M,dd)*ds(1) + ufl.dot(grav,dr)*ufl.dx
# Final weak form
residual = ufl.derivative(PSI,state,dstate) + Work
#################################################################
# CALCULATING VALUES FOR DEBUGGING #
#################################################################
drds = fem.Function(Vu)
try:
Vup = Vu.element.interpolation_points()
except TypeError:
Vup = Vu.element.interpolation_points
drds.interpolate(fem.Expression(dds(r),Vup))
a = fem.Function(Veps)
try:
Vepsp = Veps.element.interpolation_points()
except TypeError:
Vepsp = Veps.element.interpolation_points
a.interpolate(fem.Expression(norm_stretch(r,d1,d2) ,Vepsp))
psi_ = fem.assemble_scalar(fem.form(PSI))
print("Debugging: Are there any NaNs?")
print("tangent vector:",drds.x.array)
print("normal stretch:",a.x.array)
print("Strain Energy:",psi_)
#################################################################
# SETTING UP THE SOLVER AND ITERATING #
use_custom_solver = 1
if use_custom_solver:
solver = NewtonSolverCustom(domain, residual, state, bcs)
else:
problem = petsc.NonlinearProblem(residual,state,bcs)
solver = NewtonSolver(domain.comm, problem)
# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
ksp = solver.krylov_solver
ksp.setType("preonly")
ksp.setErrorIfNotConverged(True)
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
r_sol = state.sub(0).collapse()
d_sol = state.sub(1).collapse()
vtx = dolfinx.io.VTXWriter(domain.comm, "output.bp", [r_sol])
steps = 10
# incrementally apply a vertical load
for i in range(steps):
# update the value of T
T.value[1] = 3.94*i/(steps-1)
solver.solve(state)
state.x.scatter_forward()
r_sol.x.array[:] = state.x.array[Vu_to_W]
vtx.write(i)
vtx.close()
plot_pretty_picture = 1
if plot_pretty_picture:
from dolfinx import geometry
bb_tree = geometry.bb_tree(domain, domain.topology.dim)
points = np.zeros((3,N))
points[0] = np.linspace(0,L,N)
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions_points(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points.T)
for i, point in enumerate(points.T):
if len(colliding_cells.links(i)) > 0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
r_values = r_sol.eval(points_on_proc, cells)
import matplotlib.pyplot as plt
plt.figure(facecolor='black')
plt.rcParams['font.size'] = 10 # Sets the default font size to 12
ax = plt.gca()
ax.set(facecolor='black')
# set various colors
ax.spines['bottom'].set_color('white')
ax.spines['top'].set_color('white')
ax.spines['right'].set_color('white')
ax.spines['left'].set_color('white')
ax.xaxis.label.set_color('white')
ax.yaxis.label.set_color('white')
ax.tick_params(colors='white', which='both') # 'both' refers to minor and major axes
plt.plot(points[0,:],points[1,:],'--',color=(0.8,0.4,0.8))
plt.plot(r_values[:,0],r_values[:,1],'-m')
plt.axis("equal")
plt.title("N="+str(N),color='white')
plt.show()
I have confirmed that this solver works correctly and provides the same linear element convergence as the built-in solver.
When I try to run the denser meshes for quadratic or cubic functions, instead of just getting
Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 76, Error in external library
I get more specific error from MUMPS. Specifically, I either get
MUMPS error in numerical factorization: INFOG(1)=-9, INFO(2)=404
(reproduced by setting N=20*(2**4)) or
MUMPS error in numerical factorization: INFOG(1)=-10, INFO(2)=11506
(reproduced by setting N=20*(2**4)). According to the MUMPS docs, the first (INFOG(1)=-9) corresponds to “The main internal real/complex workarray S is too small.” The second (INFOG(1)=-10) corresponds to “Numerically singular matrix, or zero pivot encountered.”
Do you have suggestions for how to sniff out the singularity (if it is that) or the zero pivot? Given that the linear elements and less dense cubic elements provide the correct answer, I do not believe this is a weak form issue, but likely arising from the elements side of things. Any suggestions would be greatly appreciated! (I am also working on the plate extension of this model, and for the sake of convergence rates, I would really like to use higher order elements, but I would first like to confirm that I can demonstrate proper convergence first.)