Imposing boundary condition manually with PETScMatrix fails in parallel for -np > 2


62
views
1
5 weeks ago by
Hi, I would like to solve a problem that involves boundary conditions of the kind \(u\cdot n = g\), where \(u\) is some vectorial quantity (velocity, for instance), \(n\) is the normal vector to the boundary, and \(g\) is a prescribed function. I want to impose this boundary condition by using a rotation of the Degrees of Freedom in the boundary nodes.

My approach is to assemble the weak form into a PETSc Matrix and to multiply it afterwards by a constructed Rotation matrix. Before going into that, I have decided  to work out a much simpler case before constructing the actual Rotation matrix. I assembled the weak form of the Poisson problem into a PETScMatrix, and then multiplied it by a constructed identity matrix. The imposition of the Dirichlet boundary conditions works fine for a number of processes smaller than 2. When I run the code using mpirun -np 3 python3 testPETScPar.py I get the following error:
****************************************************************************************************************************************
Process 1: PetscDolfinErrorHandler: line '649', function 'ISLocalToGlobalMappingApply', file '/tmp/petsc-src/src/vec/is/utils/isltog.c',
: error code '63' (Argument out of range), message follows:
Process 1: ------------------------------------------------------------------------------
Process 1: Local index 18 too large 17 (max) at 3
Process 1: ------------------------------------------------------------------------------
Process 1: PetscDolfinErrorHandler: line '603', function 'ISLocalToGlobalMappingApplyIS', file '/tmp/petsc-src/src/vec/is/utils/isltog.c',
: error code '63' (Argument out of range), message follows:
Process 1: ------------------------------------------------------------------------------
Process 1:
Process 1: ------------------------------------------------------------------------------
Process 1: PetscDolfinErrorHandler: line '6116', function 'MatZeroRowsLocal', file '/tmp/petsc-src/src/mat/interface/matrix.c',
: error code '63' (Argument out of range), message follows:
Process 1: ------------------------------------------------------------------------------
Process 1:
Process 1: ------------------------------------------------------------------------------
Process 1: PETSc error in '/tmp/src/dolfin/dolfin/la/PETScMatrix.cpp', 'MatZeroRowsLocal'
Process 1: PETSc error code '63' (Argument out of range), message follows:
Process 1: ------------------------------------------------------------------------------
Process 1: Local index 18 too large 17 (max) at 3
Process 1: ------------------------------------------------------------------------------
Process 1: Elapsed wall, usr, sys time: 0.000373529, 0, 0 (DirichletBC apply)


*** Error: Unable to successfully call PETSc function 'MatZeroRowsLocal'.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /tmp/src/dolfin/dolfin/la/PETScMatrix.cpp.
*** Process: 1
****************************************************************************************************************************************


I am using Dolfin version 2018.1.0.dev0. I got inspiration from https://fenicsproject.org/qa/11633/how-to-apply-dirichlet-bcs-to-matrix-generated-with-petsc4py/ but fails for np > 2. I would like to thank you beforehand for any suggestion you could make!
from petsc4py import PETSc
from dolfin import *
import sys
import time
from mpi4py import MPI
import numpy as np

set_log_level(LogLevel.TRACE)
mesh = UnitSquareMesh(6, 6)
V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

def boundary1(x):
    return x[0] < DOLFIN_EPS
def boundary2(x):
    return x[0] > 1.0 - DOLFIN_EPS

bc1 = DirichletBC(V, Constant(0.0), boundary1)
bc2 = DirichletBC(V, Constant(1.0), boundary2)

# Define variational problem
a = inner(grad(u), grad(v))*dx
L = Constant(0.0)*v*ds
A = assemble(a)
b = assemble(L)

A_mat = as_backend_type(A).mat()
Istart, Iend = A_mat.getOwnershipRange()
print("checkAmat", rank, Istart, Iend)

b_vec = as_backend_type(b).vec()

# Create Identity Matrix
sizes = [V.dofmap().index_map().size(IndexMap.MapSize.OWNED),
         V.dofmap().index_map().size(IndexMap.MapSize.GLOBAL)]

print(rank, "dofmap1", V.dofmap().dofs())
print(rank, "sizes", sizes)

sample_mat = PETSc.Mat()
sample_mat.create(comm)
sample_mat.setSizes([sizes, sizes])
sample_mat.setType('aij')
sample_mat.setUp()

#create local to global map
lgmap = PETSc.LGMap().create(V.dofmap().dofs(), comm=comm)
sample_mat.setLGMap(lgmap, lgmap)

#check ownership of nodes
Istart, Iend = sample_mat.getOwnershipRange()
print("checksample_mat", rank, Istart, Iend)

# Fill the values
vect = interpolate(Expression('1.', degree=1), V)
vect = as_backend_type(vect.vector()).vec()
vect.set(1.)
sample_mat.setDiagonal(vect)
sample_mat.assemblyBegin()
sample_mat.assemblyEnd()

#construct modified matrix
K = sample_mat*A_mat
K.setLGMap(lgmap,lgmap)
K = PETScMatrix(K)

# Compute solution
u = Function(V)

b_vec = PETScVector(b_vec)
bc1.apply(K, b_vec)
bc2.apply(K, b_vec)
solve(K, u.vector(), b_vec, "mumps")
Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »