Hello,
I am working on a finite element project where I solve the Stokes equations with Navier slip boundary conditions using FEniCSx.
However, I am encountering a PETSc segmentation fault that seems related to memory handling. The strange part is that the error only appears for some mesh sizes, and the behavior changes when I add print() statements.
To investigate, I created a loop to iterate on several mesh sizes :
from mpi4py import MPI
import numpy as np
import time
# ----------------------------
# Spatial discretization
# ----------------------------
from dolfinx import mesh
for N in [10, 15, 20, 25, 26, 28, 35, 55, 60]:
domain = mesh.create_unit_square(MPI.COMM_WORLD, N, N, mesh.CellType.triangle)
# --------------------------------------------------
# Function spaces (Taylor-Hood)
# --------------------------------------------------
from basix.ufl import element
from ufl import MixedFunctionSpace
from dolfinx import(
fem,
default_scalar_type,
default_real_type,
)
V_el = element("Lagrange", domain.basix_cell(), 2, shape=(domain.geometry.dim,), dtype=default_real_type)
Q_el = element("Lagrange", domain.basix_cell(), 1, dtype=default_real_type)
V = fem.functionspace(domain, V_el)
Q = fem.functionspace(domain, Q_el)
W = MixedFunctionSpace(V, Q)
import ufl
from ufl import(
as_vector,
sym,
grad,
div,
SpatialCoordinate
)
def _f(x):
return as_vector((x[0], x[0]))
x = SpatialCoordinate(domain)
f = _f(x)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
left = mesh.locate_entities_boundary(
domain, fdim,
lambda x: np.isclose(x[0], 0.0)
)
right = mesh.locate_entities_boundary(
domain, fdim,
lambda x: np.isclose(x[0], 1.0)
)
top = mesh.locate_entities_boundary(
domain, fdim,
lambda x: np.isclose(x[1], 1.0)
)
bottom = mesh.locate_entities_boundary(
domain, fdim,
lambda x: np.isclose(x[1], 0.0)
)
left_dofs = fem.locate_dofs_topological(V.sub(0), fdim, left)
right_dofs = fem.locate_dofs_topological(V.sub(0), fdim, right)
top_dofs = fem.locate_dofs_topological(V.sub(1), fdim, top)
bottom_dofs = fem.locate_dofs_topological(V.sub(1), fdim, bottom)
left_right_dofs = np.concatenate([left_dofs, right_dofs])
top_bottom_dofs = np.concatenate([top_dofs, bottom_dofs])
bcx = fem.dirichletbc(default_scalar_type(0), left_right_dofs, V.sub(0))
bcy = fem.dirichletbc(default_scalar_type(0), top_bottom_dofs, V.sub(1))
bc_p = fem.dirichletbc(default_scalar_type(0), left_dofs[0:1], Q)
from ufl import(
TrialFunctions,
TestFunctions,
inner,
dx,
ds
)
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
zero = fem.Constant(domain, default_scalar_type(0))
def D(w):
return sym(grad(w))
_alpha = 1.0
_nu = 1.0
alpha = fem.Constant(domain, _alpha)
nu = fem.Constant(domain, _nu)
a = (
2 * nu * inner(D(u), D(v)) * dx
+ alpha * inner(u, v) * ds
- div(v) * p * dx
- div(u) * q * dx
+ zero * p * q * dx
)
L = inner(f, v) * dx + zero * q * dx
from dolfinx.fem.petsc import LinearProblem
petsc_options = {
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"ksp_error_if_not_converged": True,
"ksp_monitor": None,
}
problem = LinearProblem(ufl.extract_blocks(a), ufl.extract_blocks(L), bcs=[bcx, bcy, bc_p], petsc_options=petsc_options, petsc_options_prefix="NS")
time.sleep(2)
uh, ph = problem.solve()
time.sleep(2)
The program runs almost until the end of the loop, then crashes with:
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
If I add print(N) after problem.solve(), the crash happens earlier (e.g., around N=30).
If I add print(N) at the beginning of the loop, the program finishes without error.
I do not understand why memory errors appears and how adding print statements, which should not affect the solver, changes whether the program crashes.
P.S. It also happen that when I run the program for a fixed value for N (without any loop), the program sometimes crashes during solving the problem, sometimes after solving the problem and sometimes does not crash.
I use the dolfinx 0.10.0 on windows, and I use WSL as well.
Any insight would be greatly appreciated !