Generating DOLFIN matrices with PETSc and applying boundary conditions in parallel

Question

I am trying to apply Dirichlet boundary conditions to a matrix generated via petsc4py. E.g., this code,

import dolfin as dlf
from petsc4py import PETSc

mesh = dlf.UnitSquareMesh(100, 100)
W = dlf.FunctionSpace(mesh, 'CG', 1)

def boundary(x, on_boundary):
    return on_boundary

bc = dlf.DirichletBC(W, dlf.Constant(0.0), boundary)

N = mesh.num_vertices()
v = PETSc.Vec()
v.create()
v.setSizes(N) # PETSc error with np >= 3
v.setType('standard')
v.setValues(range(N), [N]*N)

B_pet = PETSc.Mat()
B_pet.createAIJ([N,N], nnz=N)
B_pet.setDiagonal(v)

lgmap = PETSc.LGMap().create(W.dofmap().dofs())
B_pet.setLGMap(lgmap, lgmap)

B_pet.assemblyBegin()
B_pet.assemblyEnd()

B = dlf.PETScMatrix(B_pet)
bc.apply(B) # DOLFIN error with np = 2

works when run in one process, but generates the following error when run on 2 processes:

*** Error:   Unable to apply boundary condition.
*** Reason:  Dimension of function space (10201) too large for application of boundary conditions to linear system (5151 rows).
*** Where:   This error was encountered inside BoundaryCondition.cpp.

and the following on 3 or more processes:

petsc4py.PETSc.Errorpetsc4py.PETSc.Error  File "PETSc/Vec.pyx", line 158, in petsc4py.PETSc.Vec.setSizes (src/petsc4py.PETSc.c:88897)
: error code 62
[0] VecSetSizes() line 1381 in /home/miguel/SourceCode/petsc/src/vec/vec/interface/vector.c
[0] Invalid argument
[0] Int value must be same on all processes, argument # 3

I know this has to do with the way the matrices/vectors are split between the different processes, but am not sure how to fix it.

Answer 1

Hi, consider adopting the following example

from petsc4py import PETSc
from dolfin import *

mesh = UnitSquareMesh(100, 100)
V = FunctionSpace(mesh, 'CG', 2)
bc = DirichletBC(V, Constant(0.0), 'on_boundary')
dofmap = V.dofmap()
# Diagonal vector
v = interpolate(Expression('pow(x[0], 2)-x[1]', degree=2), V)
v = as_backend_type(v.vector()).vec()

comm = mesh.mpi_comm().tompi4py() 
mat = PETSc.Mat()
mat.create(comm)
# Local, global
sizes = [dofmap.index_map().size(IndexMap.MapSize_OWNED), 
         dofmap.index_map().size(IndexMap.MapSize_GLOBAL)]
# Square
mat.setSizes([sizes, sizes])
# Sparse
mat.setType('aij')
mat.setUp()
# Map from local rows to gloval rows
lgmap = map(int, dofmap.tabulate_local_to_global_dofs())
lgmap = PETSc.LGMap().create(lgmap, comm=comm)
mat.setLGMap(lgmap, lgmap)

# Fill the values
mat.setDiagonal(v)
mat.assemblyBegin()
mat.assemblyEnd()

A = PETScMatrix(mat)
bc.apply(A)

# Check if we can do matvec. Just no crash
x, y = mat.createVecs()
x.setRandom(); y.setRandom()
x, y = PETScVector(x), PETScVector(y)
A.mult(x, y)

# Verify the matrix
# The diagonal entry should be coordinates of the corresponing dof.
X = V.tabulate_dof_coordinates().reshape((-1, 2))
x, y = X[:, 0], X[:, 1]
values = x**2 - y  # Just like the interpolated expression
# For a row corresponding to a dof where bc was used the value should be 1
bc_dofs = bc.get_boundary_values().keys()
first, last = dofmap.ownership_range()
for ldof in range(last-first):
    gdof = dofmap.local_to_global_index(ldof)
    if first <= gdof < last:
        indices, row_values = A.getrow(gdof)
        assert len(indices) == len(row_values) == 1
        assert indices[0] == gdof

        if ldof in bc_dofs:
            assert int(row_values[0]) == 1
        else:
            assert near(values[ldof]-row_values[0], 0)