Save the result of mpirun in a numpy array

Hi guys,

I am written a piece of code which simply solves the Laplace equation in a 3D cube and obtains the current density distribution inside it. The code runs in parallel on HPC and works perfectly. I want to save the result (the current density vector) and then call it in another code which is not in FEniCS and runs in serial. In fact, I need nodal values for each mesh coordinate in a numpy array.

First, I tried to extract the nodal values and save them by np.save. However, in this way, only a part of data was saved (based on the MPI rank).

Second, I tried to use h5py to save data. However, it was not compatible with HDF5 embedded in FEniCS 2019 and the code couldn’t be executed.

Third, I saved data in HDF5 file. Then, I tried to read it in the serial code. However, in this way, I had to define the mesh and function spaces in the serial code. This way was time-consuming and reduced code performance drastically.

Now, I am looking for a way to save the data obtained from running my code in parallel as a numpy array. Any possible help would be highly appreciated.

Best wishes,
-Hassan

Here is a simple example that gathers the data on the first process:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, "CG", 1)
u = project(Expression("x[0]", degree=1), V)
vertex_values = u.compute_vertex_values()
coordinates = mesh.coordinates()

from mpi4py import MPI
gathered_coordinates = MPI.COMM_WORLD.gather(coordinates, root=0)
gathered_values = MPI.COMM_WORLD.gather(vertex_values, root=0)
global_indices = MPI.COMM_WORLD.gather(mesh.topology().global_indices(0))
if MPI.COMM_WORLD.rank == 0:
    num_vertices = mesh.num_entities_global(0)
    all_coordinates = np.zeros((num_vertices, coordinates.shape[1]))
    all_values = np.zeros(num_vertices)
    for coord, value, indices in zip(gathered_coordinates, gathered_values, global_indices):
        all_values[indices] = value
        all_coordinates[indices] = coord

    for (coord, value) in zip(all_coordinates, all_values):
        print(coord, value)
1 Like

Dear dokken,

Thank you very much for the answer. I will check it and come back to you if there is any problem.

Dear dokken,

Your solution worked perfectly.
Since I work with a vector function space (Current density vector field), I applied a minimal change to your solution as follows:

from dolfin import *
import numpy as np
mesh = UnitCubeMesh(MPI.comm_world,4,4,4)
V = VectorFunctionSpace(mesh, "DG", 0)
u = Function(V)
u.interpolate(Expression(('3*x[0]', '2*x[1]', 'x[2]*x[0]'), degree=2))

vertex_values0 = u.sub(0, deepcopy=True).compute_vertex_values()
vertex_values1 = u.sub(1, deepcopy=True).compute_vertex_values()
vertex_values2 = u.sub(2, deepcopy=True).compute_vertex_values()
coordinates = mesh.coordinates()


gathered_coordinates = MPI.comm_world.gather(coordinates, root=0)
gathered_values0 = MPI.comm_world.gather(vertex_values0, root=0)
gathered_values1 = MPI.comm_world.gather(vertex_values1, root=0)
gathered_values2 = MPI.comm_world.gather(vertex_values2, root=0)
global_indices = MPI.comm_world.gather(mesh.topology().global_indices(0))
if MPI.comm_world.rank == 0:
    num_vertices = mesh.num_entities_global(0)
    all_coordinates = np.zeros((num_vertices, coordinates.shape[1]))
    all_values = np.zeros((num_vertices,coordinates.shape[1]))
    for coord, value0, value1, value2, indices in zip(gathered_coordinates, gathered_values0, gathered_values1, gathered_values2, global_indices):
        all_values[indices] = np.transpose([value0, value1, value2])  # dedicate each vertex value to corresponding global index
        all_coordinates[indices] = coord

    for (coord, value) in zip(all_coordinates, all_values):
        print(coord, value)

Best wishes,
-Hassan

P.S. I got help from your answer in this topic:
https://fenicsproject.discourse.group/t/how-to-get-the-solution-obatained-by-fenics-as-a-vector/5377/7