Hello,
I am working on a project where I need to save the Numpy arrays of degrees of freedom of FE functions and use them in another script. I am wondering how to get consistent results when using MPI.
Here's a minimal example of what I would like to do:
----------------------------------------
from firedrake import *
import numpy as np
mesh = UnitSquareMesh(3, 3)
V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V)
x, y = SpatialCoordinate(mesh)
f.interpolate((1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2))
a = (dot(grad(v), grad(u)) + v * u) * dx
L = f * v * dx
u = Function(V)
solve(a == L, u, solver_parameters={'ksp_type': 'cg'})
U = u.vector().gather()
if u.comm.rank == 0:
print(U)
name = 'U_global_' + str(u.comm.size) + '.npy'
np.save(name,U)
---------------------------------------------------------------------
If I run this in serial, and then with two processes in MPI, I get two files 'U_global_1.npy' and 'U_global_2.npy' but they do not agree. The first (in serial) has values:
[-0.29551776 -0.67467152 -0.67467152 -0.15155862 0.06606523 0.06606523
0.36508954 0.36508954 1.33298501 1.33298501 0.06606523 -0.15155862
0.06606522 -0.67467153 -0.67467153 -0.29551776]
while the second has values:
[-0.67467152 -0.67467152 -0.29551776 0.06606523 -0.15155862 0.06606523
0.06606523 -0.67467152 -0.15155862 -0.29551776 -0.67467152 0.06606523
1.33298501 0.36508954 0.36508954 1.33298502]
How do I get the parallel implementation to match the ordering of the d.o.f.s in the serial case? And avoiding this .gather() approach would be good, because I really only need the array gathered on a single process.
Thanks,
Peter