How to use TimeSeriesHDF to store and retrieve functions in parallel?

Question

I am running a time-dependent simulation on M processors in parallel and would like to store my simulation results in each timestep in the HDF5 format. Subsequently, I would like to retrieve the results in a different script on N processors, typically N = 1. The TimeSeriesHDF5 class looks suitable for this, but at the moment my results seem to get scrambled upon either storing or retrieval.

My questions are
1) Is this supposed to work?
2) If yes, how?
3) If no, what is the preferred alternative approach for storing simulation results and eventually retrieving them as Functions?

Please see below for a minimal example. Consider the following two scripts: store.py:

from dolfin import *

mesh = UnitSquareMesh(12,12)
V = FunctionSpace(mesh, "CG", 1)
f = Expression("x[0]", degree=1)
f = project(f, V)

series = TimeSeriesHDF5(mesh.mpi_comm(), "f")
for i in range(1):
    series.store(f.vector(), i)
    plot(f, title="f_%g" % i)

and retrieve.py:

from dolfin import *

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

series = TimeSeriesHDF5(mesh.mpi_comm(), "f")
times = series.vector_times()
for t in times:
    series.retrieve(f.vector(), t, False)
    plot(f, title="f_%g" % t)

Running

mpirun -n 2 python store.py
python retrieve.py

results in the resulting retrieved f looking scrambled.

My apologies if this question has been answered multiple times already, I couldn't find a detailed answer here in the Q&A.

Comments

To me, the plot looks exactly the same in both cases.

Sorry, I didn't run the correct commands.

Answer 1

I use a combination of HDF5File and pickle to store a time series in DOLFIN.

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 2

TimeSeriesHDF5 is not really suitable for this, as it only stores Vector and Mesh, but no dofmap information, which is essential to do what you require.

You can use HDF5File to store Functions directly, including time series. You can either write each timestep as a separate data item, or I think it also accepts series if you repeatedly write to the same name.

u = Function(Q)
hdf = HDF5File(mpi_comm, 'file.h5','w')
hdf.write(u, 'velocity', 0.0)

hdf.write(u, 'velocity', 0.1)

To read it back, you have to do something like:

hdf = HDF5File(mpi_comm, 'file.h5','r')
hdf.read(u, 'velocity/vector_1')

etc. You may need to use h5ls to check inside your .h5 file

Comments

As far as I can see, the primary use case for TimeSeries must be to store a function changing over time. As the current functionality in TimeSeriesHDF5 is rather confusing, perhaps it could be removed?

---

Probably should be addressed on the mailing list or bitbucket