Save the result of mpirun in a numpy array

I/O
1

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

2

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)
5 Likes
3

Dear dokken,

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

4

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

5

Hi dokken thank you for your quick support, I will try

6

hi dokken ,my problem is how can i store (coord,value) in csv file

7

As mentioned in the post above:

you can get it quite easily to print (and thus store it to a csv).

8

thanks ,but i already do this and i recieve a message like that:

“ZeroDivisionError: division by zero” i don’t know how can i do,i saw someone also was in the same problem like me,but yours responses i applied ,don’t work for me

9

please is there other method

10

Without a minimal example no one can help you.

11

ok that’s my problem
hello everyone, I want to store the nodal and mesh values in csv file.

ie (coor, u_nodal_values) in csv file.

but I have difficulties, I am a beginner, your suggestions will help me a lot. THANK YOU

mesh =RectangleMesh(Point(0, 0), Point(145, 45), nx, ny,‘left’)

V = FunctionSpace(mesh, ‘P’, 1)

u_nodal_values = u1.vector().get_local()

coor = mesh.coordinates()
for i in range(len(coor)):
x=coor[i][0]
y=coor[i][1]
T=u_array[i]
np.savetxt(‘dot.txt’,(x,y,T), delimiter=’,’)

but only save one line

12

Why aren’t you following my advice from:
How to get the solution obtained by fenics as a vector? - #2 by dokken (for scalar quantites)
How to get the solution obtained by fenics as a vector? - #7 by dokken (for vector quantities)

which would give you the following code:

from dolfin import *

mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, "CG", 2)
dof_coordinates = V.tabulate_dof_coordinates()

expr = Expression("x[0]", degree=1)
uh = Function(V)
uh.interpolate(expr)

out = open("data.csv", "w")
for x, val in zip(dof_coordinates, uh.vector().get_local()):
    print(f"{x[0]}, {x[1]}, {val}", file=out)
out.close()

giving you a text file:

0.0, 1.0, 0.0
0.0, 0.5, 0.0
0.5, 1.0, 0.5
0.25, 0.75, 0.25
0.25, 1.0, 0.25
0.0, 0.75, 0.0
0.0, 0.0, 0.0
0.5, 0.5, 0.5
0.25, 0.25, 0.25
0.25, 0.5, 0.25
0.0, 0.25, 0.0
0.5, 0.75, 0.5
1.0, 1.0, 1.0
0.75, 0.75, 0.75
0.75, 1.0, 0.75
0.5, 0.0, 0.5
0.5, 0.25, 0.5
0.25, 0.0, 0.25
1.0, 0.5, 1.0
0.75, 0.25, 0.75
0.75, 0.5, 0.75
1.0, 0.75, 1.0
1.0, 0.0, 1.0
1.0, 0.25, 1.0
0.75, 0.0, 0.75
1 Like
13

thank you ,this work well and help me to find answer for my problem
thank you again

14

Hi

Could you please make this example in dolfinx ?

15

This depends on what kind of function space you are using, as in DOLFINx, Nedelec/RT etc is using its proper definition being integral moments, which means that you do not have dof coordinates.

If you are using Lagrange,

dof_coordinates = V.tabulate_dof_coordinates()
uh = Function(V)

out = open("data.csv", "w")
for x, val in zip(dof_coordinates, uh.x.array):
    print(f"{x[0]} {x[1]}, {val}", file=out)
out.close()

This would yield a slight overlap (as it includes duplicate points on process boundaries in parallel).
Slicing the dof coordinates and array using the local size in the index map would resolve this (But as Im not at my laptop, i cannot write that code right now).

1 Like