Loading xdmf data back in

This may not strictly be a fenics question, but suppose I run a time dependent problem, of the sort:

output = XDMFFile("out.xdmf")

for n in range(nmax):
# computational steps
    output.write(u,t)
output.close()

If I then read this in with paraview, it’s all as expected. But I was wondering how to read the data back in with python. I saw the post here, but I don’t quite see how to adapt it to what I’m doing.

You should use write_checkpoint and read_checkpoint as described in the previous post to read things back into FEniCS.

Two follow up questions:

  1. If I am fine reading the data in with Paraview (or some similar program), is there any reason I should stop using write?
  2. I am still struggling to get the suggested approach working. Here is what a minimal example:
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

alpha = 3; beta = 1.2
f = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t', degree=2, alpha=alpha, beta=beta, t=0)

f_out = XDMFFile("test.xdmf")

for j in range(5):
    t = float(j/5)
    f.t = t
    f_out.write_checkpoint(project(f,V), "f",t)

f_out.close()

f1 = Function(V)
f_in =  XDMFFile("test.xdmf")

f_in.read_checkpoint(f1,"f",0) # this works fine
f_in.read_checkpoint(f1,"f",1) # this throws an error

The error that I get is:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-48-ab851acfbb36> in <module>
----> 1 f_in.read_checkpoint(f1,"f",1)

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to reading data from XDMF file.
*** Reason:  Storage format "" is unknown.
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  
*** -------------------------------------------------------------------------

write saves the data at the vertices, efficiently projecting the solution to a CG-1 function space. write_checkpoint saves the data with the appropriate function space, and can therefore be read into dolfin again, to the appropriate function space.
You are not saving multiple time steps to your file as you have not changed the append kwarg in write_checkpoint (see: https://fenicsproject.org/docs/dolfin/latest/cpp/de/d71/classdolfin_1_1XDMFFile.html#a9ef9031aec7bfec8f289d5210a1c811c). The following works:

from dolfin import *
nx = ny = 8 
mesh = UnitSquareMesh(nx, ny) 
V = FunctionSpace(mesh, 'P', 1) 
    
alpha = 3; beta = 1.2 
f = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t', degree=2, alpha=alpha, beta=beta, t=0) 
 
f_out = XDMFFile("test.xdmf") 
f_out.write_checkpoint(project(f, V), "f", 0, XDMFFile.Encoding.HDF5, False)  # Not appending to file
for j in range(1,5): 
    t = float(j/5) 
    f.t = t 
f_out.write_checkpoint(project(f,V), "f",t, XDMFFile.Encoding.HDF5, True)  #appending to file
    
f_out.close() 
     
f1 = Function(V) 
f_in =  XDMFFile("test.xdmf") 
   
f_in.read_checkpoint(f1,"f",0) 
f_in.read_checkpoint(f1,"f",1)   
2 Likes

Yes, that now works. One question: output in generated with write_checkpoint now doesn’t seem to be readable by Paraview. Is this to be expected, or is there some way to make it compatible with both?

then your paraview is out of data. You Need to use Paraview >=5.5.0

is there a way with the method read_checkpoint to read the actual time of the i_th time step ? let’s say for example you want to recreate the exact same xdmf file ?

from fenics import *

file_in = XDMFFile("in.xdmf")
file_out = XDMFFile("out.xdmf")

mesh = Mesh()
file_in.read(mesh)

V = FunctionSpace(mesh, "P", 1)
u = Function(V)

for i in range(0, 10):
    file_in.read_checkpoint(u, "u", i)
    t = i # ?????
    file_out.write_checkpoint(u, "u", i)
    export_file.write_checkpoint(u, "u", t)

As you can see from the documentation: https://fenicsproject.org/docs/dolfin/latest/cpp/de/d71/classdolfin_1_1XDMFFile.html#aa86a07b92188a2da32d1a0e128106c91
That is not possible, as read_checpoint only takes in a counter.

I’m using 5.7.0, and I’m seeing that some of the fields are listed as partial. The consequence seems to be that I don’t have the field at every time at which I performed a write_checkpoint. See attached screen capture.

31%20PM

Are you using the append functionality. I see none of the issues with the following code:

from dolfin import *

nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

alpha = 3; beta = 1.2
f = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t', degree=2, alpha=alpha, beta=beta, t=0)
f_out = XDMFFile("test.xdmf")

f_out.write_checkpoint(project(f,V), "f", 0.0,  XDMFFile.Encoding.HDF5, False)

for j in range(1,5):
    t = float(j/5)
    f.t = t
   
    f_out.write_checkpoint(project(f,V), "f", t,  XDMFFile.Encoding.HDF5, True)

f_out.close()

f1 = Function(V)
f_in =  XDMFFile("test.xdmf")

f_in.read_checkpoint(f1,"f",0) 
f_in.read_checkpoint(f1,"f",1) 
f_in.close()
2 Likes

Yeah, my mistake is I was appending to an existing file. Once I put False in the first call, it cleared up.

Hello Jorgen, is there still any way to read the data back in when they were already written using write instead of write_checkpoint?

No, you cannot read in data written with XDMFFile.write.