Problem with petsc ksp.solve in parallel

Linear Algebra
1

Dear all,

I am trying to solve a Laplace problem in a square domain with a hole inside and I want to impose a Dirichlet boundary condition on the border of the hole using a Lagrange multiplier.

To do that, I construct a KSP solver (with petsc4py), from a shell matrix for which I give a matrix-vec multiplication rule.

The code below works well in sequential but fails in parallel (with mpirun -np 2 python).

This the error that I get:

Traceback (most recent call last):
  File "test.py", line 140, in <module>
    l.solve(B,X)
  File "test.py", line 120, in solve
    self.ksp.solve(B, X)
  File "PETSc/KSP.pyx", line 393, in petsc4py.PETSc.KSP.solve
petsc4py.PETSc.Error: error code 73

I have found very few reference to this error 73 and none of them was usefull in my case.

Do you have any idea on how to fix this error ? Or may be it is a bad idea to use petsc4py within FEniCS and you want to suggest a better solution ?

Cheers,
Fabien

Here is the mwe:

from dolfin import *
import petsc4py.PETSc as petsc
import mshr

xc = .5
yc = .5
rc = .2
f = Constant(1.)
mu = Constant(1.)
density = 50

comm = MPI.comm_world
rank = comm.rank

# Generate mesh
domain = mshr.Rectangle(Point(0,0), Point(1,1)) - mshr.Circle(Point(xc,yc),rc)
mesh = mshr.generate_mesh(domain,density)
with XDMFFile(comm,"laplace_mesh.xdmf") as outfile:
    outfile.write(mesh)

# Boundaries
class Dirichlet(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ( near(x[0], 0., 1e-14) or near(x[0], 1., 1e-14) or near(x[1], 0., 1e-14) or near(x[1], 1., 1e-14) )

class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return near( (x[0]-xc)*(x[0]-xc) + (x[1]-yc)*(x[1]-yc), rc*rc, mesh.hmin()/10.)

facets = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Interface().mark(facets,1)

# Interface mesh
mesh_i = MeshView.create(facets, 1)

# Define Problem
class Laplace:
    def __init__(self, mesh, mesh_i):
        self.mesh = mesh
        self.mesh_i = mesh_i
        Vi = FunctionSpace(self.mesh_i, "CG", 1)
        V = FunctionSpace(self.mesh, "CG", 1)
        self.W = MixedFunctionSpace(V,Vi)
        self.dim = self.W.sub_space(1).dim()
        self.first_dof, self.last_dof = self.W.sub_space(1).dofmap().ownership_range()

    def setup_problem(self, coeff, force, boundary, boundary_condition):
        # Define coefficient
        self.mu = coeff

        # Define external volumique force
        self.f = force

        # Define boundary condition
        ff = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
        boundary().mark(ff, 1)
        Interface().mark(ff,2)
        self.bc = DirichletBC(self.W.sub_space(0), boundary_condition, ff, 1)
        self.ds = Measure("ds", domain=self.W.sub_space(0).mesh(), subdomain_data=ff)
        self.dsi = Measure("dx", domain=self.W.sub_space(1).mesh())

        # Define problem
        u = TrialFunction(self.W.sub_space(0))
        self.v = TestFunction(self.W.sub_space(0))
        self.sol = Function(self.W.sub_space(0))

        a = self.mu*inner(grad(u), grad(self.v))*dx
        self.A = PETScMatrix()    
        self.b = PETScVector()
        assemble(a, tensor=self.A)
        self.bc.apply(self.A)

        # Define Lagrange multiplier
        self.lmbda = Function(self.W.sub_space(1))

    def setup_rhs(self, B):
        L = self.f*self.v*dx
        assemble(L, tensor=self.b)
        self.bc.apply(self.b)
        solve(self.A,self.sol.vector(),self.b)
        proj = Function(self.W.sub_space(1))
        LagrangeInterpolator.interpolate(proj,self.sol)
        B[self.first_dof:self.last_dof] = proj.vector().vec()[self.first_dof:self.last_dof]

    def setup_ksp(self, solver='gmres', prec='none'):
        M = petsc.Mat().createPython([self.dim,self.dim])
        M.setPythonContext(self)
        M.setFromOptions()
        M.setUp()

        PETScOptions.set("ksp_monitor_true_residual")
        self.ksp = petsc.KSP().create()
        self.ksp.setOperators(M)
        self.ksp.setType(solver)
        pc = self.ksp.getPC()
        pc.setType(prec)
        self.ksp.setTolerances(rtol=1e-6)
        self.ksp.setFromOptions()

    def mult(self, mat, X, Y):
        self.lmbda.vector().vec()[self.first_dof:self.last_dof] = X[self.first_dof:self.last_dof]
        proj = Function(self.W.sub_space(0))
        LagrangeInterpolator.interpolate(proj, self.lmbda)
        L = proj*self.v*self.ds(2)
        assemble(L, tensor=self.b)
        self.bc.apply(self.b)
        solve(self.A,self.sol.vector(),self.b)
        proj = Function(self.W.sub_space(1))
        LagrangeInterpolator.interpolate(proj,self.sol)
        Y[self.first_dof:self.last_dof] = proj.vector().vec()[self.first_dof:self.last_dof]
    
    def solve(self, B, X):
        self.ksp.solve(B, X)
        self.lmbda.vector().vec()[self.first_dof:self.last_dof] = X[self.first_dof:self.last_dof]
        proj = Function(self.W.sub_space(0))
        LagrangeInterpolator.interpolate(proj, self.lmbda)
        L = self.f*self.v*dx - proj*self.v*self.ds(2)
        assemble(L, tensor=self.b)
        self.bc.apply(self.b)
        solve(self.A,self.sol.vector(),self.b)


# Define Problem
l = Laplace(mesh,mesh_i)
B = Function(l.W.sub_space(1)).vector().vec()
X = B.copy()

# Setup and solve problem
zeros = Constant(0.)
l.setup_problem(mu, f, Dirichlet, zeros)
l.setup_rhs(B)
l.setup_ksp()
l.solve(B,X)
with XDMFFile(comm,"laplace.xdmf") as infile:
    infile.write(l.sol)

PS: The code makes use of the MeshView function coming from the mixed-dimensional branch of FEniCS. It has be tested in the docker mixed-dimensional branch container

2

Hi Fabien,

PETSc error code 73 means : “object in argument is in wrong state, e.g. unassembled mat”.
I haven’t tested your mwe. But to me the issue might be when you do :

B[self.first_dof:self.last_dof] = proj.vector().vec()[self.first_dof:self.last_dof]

In serial, all the dofs are owned by the only processor, so no problem. But the dofs ownership is more tricky in parallel.

3

Hi Cécile,

Thank you for your reply ! You are right, the problem is with dofs ownership.
Actually, It seems that the error is in the setup_ksp function. When I create the matrix with petsc, the dofs are not distributed in the same manner than fenics does for the mesh.

Do you know how I can create a petsc matrix which has the same dofs repartition between procs ?

4

Hi again,

I managed to construct the petsc matrix by assembling a dummy matrix with fenics, getting its sizes and giving them to the petsc matrix. If you have a better idea, I am interested. The new code is given at the end of this post, for further references.

However, I still have an issue when the hole inside the domain belongs to only one proc (because of the mesh partitioner). It seems that the interpolation step done inside the mult function (see code below) has trouble to succed: I do not get any error message but the code keeps running and never stops.

Is this problem with the interpolation in parallel a known issue ?

Cheers,
Fabien

from dolfin import *
import petsc4py.PETSc as petsc

xc = .75
yc = .75
rc = .1
f = Constant(1.)
mu = Constant(1.)
density = 50
compute_mesh = False

comm = MPI.comm_world
rank = comm.rank

if compute_mesh:
    import mshr
    domain = mshr.Rectangle(Point(0,0), Point(1,1)) - mshr.Circle(Point(xc,yc),rc)
    mesh = mshr.generate_mesh(domain,density)
    with XDMFFile(comm,"laplace_mesh.xdmf") as outfile:
        outfile.write(mesh)
    exit()

# Generate mesh
if rank==0: print("Reading the mesh...")
mesh = Mesh()
with XDMFFile(comm,"laplace_mesh.xdmf") as infile:
    infile.read(mesh)

# Boundaries
class Dirichlet(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ( near(x[0], 0., 1e-14) or near(x[0], 1., 1e-14) or near(x[1], 0., 1e-14) or near(x[1], 1., 1e-14) )

class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return near( (x[0]-xc)*(x[0]-xc) + (x[1]-yc)*(x[1]-yc), rc*rc, mesh.hmin()/30.)

partition = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
for c in cells(mesh):
    partition[c] = rank
with XDMFFile(comm,"partition.xdmf") as outfile:
    outfile.write(partition)

facets = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Interface().mark(facets,1)

# Interface mesh
mesh_i = MeshView.create(facets, 1)

# Define Problem
class Laplace:
    def __init__(self, mesh, mesh_i):
        self.mesh = mesh
        self.mesh_i = mesh_i
        Vi = FunctionSpace(self.mesh_i, "CG", 1)
        V = FunctionSpace(self.mesh, "CG", 1)
        self.W = MixedFunctionSpace(V,Vi)

    def setup_problem(self, coeff, force, boundary, boundary_condition):
        # Define coefficient
        self.mu = coeff

        # Define external volumique force
        self.f = force

        # Define boundary condition
        ff = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
        boundary().mark(ff, 1)
        Interface().mark(ff,2)
        self.bc = DirichletBC(self.W.sub_space(0), boundary_condition, ff, 1)
        self.ds = Measure("ds", domain=self.W.sub_space(0).mesh(), subdomain_data=ff)
        self.dsi = Measure("dx", domain=self.W.sub_space(1).mesh())

        # Define problem
        u = TrialFunction(self.W.sub_space(0))
        self.v = TestFunction(self.W.sub_space(0))
        self.sol = Function(self.W.sub_space(0))

        a = self.mu*inner(grad(u), grad(self.v))*dx
        self.A = PETScMatrix()    
        self.b = PETScVector()
        assemble(a, tensor=self.A)
        self.bc.apply(self.A)

        # Define Lagrange multiplier
        self.lmbda = Function(self.W.sub_space(1))

    def setup_rhs(self, B):
        L = self.f*self.v*dx
        assemble(L, tensor=self.b)
        self.bc.apply(self.b)
        solve(self.A,self.sol.vector(),self.b)
        proj = interpolate(self.sol, self.W.sub_space(1))   
        cproj = proj.vector().vec().getArray(readonly=1)
        print("proj", rank,proj.vector().vec().getOwnershipRange())
        BB = B.getArray(readonly=0)
        BB[:] = cproj[:]

    def setup_ksp(self, solver='gmres', prec='none'):
        fM = assemble(TrialFunction(self.W.sub_space(1))*TestFunction(self.W.sub_space(1))*dx)
        M = petsc.Mat()
        M.create(comm=petsc.COMM_WORLD)
        M.setSizes(as_backend_type(fM).mat().getSizes())
        M.setType(petsc.Mat.Type.PYTHON)
        M.setPythonContext(self)
        M.setUp()

        PETScOptions.set("ksp_monitor_true_residual")
        self.ksp = petsc.KSP().create()
        self.ksp.setType(solver)
        self.ksp.setOperators(M)
        pc = self.ksp.getPC()
        pc.setType(prec)
        self.ksp.setTolerances(rtol=1e-6)
        self.ksp.setFromOptions()

        x, b = M.getVecs()
        return x, b

    def mult(self, mat, X, Y):
        clmbda = self.lmbda.vector().vec().getArray(readonly=0)
        XX = X.getArray(readonly=1)
        clmbda[:] = XX[:]
        proj = Function(self.W.sub_space(0))
        LagrangeInterpolator.interpolate(proj, self.lmbda) # <---- Here is the issue when the hole belongs to only one proc
        L = proj*self.v*self.ds(2)
        assemble(L, tensor=self.b)
        self.bc.apply(self.b)
        solve(self.A,self.sol.vector(),self.b)
        proj = interpolate(self.sol,self.W.sub_space(1))
        cproj = proj.vector().vec().getArray(readonly=1)
        YY = Y.getArray(readonly=0)
        YY[:] = cproj[:]
    
    def solve(self, B, X):
        self.ksp.solve(B, X)
        clmbda = self.lmbda.vector().vec().getArray(readonly=0)
        XX = X.getArray(readonly=1)
        clmbda[:] = XX[:]
        proj = Function(self.W.sub_space(0))
        LagrangeInterpolator.interpolate(proj, self.lmbda)
        L = self.f*self.v*dx - proj*self.v*self.ds(2)
        assemble(L, tensor=self.b)
        self.bc.apply(self.b)
        solve(self.A,self.sol.vector(),self.b)


# Define Problem
l = Laplace(mesh,mesh_i)

# Setup and solve problem
zeros = Constant(0.)
if rank==0: print("setting up the problem...")
l.setup_problem(mu, f, Dirichlet, zeros)
if rank==0: print("setting up the KSP...")
x, b = l.setup_ksp()
if rank==0: print("setting up the right hand side of the linear problem...")
l.setup_rhs(b)
l.solve(b,x)
with XDMFFile(comm,"laplace.xdmf") as infile:
    infile.write(l.sol)