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.

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]

[-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.