# Setting SNES solver with petsc4py segfaults

**URL:** <https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149>\
**Category:** Linear Algebra\
**Created:** [February 24, 2021, 3:08pm UTC](https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149 "2021-02-24T15:08:07Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\
**Post date:** [February 24, 2021, 3:08pm UTC](https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149/1 "2021-02-24T15:08:07Z")

</div>

Hi,  
can somebody confirm that the snippet below produces a SEGFAULT ? It seems that removing the `setVariableBounds` works fine. Does anyone has any idea on the origin of this problem ?

```python
from dolfin import *
from petsc4py import PETSc

N = 5
mesh = UnitSquareMesh(N, N)

V = FunctionSpace(mesh, "CG", 1)

d = Function(V)
dlb = Function(V)
dub = Function(V)
dub.interpolate(Constant(1.))

d_ = TestFunction(V)
dd = TrialFunction(V)

energy = ((d-0.5)**2)*dx
F = derivative(energy, d, d_)
J = derivative(F, d, dd)

b = PETScVector()  
A = PETScMatrix()

def function(solver, x):
    return assemble(energy)
def gradient(solver, x, b):
    b = PETScVector(b)
    return assemble(F, tensor=b)
def hessian(solver, x, A, P):
    A = PETScMatrix(A)
    return assemble(J, tensor=A)

snes_solver = PETSc.SNES().create(PETSc.COMM_WORLD)
PETScOptions.set("snes_monitor")
snes_solver.setFromOptions()
snes_solver.setType(PETSc.SNES.Type.VINEWTONSSLS)
snes_solver.setFunction(gradient, b.vec())
snes_solver.setJacobian(hessian, A.mat())
snes_solver.setVariableBounds(dlb.vector().vec(), dub.vector().vec())
snes_solver.solve(None, d.vector().vec())
print("Succeeded SNES Solve")

```

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [February 24, 2021, 3:15pm UTC](https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149/2 "2021-02-24T15:15:57Z")

</div>

I can confirm that this gives me a segmentation fault. (I’m using FEniCS on Ubuntu, installed through the package manager.)

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [February 24, 2021, 3:17pm UTC](https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149/3 "2021-02-24T15:17:08Z")

</div>

Yep, segfault for me with native build in `snes_solver.setVariableBounds`.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 24, 2021, 3:40pm UTC](https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149/4 "2021-02-24T15:40:30Z")

</div>

Maybe not the answer you would like, but this does not segfault with dolfinx:

```auto
import dolfinx
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np
import ufl
N = 5
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, N, N)

V = dolfinx.FunctionSpace(mesh, ("CG", 1))

def one(x):
    return np.ones(x.shape[1])

d = dolfinx.Function(V)
dlb = dolfinx.Function(V)
dub = dolfinx.Function(V)
dub.interpolate(one)

d_ = ufl.TestFunction(V)
dd = ufl.TrialFunction(V)

energy = ((d-0.5)**2)*ufl.dx
F = ufl.derivative(energy, d, d_)
J = ufl.derivative(F, d, dd)

b = dolfinx.fem.create_vector(F)  
A = dolfinx.fem.create_matrix(J)

def function(solver, x):
    return dolfinx.fem.assemble_scalar(energy)
def gradient(solver, x, b):
    b = b.copy()
    with b.localForm() as loc:
        loc.set(0)
    return dolfinx.fem.assemble_vector(b, F)
def hessian(solver, x, A, P):
    A = A.copy()
    A.setZeroEntries()
    return dolfinx.fem.assemble_matrix(A, J)

snes_solver = PETSc.SNES().create(MPI.COMM_WORLD)
opts = PETSc.Options()
opts["snes_monitor"] = None
snes_solver.setFromOptions()
snes_solver.setType(PETSc.SNES.Type.VINEWTONSSLS)
snes_solver.setFunction(gradient, b)
snes_solver.setJacobian(hessian, A)
snes_solver.setVariableBounds(dlb.vector, dub.vector)
snes_solver.solve(None, d.vector)
print("Succeeded SNES Solve")

```

Please let me know if there where any typos in the conversion. dolfinx uses version `3.14.0` of petsc4py/PETSc

---

<div class="post-metadata">

**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\
**Post date:** [February 24, 2021, 3:41pm UTC](https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149/5 "2021-02-24T15:41:01Z")

</div>

Thanks, what is strange is that setting up a similar problem (adapted from here [Re: [petsc-users] Issues with the Variational Inequality solver](https://www.mail-archive.com/petsc-users@mcs.anl.gov/msg27359.html)) without going through FEniCS works:

```auto
import petsc4py
from petsc4py import PETSc
import numpy as np
psiA = 9.400282e-01
psiC = 6.045961e-09
k1 = 0.5
opt_VI = 1

#========================
# Variational inequality
#========================
class Reaction(object):
  def __init__ (self):
    self.psiA = 0
    self.psiC = 0
    self.k1 = 0
  def formData(self, x, a, c, k):
    X = x.getArray()
    X[0] = 0
    X[1] = 0
    X[2] = 0
    self.psiA = a
    self.psiC = c
    self.k1 = k
  def formFunction(self, snes, x, f):
    X = x.getArray(readonly=True)
    F = f.getArray()
    F[0] = X[0] - 2*X[1] - self.psiA
    F[1] = X[2] + 2*X[1] - self.psiC
    F[2] = pow(X[2],2) - self.k1*pow(X[0],2)*X[1]
  def formJacobian(self, snes, x, J, P):
    X = x.getArray(readonly=True)
    J.zeroEntries()
    J.setValues([0,1,2],[0,1,2], [1,-2,0,0,2,1,-2*self.k1*X[0]*X[1],-self.k1*pow(X[0],2),2*X[2]])
    J.assemble()
    return PETSc.Mat.Structure.SAME_NONZERO_PATTERN

pde = Reaction()
snes = PETSc.SNES().create(comm=PETSc.COMM_SELF)
A = PETSc.Mat().createDense(size=[3,3],comm=PETSc.COMM_SELF)
A.setUp()
J = A
F,X = A.createVecs()
snes.setFunction(pde.formFunction,F)
snes.setJacobian(pde.formJacobian,J,J)

vilb = X.duplicate()
viub = X.duplicate()
vilb.set(0)
viub.set(10)
snes.setType('vinewtonrsls')
snes.setVariableBounds(vilb,viub)

pde.formData(X,psiA,psiC,k1)
snes.solve(None, X)
print(X.array)

```

---

<div class="post-metadata">

**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\
**Post date:** [February 24, 2021, 3:43pm UTC](https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149/6 "2021-02-24T15:43:03Z")

</div>

Thanks Jorgen, that might be a good motivation for me to finally test dolfin-x

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 24, 2021, 3:55pm UTC](https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149/7 "2021-02-24T15:55:11Z")

</div>

No problem.  
If you want a summary of the syntax changes, have a look at my tutorial (under development): [The FEniCS-X tutorial — FEniCS-X tutorial](https://jorgensd.github.io/dolfinx-tutorial/)

---

<div class="post-metadata">

**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\
**Post date:** [February 25, 2021, 1:00pm UTC](https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149/8 "2021-02-25T13:00:41Z")

</div>

FYI, after some fiddling around, here is a version that does not segfaults. It seems that the origin is due to `b=PETScVector()` which has to be replaced by `b=Function(V).vector()`, although I don’t know exactly why. I also forgot some update on the d variable to obtain the good solution.

```python
from dolfin import *
from petsc4py import PETSc

N = 5
mesh = UnitSquareMesh(N, N)

V = FunctionSpace(mesh, "CG", 1)

d = Function(V)
dlb = Function(V)
dub = Function(V)
dub.interpolate(Constant(1.))

d_ = TestFunction(V)
dd = TrialFunction(V)

energy = ((d-0.5)**2)*dx
F = derivative(energy, d, d_)
J = derivative(F, d, dd)

b = Function(V).vector()
A = PETScMatrix()

def function(solver, x):
    return assemble(energy)
def gradient(solver, x, b):
    d.vector()[:] = x
    b = PETScVector(b)
    return assemble(F, tensor=b)
def hessian(solver, x, A, P):
    A = PETScMatrix(A)
    return assemble(J, tensor=A)

snes_solver = PETSc.SNES().create(PETSc.COMM_WORLD)
PETScOptions.set("snes_monitor")
snes_solver.setFromOptions()
snes_solver.setType(PETSc.SNES.Type.VINEWTONSSLS)
snes_solver.setFunction(gradient, b.vec())
snes_solver.setJacobian(hessian, A.mat())
snes_solver.setVariableBounds(dlb.vector().vec(), dub.vector().vec())
snes_solver.solve(None, d.vector().vec())
print("Succeeded SNES Solve", d.vector()[:])

```

---

<div class="post-metadata">

**Author:** ![manish\_07](https://avatars.discourse-cdn.com/v4/letter/m/e56c9b/32.png) [@manish\_07](https://fenicsproject.discourse.group/u/manish_07)\
**Post date:** [July 5, 2023, 8:26pm UTC](https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149/9 "2023-07-05T20:26:09Z")

</div>

I am trying to run the provide code in dolfinx 0.4 version. The `create_vector` and `create_matrix` are not the part of module and geting the following error

`module 'dolfinx.fem' has no attribute 'create_vector'`

so, I make some changes. The updated and working code for dolfinx 0.4 version as follows:

```auto
import dolfinx
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import mesh
N = 5
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

V = dolfinx.fem.FunctionSpace(domain, ('CG', 1))

def one(x):
    return np.ones(x.shape[1])

d = dolfinx.fem.Function(V)
dlb = dolfinx.fem.Function(V)
dub = dolfinx.fem.Function(V)
dub.interpolate(one)

d_ = ufl.TestFunction(V)
dd = ufl.TrialFunction(V)

energy = ((d-0.5)**2)*ufl.dx
F = ufl.derivative(energy, d, d_)
J = ufl.derivative(F, d, dd)

b = dolfinx.fem.petsc.create_vector(dolfinx.fem.form(F))
A = dolfinx.fem.petsc.create_matrix(dolfinx.fem.form(J))

def function(solver, x):
    return dolfinx.fem.petsc.assemble_scalar(energy)
def gradient(solver, x, b):
    b = b.copy()
    with b.localForm() as loc:
        loc.set(0)
    return dolfinx.fem.petsc.assemble_vector(b,dolfinx.fem.form(F))
def hessian(solver, x, A, P):
    A = A.copy()
    A.setZeroEntries()
    return dolfinx.fem.petsc.assemble_matrix(A, dolfinx.fem.form(J))

snes_solver = PETSc.SNES().create(MPI.COMM_WORLD)
opts = PETSc.Options()
opts["snes_monitor"] = None
snes_solver.setFromOptions()
snes_solver.setType(PETSc.SNES.Type.VINEWTONSSLS)
snes_solver.setFunction(gradient, b)
snes_solver.setJacobian(hessian, A)
snes_solver.setVariableBounds(dlb.vector, dub.vector)
snes_solver.solve(None, d.vector)
print("Succeeded SNES Solve")

```

Thank you,  
Manish