Hi everyone,

I have a problem when using PETSc’s SNES directly for solving nonlinear problems with a direct (LU) solver for the linear system, at least for newer PETSc versions. I have created the following MWE, showing the problem

```
from fenics import *
from petsc4py import PETSc
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
v = TestFunction(V)
F = dot(grad(u), grad(v)) * dx - Constant(1.0) * v * dx
def boundary(x, on_boundary):
return on_boundary
bcs = DirichletBC(V, Constant(0.0), boundary)
du = TrialFunction(V)
dF = derivative(F, u, du)
def ksp_solve():
A = PETScMatrix()
b = PETScVector()
assemble(dF, tensor=A)
assemble(-F, tensor=b)
bcs.apply(A)
bcs.apply(b)
A = as_backend_type(A).mat()
b = as_backend_type(b).vec()
opts = PETSc.Options()
opts["ksp_type"] = "cg"
opts["pc_type"] = "lu"
opts["ksp_monitor_true_residual"] = None
ksp = PETSc.KSP().create()
ksp.setOperators(A)
ksp.setFromOptions()
x = u.vector().vec()
ksp.solve(b, x)
ksp.destroy()
def snes_solve():
def assemble_function(snes, x, f):
u.vector().vec().setArray(x)
u.vector().apply("")
b = PETScVector(f)
assemble(F, tensor=b)
bcs.apply(b)
def assemble_jacobian(snes, x, J, P):
u.vector().vec().setArray(x)
u.vector().apply("")
A = PETScMatrix(J)
assemble(dF, tensor=A)
bcs.apply(A)
AP = PETScMatrix(P)
assemble(dF, tensor=AP)
bcs.apply(AP)
f = PETScVector()
J = PETScMatrix()
P = PETScMatrix()
opts = PETSc.Options()
opts["snes_type"] = "newtonls"
opts["snes_monitor"] = None
opts["ksp_type"] = "cg"
opts["pc_type"] = "none" # when changing this to "lu", things break
opts["ksp_monitor_true_residual"] = None
snes = PETSc.SNES().create()
snes.setFunction(assemble_function, f.vec())
snes.setJacobian(assemble_jacobian, J.mat(), P.mat())
snes.setFromOptions()
x = u.vector().vec()
snes.solve(None, x)
snes.destroy()
# ksp_solve() # works fine
snes_solve() # does not work properly for "lu"
```

As you can see, I am using a linear model problem, which is to be solved. The problem is not present, when I am using the KSP directly, as can be seen from my example. However, when I am using the `snes_solve`

function and set

```
opts["pc_type"] = "lu"
```

I am getting the error message

```
Segmentation fault (core dumped)
```

with PETSc 3.21 and 3.20, however, the code works fine for my old PETSc version 3.17.4. I have also tested several linear solvers (petsc, umfpack, mumps), but the error is present regardless of which solver I am trying to use.

Also, I have added a “dummy” computation for the preconditioner matrix in `assemble_jacobian`

, but this does not change anything as far as I noticed.

Does anyone have an idea / solution for this strange problem?

Thanks a lot in advance,

Sebastian