MPI Error: 'Invalid rank, error stack' when reading mesh

Hi everyone,
I am trying to get my first 3D project to work. I have this simple script to just read a mesh I created from a .geo file in Gmsh 4.13.1 (.geo file at the very end):

from mpi4py import MPI
import dolfinx.io.gmshio
mesh, cell_tags, facet_tags = dolfinx.io.gmshio.read_from_msh("3D-chamber4.msh", MPI.COMM_WORLD, 0, gdim=3)
print(mesh, cell_tags, facet_tags)

However, I keep running into this error that I have read a couple of posts about but couldn’t figure it out:

Info : Done reading ‘3D-chamber4.msh’
Invalid rank, error stack:
internal_Issend(61644): MPI_Issend(buf=0x5632a62ba6e1, count=1, MPI_BYTE, 1, 1, comm=0x84000005, request=0x5632a62b9024) failed
internal_Issend(61605): Invalid rank has value 1 but must be nonnegative and less than 1
Abort(1010434822) on node 0 (rank 0 in comm 416): application called MPI_Abort(comm=0x84000005, 1010434822) - process 0

Trying to run the .py with mpirun -n 2 python3 read-mesh.py did not help. From other posts I read that it is usually a problem with unassigned physical surfaces. However, if I comment out all the physical surfaces in the .geo file, the read-mesh.py suddenly works and the output is

Info : Done reading ‘3D-chamber4.msh’
<dolfinx.mesh.Mesh object at 0x7fef52a31fd0> <dolfinx.mesh.MeshTags object at 0x7fef52a31e80> <dolfinx.mesh.MeshTags object at 0x7fef52bfafd0>

I’ve also tried to add Recombine; in the extrude command, but that just led to mixed element types and a different error message. The current .msh file only has Tetrahedas for 3D elements.

SetFactory("OpenCASCADE");

// Mesh sizing
h=1;
a=0.01;
Mesh.CharacteristicLengthMin = a/2*h;
Mesh.CharacteristicLengthMax = a*h;//+

Circle(1) = {0, 0, 0, 0.035, 0, 2*Pi};
Circle(2) = {0, 0, 0, 0.01, 0, 2*Pi};
Circle(3) = {0, 0, 0, 0.0127, 0, 2*Pi};

// Chamber inlet annulus
Curve Loop(1) = {1};
Curve Loop(2) = {3};
Plane Surface(1) = {1, 2};

// Fuel inlet
Curve Loop(3) = {3};
Curve Loop(4) = {2};
Plane Surface(2) = {3, 4};
Physical Surface(12) = {2};

// Lox inlet
Curve Loop(5) = {2};
Plane Surface(3) = {5};
Physical Surface(11) = {3};

//Extrude
Extrude {0, 0, 0.2} {
  Curve{1}; Layers{10};
}

// Outlet
Curve Loop(7) = {5};
Plane Surface(5) = {7};
Physical Surface(13) = {5};

// Walls
Physical Surface(14) = {4,1};

// Volume
Surface Loop(1) = {4, 5, 1, 2, 3};
Volume(1) = {1};
Physical Volume(15) = {1};

This is due to the fact that you are not tagging the correct surfaces in your geo file.
This is a common mistake when using the GMSH geo format, rather than their C++ or Python API.

This can be illustrated with:

import meshio
import numpy as np

mesh = meshio.read("mesh.msh")

tetrahedrons = mesh.cells_dict["tetra"]
triangles = mesh.cells_dict["triangle"]

tetra_vertices = np.unique(tetrahedrons.flatten())
tri_vertices = np.unique(triangles.flatten())
print(np.isin(tri_vertices, tetra_vertices))

for your mesh, you will observe that many of the vertices in the triangles are not part of any tetrahedra.

You should use coherence to merge duplicate nodes:

Hi there!

I’m getting a similar error while generating a simple 3D cylinder through the Python API.

gmsh.initialize()

L = 2.2
D = 0.41
gdim = 3
mesh_comm = MPI.COMM_WORLD
model_rank = 0

if mesh_comm.rank == model_rank:
    cylinder = gmsh.model.occ.addCylinder(
        0, 0, 0,
        L, 0, 0,
        D,
        tag=1,
    )
    gmsh.model.occ.synchronize()

fluid_marker = 1
if mesh_comm.rank == model_rank:
    volumes = gmsh.model.getEntities(dim=gdim)
    assert len(volumes) == 1
    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
    gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")

inlet_marker, outlet_marker, wall_marker = 2, 3, 4
inflow, outflow, walls = [], [], []
if mesh_comm.rank == model_rank:
    boundaries = gmsh.model.getBoundary(volumes, oriented=False)
    for boundary in boundaries:
        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
        if np.allclose(center_of_mass, [0, 0, 0]):
            inflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L, 0, 0]):
            outflow.append(boundary[1])
        else:
            walls.append(boundary[1])
    gmsh.model.addPhysicalGroup(2, walls, wall_marker)
    gmsh.model.setPhysicalName(2, wall_marker, "Walls")
    gmsh.model.addPhysicalGroup(2, inflow, inlet_marker)
    gmsh.model.setPhysicalName(2, inlet_marker, "Inlet")
    gmsh.model.addPhysicalGroup(2, outflow, outlet_marker)
    gmsh.model.setPhysicalName(2, outlet_marker, "Outlet")

if mesh_comm.rank == model_rank:
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.08)
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(2)
    gmsh.model.mesh.optimize("Netgen")
    gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
    gmsh.write("/output/test.msh")

mesh_data = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)

Executing this code leads to the following error:

Invalid rank, error stack:
internal_Issend(117): MPI_Issend(buf=0x1942ec91, count=1, MPI_BYTE, 1, 1202, comm=0x84000001, request=0x193dc524) failed
internal_Issend(78).: Invalid rank has value 1 but must be nonnegative and less than 1
Abort(943326982) on node 0 (rank 0 in comm 496): application called MPI_Abort(comm=0x84000001, 943326982) - process 0

Note that I adapted the code from Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark).

Any help would be greatly appreciated!

Thanks :slight_smile:

Hi guys!

It looks like my issue is coming from the use of gmsh.model.mesh.optimize("Netgen"), which erases the second order nature of the tetra elements:

<meshio mesh object>
  Number of points: 6452
  Number of cells:
    triangle6: 2612
    tetra: 10921
  Cell data: gmsh:physical, gmsh:geometrical
  Field data: Inlet, Outlet, Walls, Fluid

Deleting this line leads to the following output with no error:

<meshio mesh object>
  Number of points: 17966
  Number of cells:
    triangle6: 2612
    tetra10: 11547
  Cell data: gmsh:physical, gmsh:geometrical
  Field data: Inlet, Outlet, Walls, Fluid

Good to know that Netgen optimization should not be used when dealing with second order volumetric elements (at least for tetras).

You can always use the order increase after the optimization to make it curved again

That’s true :smiley:

But in the original example I used about the DFG 2D-3 benchmark, the Netgen optimization is done after the order increase, without any issue.

So I guess the issue arise only on 3D elements.

Might be useful to update the example by performing the Netgen optimization before the order increase, so other people like me won’t get confused when trying to adapt the example in 3D.

Anyway thanks for you feedback @dokken!

You people are doing an amazing job with the FEniCS project, and I’m looking forward to deepening more on this technology.