Loading an initial condition from a .vtu file

Question

I was wondering for an efficient way of loading the output from a solution as an initial condition for a new problem. Ideally this would work for loading onto meshes with different coarseness, but I haven't found an easy way to load it onto the same mesh.

(edit) It seems that the TimeSeries format does exactly what I'm looking for (when using the same mesh)

Answer 1

You can only load data from XML file. You can easily produce these

File('saved_mesh.xml') << mesh
File('saved_u.xml') << u

If you need to set initial condition on different mesh, you need to load both mesh and data

mesh_old = Mesh('saved_mesh.xml')
V_old = FucntionSpace(mesh_old, ...)
u_old = Function(V_old, 'saved_u.xml')

and project or interpolate to new mesh

mesh_new = ...
V_new = FunctionSpace(mesh_new, ...)
u_new = project(u_old, V_new)

This projection/interpolation from non-matching mesh does not work in parallel in DOLFIN. (There is workaround bypassing this.)

Comments

Great! Thanks a lot for the help. The time_series.retrieve() command was almost there, but it crashes (seg-fault) when loading a mesh. Loading the xml file directly takes care of that.

Regarding projection/interpolation in non-matching meshes, I had read about the issue in parallel, but at the moment I am running on a single process and thought that the parallel features were off by default. I'll give your branch a try.

---

Yes, projection/interpolation on non-matching meshes works fine when run on sigle proccess!

Note that I'm not sure if `serial-meshes` branch compiles with current master of UFL, UFC, FFC and Instant. But it should definitely compile with masters at 31/05/2013. I need to resolve this dependency issue in future - probably by cloning compatible commits of these dependencies to my repo.

Answer 2

If speed and memory is a concern, I suggest you use scipy routines to save and load your data.

I have found that saving and loading a state u via converting u into a numpy array, and using scipy.io routines to save and load the data is much faster than using << and xml files.

Going via scipy appears to be 30 times faster, while consuming only 10% of the memory, necessary for the xml files. These numbers hold for a scalar function on a 512x512 UnitSquareMesh for pw. linear 'CG' functions. The speedup increases with more DOFs.

The code for both procedures is given below. The speed test can be found and inspected on my github account .

from dolfin import UnitSquare, FunctionSpace, Expression, Function, File, project
from scipy.io import loadmat, savemat


mesh = UnitSquareMesh(512, 512)
V = FunctionSpace(mesh, "CG", 2)

u = Expression('x[0]*x[0]')
ufun = project(u, V)

# generic way as Jan Blechta has pointed out
File('saved_u.xml') << ufun
u_old = Function(V, 'saved_u.xml') 

# conversions and scipy io functions
uvec = ufun.vector().array().reshape(len(ufun.vector()),1)
savemat('saved_u', { 'uvec': uvec }, oned_as='column')
uvec_old = loadmat('saved_u')['uvec']
u_old = Function(V)
u_old.vector().set_local(uvec_old[:,0])