MPI crashes with vector() for mesh loaded from file

Question

Hi,

I am new to MPI, it looks seductively easy to just run a script with mpirun but I am running into (in my view) strange errors:

I am working in python on Ubuntu, with Fenics 1.5 from the fenics ppa. The program itself is quite long but I have managed to reproduce the error with a short snippet (below).

Basically, if I define a RectangularMesh in the code, things execute with MPI. But if I load an external mesh (which I really need to) the program crashes if I run with mpirun -n X with X > 1. If X=1 all it executes nicely.
I could add that the external mesh has significantly fewer cells than the rectangle in the snippet.

The offender seems to be u.vector(). I am lost!

Advice would be much appreciated!

The snippet:

from dolfin import *
from dolfin.cpp.mesh import RectangleMesh

mesh = Mesh('external.xml')
#mesh = RectangleMesh(0.0,0.0, 1.0, 1.0, 1000, 1000)
print mesh

V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)

u = Function(V)
t = u.vector()

print "OK til end"

The error thrown:

  File "test.py", line 21, in <module>
    t = u.vector()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to access vector of degrees of freedom.
*** Reason:  Cannot access a non-const vector from a subfunction.
*** Where:   This error was encountered inside Function.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.5.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

Looking at the code of Function.cpp, line 339, a check is done

if (_vector->size() != _function_space->dofmap()->global_dimension())

Perhaps I need to help the program along a bit to set up the dofmap for example?
Something's missing... — gobetween, Apr 27, 2015
I have figured out that this problem does not occur with all input meshes.

If I save a simple RectangularMesh generated by Fenics to XML and instead load that everything works fine. I.e. it is perhaps not the xml loading per se that is the culprit but could be a "bad" mesh (?)

Any ideas what could be causing this, or how I can analyse the mesh and/or code to find where things go wrong? — gobetween, May 1, 2015
More clues to add to the puzzle:

the same input mesh causes PETSc error 56 if I switch to KrylovSolver CG/ILU. — gobetween, May 2, 2015

MPI crashes with vector() for mesh loaded from file

1 Answer