Saving multiple functions to a File gives incorrect results without any error message

[Corrado Maurini]

It is not clear to me if saving multiple functions to a file through the File interface is intended to be supported. However the current behavior is not acceptable.

The following example shows that one can save multiple functions to a file without getting any error message. However, when trying to visualize the result with Paraview one gets that:

• with pvd format all the fields appear in the Paraview menu but only the last saved is non-null.

• with xdmf format one gets an error message when dying to open the file in Paraview.

from dolfin import *

mesh = UnitSquareMesh(4,4)

# Create named functions 
U = FunctionSpace(mesh, "CG", 1)
u1 = Function(U, name = "u1")
u2 = Function(U, name = "u2")
u1.interpolate(Expression("x[0]"))
u2.interpolate(Expression("x[1]"))

# Save to file  
file_xdmf = File("my_file.xdmf")
file_pvd = File("my_file.pvd")
file_xdmf << (u1,0.)
file_xdmf << (u2,0.)
file_pvd << (u1,0.)
file_pvd << (u2,0.)

##
## Uncomment to test also vector-valued functions
##
#V = VectorFunctionSpace(mesh, "CG", 1)
#v = Function(V, name = "v") 
#v.interpolate(Expression(("x[0]","x[1]")))  
#file_xdmf << (v,0.)
#file_xdmf << (v,0.) 
#file_pvd << (v,0.)  

del file_pvd
del file_xdmf 

[Chris Richardson]

It is probably possible to divert each Function into a different block of XML in the XDMF format, e.g. by putting the Function name as part of the Grid Name attribute for the TimeSeries.

[Chris Richardson]

New functionality now in branch chris/xdmf-multiple-function

[Chris Richardson]

XDMF now supports multiple functions per file, though it is not clear if ParaView accepts them correctly. More testing is probably needed.

PVD file support will be more difficult. I'm going to postpone any more changes to 1.6

[Tianyi Li]

Using the example given by @cmaurini, the XDMF output obtained gives in ParaView a correct u2(x,y)=y field but produces an uniform u1(x,y)=1. Same case for the PVD format.

[Chris Richardson]

@tianyikillua - my experience is that it is possible to get ParaView to display both fields, but it is not easy. You have to fiddle around with the Properties panel, under "Point Arrays" and "Hierarchy" and "Coloring". If you can find a flaw in the produced XDMF and suggest a fix, I'm happy to include it...

To be honest, I think it is better just to keep each field in a separate file.

[Chris Richardson]

Another example of this problem.

It is not clear what a suitable user interface should look like. In most cases, what is required is a way to save multiple Functions at the same timestep.

We could implement XDMFFile::write(std::vector<std::shared_ptr<const Function> >, double time) I suppose.

[Allan Leal]

A python script can be found here, which will parse the xdmf file and transform it so that Paraview can properly understand its contents.

I would say this should be only a temporary solution. Perhaps that script could be transformed into a utility method in dolfin.

As to the interface, I would stick with the practice:

file << (u1, t)
file << (u2, t)

because u2 above can in fact be the same Function u1, but with updated values. If vector<const Function*> is required, then this is no longer possible to do.

[Chris Richardson]

@allanmulin - I'm not sure what you mean about the interface. My suggestion would result in:

file << ([u1, u2, u3], 0.0)
file << ([u1, u2, u3], 1.0)

etc. The problem at the moment is that there is no guarantee of two functions having the same set of times in their time series.

[Allan Leal]

@chris_richardson Consider the example:

u = Function(U)

while t < T:
  for i in range(10):
      u.interpolate(Expression("i*x[0] + x[1]/i", i=i))
      u.rename('u' + str(i))
      file << (u, t)

Note that only one instance of Function is necessary. Now, if one has to do file << ([u1, u2, u3, ..., u10], t), then this becomes a burden to the user - in that specific case, for example, there would be 9 unnecessary Function instances to manage.

I believe that reading one pair (Function, double) at a time like in file << (u, t) is the most flexible solution.

[Prof Garth Wells]

We should stop using the << syntax. Functions are not streams, and we should be able to stick any data we want in a HDF5 file (with a tag) and pull it out with a tag.

[Jan Blechta]

Generally agree with Garth. But I'd remind that a tag could be derived from a combination of object name/id and/or time it it helps anywhere.

[Allan Leal]

@chris_richardson, what about just modifying the destructor of class XDMFFile? You keep everything else as it is for now, and when the file is no longer going to be written, you apply this python script to the created file xdmf file. One would only need to use the C++ library pugixml instead of the python library ElementTree.

This way, we can save one function at a time (with a tag if needed) like how I described earlier.

As to @garth-wells, a double, say 3.14, isn't a stream like a Function, and still it is quite obvious and standard what std::cout << 3.14; means. I find the functionality provided by operator<< in class File quite intuitive. However, I would not see a problem to use file.write(u, tag) and file.read(u, tag).

[Chris Richardson]

Thanks for the feedback. I'm sure we'll find a more suitable fix for the next release.

[Prof Garth Wells]

@allanleal It's not that Function is a stream either. The design of File and File::operator<< is flawed, which we have to live with now but we shouldn't perpetuate the design mistake, because we end up in even more of a knot.

HDF5 (and to a lesser extend XDMF) is a very powerful format and shouldn't go through the File interface. The File interface 'works' if all file backend have the same functionality and are suited to the same interface, which they aren't. The File interface is one of the reasons DOLFIN has never properly supported hierarchical file formats, despite some formats being designed to do this. e.g. XML and HDF5.

[Martin Sandve Alnæs]

What Garth said. Let's not use << for anything else than legacy simple writing of an object to file.

[Nico Schlömer]

What is the currently recommended way for storing data then? (Most (all?) examples from documtation use File.)

[Prof Garth Wells]

@nschloe You can use File. In developing DOLFIN, we just shouldn't make the File interface 'support' features that are not appropriate for the interface.

[Tormod Landet]

Just a note about a previous comment:

Using the example given by Corrado Maurini, the XDMF output obtained gives in ParaView a correct u2(x,y)=y field but produces an uniform u1(x,y)=1. Same case for the PVD format.

What is happening is that the pvd file points to multiple vtu files per time step. Since these vtu files all contain the mesh there will be multiple meshes per time step. If you see a uniform field then you are looking at a mesh belonging to another function. You can right click on the mesh and select "Hide Block" to see what is below. You can play with "Hide Block" and "Show All Blocks" to see all fields. This is inconvenient, but it may be how it must be unless multiple functions end up on the same vtu file.

[Chris Richardson]

Using master branch, the XDMF now produces the correct output, and uses the same mesh for both u1 and u2 (at least, it does for me with the code above).

[Tormod Landet]

I can confirm this. XDMF works, PVD does not (without hiding and showing meshes).

[Chris Richardson]

Closing this, as it now broadly works for XDMF. If PVD is an issue, then it can be reopened as a new issue, but I suggest deprecating PVD where possible.

[Prof Garth Wells]

Does PVD throw an error?

[Chris Richardson]

I guess it doesn't throw an error, but the usage is a bit strange, anyway. This is could be considered a symptom of <<.

[Prof Garth Wells]

Could it throw an error, or does the interface prevent this? First write could store the reference to the Function.