### project gives wrong results when run on 2 parallel processes

193
views
0
8 months ago by
I'm trying, for the first time, to get FEniCS working under MPI. In the following code I generate a Function that is 1 everywhere then project it onto a higher-degree FunctionSpace. I begin by generating the DOF vector in one local vector rho0vals, then scatter it to a distributed vector rho0. To confirm that this is working correctly, I gather rho0 into a local vector lrho0, and it is indeed what it should be: all 1s. Then I project the Function onto a degree 2 FunctionSpace.

import mpi4py.MPI as MPI
import numpy as np
from os import path
from datetime import datetime
import petsc4py, sys
#
# this needs to be done before importing PETSc
#
petsc4py.init(sys.argv)
from petsc4py import PETSc
import fenics as fe
from fenics import (UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh,
FiniteElement, MixedElement, VectorElement,
TestFunctions, FunctionAssigner, GenericFunction,
Expression, FunctionSpace, VectorFunctionSpace,
TrialFunction, TestFunction, Function, Constant,
Measure, FacetNormal, CellSize, PETScVector,
PETScMatrix)

def log(*args, **kwargs):
rank = MPI.COMM_WORLD.Get_rank()
print('rank = ', rank, *args, **kwargs, flush=True)

def main():
log('', '*'*60)
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
nelements = 8
degree = 2
murho0 = 1.0
mesh = UnitIntervalMesh(nelements)
fe.parameters['linear_algebra_backend'] = 'PETSc'
Cd1 = fe.FunctionSpace(mesh, 'CG', 1)
rho0 = Function(Cd1)
size = rho0.vector().size()
log('rho0.vector().size()', rho0.vector().size())
prho0 = fe.as_backend_type(rho0.vector()).vec()
rho0vals = PETSc.Vec().createSeq(size)
if (rank == 0):
vals = np.array([1.0]*size)
rho0vals[:] = vals
else:
rho0vals[:] = np.zeros(size)
log('rho0vals.array', rho0vals.array)
log('prho0.array (before scatter)', prho0.array)
indexes = PETSc.IS().createStride(size=size, step=1)
log('indexes.array', indexes.array)
scatter = PETSc.Scatter().create(rho0vals, indexes, prho0, indexes)
prho0.assemble()
log('prho0.array (after scatter)', prho0.array)
log('rho0.vector().array().shape', rho0.vector().array().shape)
log('rho0.vector().array()', rho0.vector().array())
lrho0 = fe.Vector(fe.mpi_comm_self())
rho0.vector().gather(lrho0, np.array(range(size), 'intc'))
log('lrho0.array()', lrho0.array())
CE = FiniteElement('CG', mesh.ufl_cell(), degree)
CS = FunctionSpace(mesh, CE)
f = fe.project(rho0, CS, solver_type='gmres')
allf = fe.Vector(fe.mpi_comm_self())
log('CS.dim())', CS.dim())
f.vector().gather(allf, np.array(range(CS.dim()), 'intc'))
log('f.vector().array()', f.vector().array())
log('allf.array()', allf.array())

if __name__ == "__main__":
# execute only if run as a script
main()
​

Run as a single process, this does about what I expect. The projected Function, f, ends up close to 1 everywhere. But when I run it on two processes, as follows:

mpiexec -n 2 python min_project_error.py

here's the result:

rank =  0  ************************************************************
Building mesh (dist 0a)
rank =  1  ************************************************************
Building mesh (dist 1a)
rank =  0 rho0.vector().size() 9
rank =  1 rho0.vector().size() 9
rank =  1 rho0vals.array [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
rank =  0 rho0vals.array [ 1.  1.  1.  1.  1.  1.  1.  1.  1.]
rank =  1 prho0.array (before scatter) [ 0.  0.  0.  0.]
rank =  0 prho0.array (before scatter) [ 0.  0.  0.  0.  0.]
rank =  0 indexes.array [0 1 2 3 4 5 6 7 8]
rank =  1 indexes.array [0 1 2 3 4 5 6 7 8]
rank =  1 prho0.array (after scatter) [ 1.  1.  1.  1.]
rank =  0 prho0.array (after scatter) [ 1.  1.  1.  1.  1.]
rank =  1 rho0.vector().array().shape (4,)
rank =  0 rho0.vector().array().shape (5,)
rank =  0 rho0.vector().array() [ 1.  1.  1.  1.  1.]
rank =  1 rho0.vector().array() [ 1.  1.  1.  1.]
rank =  0 lrho0.array() [ 1.  1.  1.  1.  1.  1.  1.  1.  1.]
rank =  1 lrho0.array() [ 1.  1.  1.  1.  1.  1.  1.  1.  1.]
rank =  0 CS.dim()) 17
rank =  1 CS.dim()) 17
rank =  1 f.vector().array() [ 0.4267762   1.08578788  1.01473136  0.98743545  1.00260119  0.99783115
1.00087236  0.99956307]
rank =  0 f.vector().array() [ 0.50000003  0.91420872  1.07322352  0.98526743  1.01257283  0.99740252
1.00216333  0.9991271   1.0004357 ]
rank =  1 allf.array() [ 0.50000003  0.91420872  1.07322352  0.98526743  1.01257283  0.99740252
1.00216333  0.9991271   1.0004357   0.4267762   1.08578788  1.01473136
0.98743545  1.00260119  0.99783115  1.00087236  0.99956307]
rank =  0 allf.array() [ 0.50000003  0.91420872  1.07322352  0.98526743  1.01257283  0.99740252
1.00216333  0.9991271   1.0004357   0.4267762   1.08578788  1.01473136
0.98743545  1.00260119  0.99783115  1.00087236  0.99956307]
​

These numbers are wrong. At the left end of each local vector, the number is about half what it ought to be. What am I doing wrong?

I'm using the anaconda fenics distribution on a Mac mini with MacOS Sierra.The backend is PETSc.  Here's what conda has to say about the version:

fenics 2017.1.0 py36_4
----------------------
file name   : fenics-2017.1.0-py36_4.tar.bz2
name        : fenics
version     : 2017.1.0
build string: py36_4
build number: 4
channel     : conda-forge
size        : 5.1 MB
arch        : x86_64
has_prefix  : True
md5         : af9cb33f8e488bf5fe826656708278dd
noarch      : None
platform    : darwin
requires    : ()
subdir      : osx-64
url         : https://conda.anaconda.org/conda-forge/osx-64/fenics-2017.1.0-py36_4.tar.bz2
​

Thanks for any insight.
Community: FEniCS Project
Thinking the problem might have something to do with ghost values, I tried it with fe.parameters['ghost_mode'] = 'shared_vertex' and fe.parameters['ghost_mode'] = 'shared_facet'. This had no effect.
written 8 months ago by Leon Avery

2
7 months ago by
OK, I have an answer, more or less. The part of the incantation I was missing is apply. Specifically, the following modification:

    rho0.vector().apply('')

done just before the project causes sensible answers to be produced, both running as a single process and as two (or even three).

I'm still a bit frustrated by the (as far as I can locate) total absence of documentation of the apply method. It apparently requires a single string argument, but what this string signifies I have no idea. In the hope of provoking an informative error message, I've tried '', 'xxx', and 'rong', with identical results.
apply finalise the vector (off-process vector entries are exchanged). You also need to take care that qhost entries in the vector are updated.
written 7 months ago by Garth Wells
Thank you. What do I need to do to ensure that ghost entries are updated? I noticed in PETScVector.cpp that apply calls update_ghost_values, which in turn calls PETsc functions to check if the vector is ghosted and update it. Is anything more necessary?
written 7 months ago by Leon Avery
1
Yes, that should do it.
written 7 months ago by Garth Wells