# Petsc Snes solver in dolfinx

**URL:** <https://fenicsproject.discourse.group/t/petsc-snes-solver-in-dolfinx/11677>\
**Category:** dolfinx\
**Created:** [July 6, 2023, 11:21am UTC](https://fenicsproject.discourse.group/t/petsc-snes-solver-in-dolfinx/11677 "2023-07-06T11:21:41Z")\
**Posts on this page:** 2\
**Page:** 1

<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 6, 2023, 11:21am UTC](https://fenicsproject.discourse.group/t/petsc-snes-solver-in-dolfinx/11677/1 "2023-07-06T11:21:41Z")

</div>

Hello everyone,

I am trying to use the petsc snes solver in the dolfinx. I am very new to the dolfinx. I am taking the help from the old discussion

> [@Setting SNES solver with petsc4py segfaults](https://fenicsproject.discourse.group/t/setting-snes-solver-with-petsc4py-segfaults/5149):
>
> 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 ? 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 = derivativ…

The code for petsc snes in dolfin is working and provides the desired results. Here is the code I am using for that

```auto
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()[:])

```

In similar way, the code for dolfinx is developed and the code is 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):
# x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        #print('q')
 # x.copy(d.vector)
  # d.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
    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")

```

The code in dolfinx is running without any error but the snes solver is giving zero norm in the first iteration which is not correct. Please, guide me where is making mistake.

Thank you,  
Manish

---

<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:** [July 6, 2023, 2:02pm UTC](https://fenicsproject.discourse.group/t/petsc-snes-solver-in-dolfinx/11677/2 "2023-07-06T14:02:41Z")

</div>

See for example [dolfinx/python/test/unit/nls/test\_newton.py at main · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/main/python/test/unit/nls/test_newton.py#L59-L86) and corresponding [dolfinx/python/test/unit/nls/test\_newton.py at main · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/main/python/test/unit/nls/test_newton.py#L179-L202).