I am encountering an error when running this minimal working example (see below) of a nonlinear mixed FEniCSx problem inside an HPC environment using Apptainer. I built the FEniCSx image inside the container, but I am unsure whether the issue is related to PETSc configuration, MPI, or how dolfinx was built inside the container.
Issue : I repeatedly encounter Runtime Error PETSc code 63, and AttributeError of ‘NonlinearProblem’ object having no attribute ‘_snes’. This is when trying to initialize the nonlinear problem instance in the code line
problem = NonlinearProblem(F, U, J=fem.form(J), bcs=[bc], petsc_options_prefix=petsc_options)
The FEniCSx configuration I have inside the container is :
Python 3.12.3
/dolfinx-env/bin/python
fenics-basix 0.10.0
fenics-dolfinx 0.10.0
fenics-ufl 2025.2.0
mpi4py 4.1.1
petsc4py 3.24.0
Here is the MWE :
from petsc4py import PETSc
from mpi4py import MPI
import ufl
import numpy as np
from dolfinx import mesh, fem, io, nls, log, io
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import os
import basix
a = 10 # domain size
b = 6 # height of initial gaussian
c = 5e-2 # stdev of initial gaussian
kappa = 0.1 # diffusion coefficient
Pe = 0.5 # Desired Peclet Number (to compute mesh size, always <1)
cfl = 0.5
alpha = 1e-3 # reaction term
q = 1 # constant in the beta-term, equal to q/epsilon where q=1
epsilon = 0.3 # constant in the reaction term
u_pc = 3 # limit in the step-function
t = 0 # Start time
T = 30 # End time
max_wind = 1
nx, ny = 200,200
h = (a/200)
dt = 0.1
num_steps = np.int64(T//dt)
# initialize domain
bottom_left = np.array([-a,-a])
top_right = np.array([a,a])
domain = mesh.create_rectangle(MPI.COMM_WORLD, [bottom_left, top_right], [nx, nx], mesh.CellType.triangle)
el_u = basix.ufl.element("Lagrange", domain.basix_cell(), 1)
el_beta = basix.ufl.element("Lagrange", domain.basix_cell(), 1)
el_mixed = basix.ufl.mixed_element([el_u, el_u])
mixed_space = fem.functionspace(domain, el_mixed)
# set homogenous dirichlet bcs
w_D = fem.Function(mixed_space)
w_D.x.array[:] = 0.0
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
wh = fem.Function(mixed_space)
dofs = fem.locate_dofs_topological(mixed_space, fdim, boundary_facets)
bc = fem.dirichletbc(w_D, dofs)
bc.set(wh.x.array)
wh.x.scatter_forward()
# setup variational form
# Test functions and unknown
V0 = mixed_space.sub(0).collapse()[0]
V1 = mixed_space.sub(1).collapse()[0]
v_, w_ = ufl.TestFunctions(mixed_space)
U = fem.Function(mixed_space) # unknown solution
u_, beta_ = ufl.split(U)
U_prev = fem.Function(mixed_space)
a, b, c = 20, 6, 5e-2
def ic_func_u(x):
y = 6 * np.exp( -0.05* ((x[0]+9)**2 + (x[1]+9)**2))
return y
def ic_func_beta(x):
# x has shape (gdim, npoints)
npoints = x.shape[1]
return np.random.uniform(0.0, 1.0, size=npoints)
U_prev.sub(0).interpolate(ic_func_u)
U_prev.sub(1).interpolate(ic_func_beta)
u_prev_fn, beta_prev_fn = U_prev.split()
u_pc_const = fem.Constant(domain, PETSc.ScalarType(u_pc))
delta_const = fem.Constant(domain, PETSc.ScalarType(0.05))
t = (u_ - (u_pc - delta_const)) / (2.0 * delta_const)
# polynomial ramp
ramp = 3*t**2 - 2*t**3
H = ufl.conditional(ufl.lt(u_, u_pc - delta_const), 0.0, ufl.conditional(ufl.gt(u_, u_pc + delta_const), 1.0, ramp))
dt_const = fem.Constant(domain, PETSc.ScalarType(dt))
epsilon_const = fem.Constant(domain, PETSc.ScalarType(epsilon))
q_const = fem.Constant(domain, PETSc.ScalarType(q))
kappa_const = fem.Constant(domain, PETSc.ScalarType(kappa))
alpha_const = fem.Constant(domain, PETSc.ScalarType(alpha))
u_term = u_ / (1 + epsilon_const * u_)
F1 = ((-H * (epsilon_const / q_const) * ufl.exp(u_term) * beta_) * dt_const * v_
- beta_ * v_
+ beta_prev_fn * v_) * ufl.dx
val_c = fem.Constant(domain, PETSc.ScalarType(float(np.sqrt(2)/(2))))
wind_vec = ufl.as_vector([val_c, val_c])
F2 = ((kappa_const * ufl.dot(ufl.grad(u_), ufl.grad(w_))
+ w_ * ufl.dot(wind_vec, ufl.grad(u_))
- H * beta_ * ufl.exp(u_term) * w_
+ alpha_const * u_ * w_) * dt_const
+ u_ * w_ - u_prev_fn * w_) * ufl.dx
F = F1 + F2
u_sol = fem.Function(V0)
beta_sol = fem.Function(V1)
# setup solver
problem = NonlinearProblem(F, U, bcs=[bc])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
# configure linear solver
solver.report = False
solver.convergence_criterion = "incremental"
solver.error_on_nonconvergence = True
solver.rtol = np.sqrt(np.finfo(float).eps) * 1e-2
ksp = solver.krylov_solver
opts = PETSc.Options()
prefix = ksp.getOptionsPrefix()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{prefix}ksp_rtol"] = 1e-8
opts[f"{prefix}ksp_max_it"] = 1000
sys = PETSc.Sys()
if sys.hasExternalPackage("superlu_dist"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "superlu_dist"
elif sys.hasExternalPackage("mumps"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
solver.solve(U)