Write a function to HDF5 in parallel - garbled data when read back in?

Question

I'm trying to write a function to an hdf5 file viz Johan's example in:
here

When I try this in 1.4.0 I still get garbled data if I run the first program with more than one core.
With the first with one core, I can run the second, reading program with as many cores as I want.

(Code copied from above link)

from dolfin import *
import cPickle

mesh = UnitSquareMesh(12,12)
V = FunctionSpace(mesh, "CG", 1)
e = Expression("x[0]*(t*0.1+1)", t=0., degree=1)
u = project(e, V)
f = HDF5File(mesh.mpi_comm(), "u.h5", "w")
f.write(u, "/initial")
times = []
t = 0; dt=0.5
while t < 5:
    t += dt
    e.t = t
    u.interpolate(e)
    times.append(t)
    f.write(u.vector(), "/values_{}".format(len(times)))
cPickle.dump(times, open("times.cpickle", "w"))

Then when I load the values I do something like:

from dolfin import *
import cPickle

mesh = UnitSquareMesh(12,12)
V = FunctionSpace(mesh, "CG", 1)
u = Function(V)

times  = cPickle.load(open("times.cpickle"))
f = HDF5File(mesh.mpi_comm(), "u.h5", "r")
f.read(u, "/initial")
plot(u, title="u_{}".format(0))
for ind, t in enumerate(times):
    f.read(u.vector(), "/values_{}".format(ind+1), True)
    plot(u, title="u_{}".format(t))

Answer 1

Doing it this way assumes the mapping from Vector to dof does not change.
Try using the hdf5 function read.In your example, for instance, the initial Function should be ok

Comments

Thanks Chris, you are right, the initial function does look correct. I tried your suggestion on the same link (writing the function to the file with a time index) but couldn't work out how to access other time steps. Specifically, my attempt to read the 3rd result vector with:

f.read(u, "/test/vector_3")

came up with an error:

RuntimeError: *** Error: Group with name "/test/vector_3" does not exist

My h5ls -r returns:

/                        Group
/mesh                    Group
/mesh/cell_indices       Dataset {288}
/mesh/coordinates        Dataset {169, 2}
/mesh/topology           Dataset {288, 3}
/test                    Group
/test/cell_dofs          Dataset {864}
/test/cells              Dataset {288}
/test/vector             Dataset {169}
/test/vector_1           Dataset {169}
/test/vector_10          Dataset {169}
/test/vector_2           Dataset {169}
/test/vector_3           Dataset {169}
/test/vector_4           Dataset {169}
/test/vector_5           Dataset {169}
/test/vector_6           Dataset {169}
/test/vector_7           Dataset {169}
/test/vector_8           Dataset {169}
/test/vector_9           Dataset {169}
/test/x_cell_dofs        Dataset {289}

It seems like in trying to read to a function, it is looking for a group with that name (and only getting a vector). As is, it seems like my best option is to combine the techniques and write new functions for each timestep. Unless you have another suggestion?
Thanks

---

For those interested, my current working files are:

from dolfin import *
import cPickle

mesh = UnitSquareMesh(12,12)
V = FunctionSpace(mesh, "CG", 1)
e = Expression("x[0]*(t*0.1+1)", t=0., degree=1)
u = project(e, V)
f = HDF5File(mesh.mpi_comm(), "u.h5", "w")
f.write(u, "/initial")
times = []
t = 0; dt=0.5
while t < 5:
    t += dt
    e.t = t
    u.interpolate(e)
    times.append(t)
    f.write(u, "/values_{}".format(len(times)))
cPickle.dump(times, open("times.cpickle", "w"))

and

from dolfin import *
import cPickle

mesh = UnitSquareMesh(12,12)
V = FunctionSpace(mesh, "CG", 1)
u = Function(V)

times  = cPickle.load(open("times.cpickle"))
f = HDF5File(mesh.mpi_comm(), "u.h5", "r")
f.read(u, "/initial")
file = File('test.pvd')
file << (u,0.)
#plot(u, title="u_{}".format(0))
for ind, t in enumerate(times):
    f.read(u, "/values_{}".format(ind+1))
    file << (u,t)
    #plot(u, title="u_{}".format(t))

The increase in file size is 34kB (vectors only) vs 134kB (functions: current examples). Works as a stopgap on multiple cores writing and reading.

---

I thought the functionality of reading a vector as a Function (as above) worked, but as it turns out, there is a typo in the source code... I'll raise an issue to get it fixed for the next release...

Answer 2

You could consider trying cbcpost (that's what I ended up doing):

http://cbcpost.readthedocs.org/en/latest/