MPI/gathering function vector data for numpy operations

Question

Hello,
I would like to access a Function vector and perform numpy operations on values assembled from all MPI processes. Specifically, I'd like to be able to store data (interpolating the vector data onto a line) into a text/pickle file that I can load

import numpy as np 
newX = np.asarray(x.vector()[:])
print newX*1. 

Is there a recommended way to do this? I've attached a minimal(ish) example that fails for more than 1 processor? Thanks in advance
Pete

from dolfin import *
dolf=1

# Define Dirichlet boundary (x = 0 or x = 1)
class leftBoundary(SubDomain):
  def inside(self,x,on_boundary):
    outer  = x[0] < DOLFIN_EPS 
    return outer and on_boundary

# Define Dirichlet boundary (x = 0 or x = 1)
class elseBoundary(SubDomain):
  def inside(self,x,on_boundary):
    outer  = x[0] >=DOLFIN_EPS 
    return outer and on_boundary

class rightBoundary(SubDomain):
  def inside(self,x,on_boundary):
    outer  = x[0] >=(1-DOLFIN_EPS)
    return outer and on_boundary

# Create mesh and define function space
def doit(filename="none"):

  # Create mesh and define function space
  if(filename=="none"):
    #mesh = UnitCube(8,8,8)
    mesh = UnitCubeMesh(16,16,16)
  V = FunctionSpace(mesh, "Lagrange", 1)
  dims = mesh.ufl_cell().geometric_dimension()

  # Define boundary condition
  subdomains = MeshFunction("size_t",mesh,dims)
  subsurfaces = MeshFunction("size_t",mesh,dims-1)
  boundary = leftBoundary()
  leftmarker = 2
  boundary.mark(subsurfaces,leftmarker)
  boundary = rightBoundary()
  elsemarker = 3
  boundary.mark(subsurfaces,elsemarker)
  ds = Measure("ds")[subsurfaces] 
  dx = Measure('dx')[subdomains]

  # save 

  bcs = []
  u0 = Constant(1.0)
  bcs.append(DirichletBC(V, u0, subsurfaces,leftmarker))
  u0 = Constant(0.0)
  bcs.append(DirichletBC(V, u0, subsurfaces,elsemarker))


  # Define variational problem
  u = TrialFunction(V)
  v = TestFunction(V)

  # Solve 
  pv = Constant(1.)
  if(dolf==1):
    form = pv * inner(grad(u), grad(v))*dx(domain=mesh)
    g = Constant(0.)
    form += g*v*dx 

  # Compute solution
  a = lhs(form)
  L = rhs(form)
  x = Function(V)
  solve(a == L, x, bcs)

  tot = assemble(x*dx)
  vol = assemble(Constant(1.)*dx(domain=mesh))
  conc = tot/vol
  print "Tot/Vol/Conc %f %f %f" %(tot,vol,conc)

  # NOT MPI SAFE 
  #print x([0.5,0.5,0.5])
  if MPI.rank(mpi_comm_world())==0: 
    # still unsafe 
    #print x([0.5,0.5,0.5])
    1
    import numpy as np 
    newX = np.asarray(x.vector()[:])
    print newX*1.



  return x 

doit()

Answer 1

It may not fulfil your requirements, but one option is to save the Function using HDF5File which correctly sorts out the mapping in parallel. You can explore the file using things like h5py or just read back in in serial to postprocess.

Hdf=HDF5File(mesh.mpi_comm(), "a.h5", "w")
Hdf.write(u, "fun_u")

I think

Comments

Hi Chris - thanks for your response. I am currently running 1.4, which apparently doesn't have HDF5File. Is there a similar function in 1.4? I had trouble building 1.5, so I'd like to see what 1.4 could do before reimmersing myself in installing 1.5
Thanks!
pete

---

It should be in 1.4, but it us an optional dependency - you need to have a parallel hdf5 library installed too.

---

Thanks! Actually, I was able to install v1.5 via the ubuntu repo without difficulty, so I'm able to use the HDF5File writer.

One minor follow-up question - I'm trying to read/write the mesh as part of my routine, but FEniCS complains about my input argument 'mesh'. Is there a recommended way of loading a mesh object after its written to hdf?

Either way, thanks for providing a solution to this problem!!

Pete

      if 1:
    Hdf=HDF5File(mesh.mpi_comm(), "a.h5", "w")
    Hdf.write(mesh, "mesh")
    Hdf.write(x, "x")



  return x

def readSingle():
  #mesh = Mesh()
  mesh = UnitCubeMesh(16,16,16)
  f = HDF5File(mesh.mpi_comm(),'a.h5','r')
  mesh = Mesh()
  f.read(mesh,"mesh")
  V = FunctionSpace(mesh,"CG",1)
  x = Function(V)
  f.read(x,'x')

/////////////////

Traceback (most recent call last):
  File "parallelTest.py", line 162, in <module>
    readSingle()
  File "parallelTest.py", line 98, in readSingle
    f.read(mesh,"mesh")
TypeError: in method 'HDF5File_read', argument 2 of type 'dolfin::MeshValueCollection< bool >

---

That looks like a bug in dolfin, worryingly. I'll check it oout. You can try using xdmf to save and load meshes. That uses h5 under the hood, but has the simple File interface

---

OK, I found the reason for your error message.
You need to do this when reading:

f.read(mesh, "mesh", False)

The reason for the final bool argument is to decide whether or not to reuse the partitioning
information which is saved with the mesh. It is irrelevant when reading in serial, but can speed up the load time when reading back on the same number of nodes as it was saved on.

---

You can also use mpi_comm_world() when you don't have an existing mpi_comm.

---

Many thanks for your help. All works well now!
Pete