Recommended way to combine functions from different checkpoint files based on the same mesh

[stephankramer]

My question is basically how to best fix checkpointing for Thetis (see also issue thetisproject/thetis#336)

The current (broken) pattern is that the user provides a mesh, Thetis creates a solver object with various function spaces and necessary functions, and also the user may specify various additional user created functions, all based on that same mesh/domain, and then after that the user can decide to apply some user defined initial condition, or load the initial condition from a checkpoint created in a previous run. Currently loading from the checkpoint involves loading the function from the checkpoint file using the mesh that's stored in the checkpoint file (see https://github.com/thetisproject/thetis/blob/master/thetis/exporter.py#L222) and then we assign its values to the already previously created function which was based on the user supplied mesh. That last step is currently broken as firedrake (ufl?) has gotten stricter on not allowing assignment between functions on (what it thinks are) different domains.

One hacky solution would be to simply copy .dat.data[:] values instead of calling assign. Is this safe in parallel? What about if I restart on a different number of cores?

Alternatively we would have to force the user to read the mesh from one of the checkpoint files (possibly with some helper function) and then build on top of that. Note that there are typically multiple checkpointfiles involved, each with different functions that we want to use: say some preprocessed bathymetry, then an initial condition for velocity and an initial condition for elevation. But am I right in saying that we are allowed to load the mesh from one checkpoint file A, and then load functions from another checkpoint file B using the mesh from checkpoint A instead of loading the mesh from B?

[pefarrell]

Not an answer, but in defcon I do function_to_return = Function(function_space_on_mesh_I_want, val=function_from_checkpoint.dat):

https://bitbucket.org/pefarrell/defcon/src/dc77cbc7d321240b122dfc212364bf08c1631433/defcon/iomodule.py#lines-187

[ksagiyam]

I think you are right. It is expected that we can load a mesh from fileA, using which we can load a function from fileB. The following simple example works.

from firedrake import *
import numpy as np


mesh = UnitSquareMesh(2, 2, quadrilateral=True, name="m")
V = FunctionSpace(mesh, "CG", 3)
x, y = SpatialCoordinate(mesh)
f = Function(V, name="f").interpolate(exp(x + 2 * y))
with CheckpointFile("fileA.h5", "w") as cf:
    cf.save_mesh(mesh)
with CheckpointFile("fileB.h5", "w") as cf:
    cf.save_function(f)
with CheckpointFile("fileA.h5", "r") as cf:
    mesh = cf.load_mesh("m")
with CheckpointFile("fileB.h5", "r") as cf:
    g = cf.load_function(mesh, "f")
assert np.allclose(f.dat.data, g.dat.data)

With this in mind, I think, the alternative approach should be used as we do not want to rely on that the user has provided correct mesh with correct distribution/overlap and that the user has constructed correct function spaces on the mesh.