Setting SNES solver with petsc4py segfaults

Linear Algebra
1

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 = 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")
2

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

3

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

4

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

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

1 Like
5

Thanks, what is strange is that setting up a similar problem (adapted from here Re: [petsc-users] Issues with the Variational Inequality solver) without going through FEniCS works:

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)
6

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

7

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

8

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.

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

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:

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