This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Import HDF5 mesh without partition it (MPI). Easy with an xml file, but with a HDF5?

+1 vote

Hello everyone,

Please consider my issue:

In a first file, i save a mesh with the format HDF5 (the code uses MPI)

hdf = HDF5File(mesh_phase_A.mpi_comm(), 'mesh_A.h5', "w")
    hdf.write(mesh_phase_A, "mesh")

In a second file (with MPI), i load the mesh

mesh_partitioned = dolfin.Mesh()
hdf5 = dolfin.HDF5File(dolfin.mpi_comm_world(), 'mesh_A.h5, 'r')
hdf5.read(mesh_partitioned, '/mesh', False)

The mesh is paritioned between all the process with this syntax.
I would like also to have the mesh non-partitioned.

If the mesh to be loaded was in a xml format, i could have done simply:

if rank_process==0:
    mesh_nonpartitioned = Mesh(mpi_comm_self(),'mesh_A.xml.gz') 

Then: what is the good syntax to import (in a MPI program) a mesh in a HDF5 format, without partition it ?

Regards.

asked Apr 28, 2017 by fussegli FEniCS Novice (700 points)

1 Answer

+1 vote

What about

meshname = meshname
mesh = Mesh(meshname + ".xml")
hdf = HDF5File(mesh.mpi_comm(), meshname+".h5", "w")
hdf.write(mesh, "/mesh")
hdf.close()

and then

meshname = meshname
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(),  meshname+".h5", "r")
hdf.read(mesh, "/mesh", False)

?

answered May 1, 2017 by Claire L FEniCS User (2,120 points)

This is the solution. A little addition: If you don't want to load the same mesh on each process, you can do something like this:

import dolfin as df
mesh_name = 'fullspace'

mpi_comm = df.mpi_comm_world()
my_rank = df.MPI.rank(mpi_comm)

if my_rank==0:
    mesh = df.Mesh()
    hdf = df.HDF5File(df.mpi_comm_self(),'../meshes/_h5/' + mesh + '.h5', 'r')
    hdf.read(mesh, "/mesh", False)
else:
    mesh= None

df.MPI.barrier(df.mpi_comm_world())
print(my_rank, mesh_np)

The output for mpirun -n 4 would be something like:

(0, dolfin.cpp.mesh.Mesh *' at 0x7fab0f478150> >)
(1, None)
(2, None)
(3, None)

instead of the following one without the if / else statements.

(2, dolfin.cpp.mesh.Mesh *' at 0x7f4982f68150> >)
(3, dolfin.cpp.mesh.Mesh *' at 0x7fcb0726f150> >)
(0, dolfin.cpp.mesh.Mesh *' at 0x7f125c34a150> >)
(1, dolfin.cpp.mesh.Mesh *' at 0x7f4c612ac150> >)

Indeed i can read it.
However, if i try to assemble or plot the mesh, then the code stops (there is no error message, it just happen nothing, through CPU is being used, and i have to manually stop the program).
So, with this method, yes i can 'read' it, but i can't manipulate it.

Are you using something like the following syntax?:

if my_rank==0:
    do something with the mesh
else
    pass

Otherwise it won't work, since the other processes "don't know" the mesh.

Alternatively, you can think about broadcasting (mpi broadcast) the mesh to all the processes or spawn a single process doing some work with help of mpi spawn (from mpi4py package).

...