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

239

views

0

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`.

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:

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:

Thanks for any insight.
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.

```
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)
scatter(rho0vals, prho0, addv=True, mode=PETSc.ScatterMode.SCATTER_FORWARD)
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
license : LGPL 3.0
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

written
10 months ago by
Leon Avery

### 1 Answer

2

OK, I have an answer, more or less. The part of the incantation I was missing is `apply`. Specifically, the following modification:

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.

```
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
10 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
10 months ago by
Leon Avery

1

Yes, that should do it.

written
10 months ago by
Garth Wells

Please login to add an answer/comment or follow this question.