Petsc Snes solver in dolfinx

dolfinx
1

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

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

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:

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

2

See for example dolfinx/python/test/unit/nls/test_newton.py at main · FEniCS/dolfinx · GitHub and corresponding dolfinx/python/test/unit/nls/test_newton.py at main · FEniCS/dolfinx · GitHub.

1 Like