This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Time series of many Function instances in a single file

0 votes

My problem requires the output of many state variables - these are mainly concentrations of chemical species in a reactive fluid flow problem.

Ideally, every Function instance for each of these state variables should be output to a single xdmf file.

I have tried the following:

file = XDMFFile(mesh.mpi_comm(), 'result.xdmf')
while t < T:
    for species in species_list:
        file << (n[species], t)
    t += dt

where n[species] returns the Function instance for a species named species (e.g., H2O, CO2).

When loading result.xdmf in Paraview, only the time series for the last chemical species is available - a spreadsheet view of this dataset shows nan for all others.

Having one file for every species requires a lot more effort to setup the visualization in Paraview. Also, a single spreadsheet view of all these data is not easily achievable, which I often need to create plots in gnuplot.

What am I doing wrong?

asked Feb 13, 2015 by allanleal FEniCS Novice (280 points)

Here is the xdmf file for two chemical species (Ca+2, and Mg+2). There are only two time steps.

<?xml version="1.0"?>
<Xdmf Version="2.0" xmlns:xi="http://www.w3.org/2001/XInclude">
  <Domain>
    <Grid Name="TimeSeries_Ca+2" GridType="Collection" CollectionType="Temporal">
      <Time TimeType="List">
        <DataItem Format="XML" Dimensions="2"> 0 100</DataItem>
      </Time>
      <Grid Name="Ca+2_0" GridType="Uniform">
        <Topology NumberOfElements="20" TopologyType="Triangle">
          <DataItem Format="HDF" Dimensions="20 3">data.h5:/Mesh/0/topology</DataItem>
        </Topology>
        <Geometry GeometryType="XY">
          <DataItem Format="HDF" Dimensions="22 2">data.h5:/Mesh/0/coordinates</DataItem>
        </Geometry>
        <Attribute Name="Ca+2" AttributeType="Scalar" Center="Node">
          <DataItem Format="HDF" Dimensions="22 1">data.h5:/VisualisationVector/0</DataItem>
        </Attribute>
      </Grid>
      <Grid Name="Ca+2_2" GridType="Uniform">
        <Topology NumberOfElements="20" TopologyType="Triangle">
          <DataItem Format="HDF" Dimensions="20 3">data.h5:/Mesh/2/topology</DataItem>
        </Topology>
        <Geometry GeometryType="XY">
          <DataItem Format="HDF" Dimensions="22 2">data.h5:/Mesh/2/coordinates</DataItem>
        </Geometry>
        <Attribute Name="Ca+2" AttributeType="Scalar" Center="Node">
          <DataItem Format="HDF" Dimensions="22 1">data.h5:/VisualisationVector/2</DataItem>
        </Attribute>
      </Grid>
    </Grid>
    <Grid Name="TimeSeries_Mg+2" GridType="Collection" CollectionType="Temporal">
      <Time TimeType="List">
        <DataItem Format="XML" Dimensions="2"> 0 100</DataItem>
      </Time>
      <Grid Name="Mg+2_1" GridType="Uniform">
        <Topology NumberOfElements="20" TopologyType="Triangle">
          <DataItem Format="HDF" Dimensions="20 3">data.h5:/Mesh/1/topology</DataItem>
        </Topology>
        <Geometry GeometryType="XY">
          <DataItem Format="HDF" Dimensions="22 2">data.h5:/Mesh/1/coordinates</DataItem>
        </Geometry>
        <Attribute Name="Mg+2" AttributeType="Scalar" Center="Node">
          <DataItem Format="HDF" Dimensions="22 1">data.h5:/VisualisationVector/1</DataItem>
        </Attribute>
      </Grid>
      <Grid Name="Mg+2_3" GridType="Uniform">
        <Topology NumberOfElements="20" TopologyType="Triangle">
          <DataItem Format="HDF" Dimensions="20 3">data.h5:/Mesh/3/topology</DataItem>
        </Topology>
        <Geometry GeometryType="XY">
          <DataItem Format="HDF" Dimensions="22 2">data.h5:/Mesh/3/coordinates</DataItem>
        </Geometry>
        <Attribute Name="Mg+2" AttributeType="Scalar" Center="Node">
          <DataItem Format="HDF" Dimensions="22 1">data.h5:/VisualisationVector/3</DataItem>
        </Attribute>
      </Grid>
    </Grid>
  </Domain>
</Xdmf>

As you can see, there are two time series. You can probably get this file to work by doing a bit of manual editing on it, e.g. move the <Attribute Name="Mg+2"> sections into the <Grid> sections of the Ca+2 section, i.e. something like this:

<?xml version="1.0"?>
<Xdmf Version="2.0" xmlns:xi="http://www.w3.org/2001/XInclude">
  <Domain>
    <Grid Name="TimeSeries_Ca+2" GridType="Collection" CollectionType="Temporal">
      <Time TimeType="List">
        <DataItem Format="XML" Dimensions="2"> 0 100</DataItem>
      </Time>
      <Grid Name="Ca+2_0" GridType="Uniform">
        <Topology NumberOfElements="20" TopologyType="Triangle">
          <DataItem Format="HDF" Dimensions="20 3">data.h5:/Mesh/0/topology</DataItem>
        </Topology>
        <Geometry GeometryType="XY">
          <DataItem Format="HDF" Dimensions="22 2">data.h5:/Mesh/0/coordinates</DataItem>
        </Geometry>
        <Attribute Name="Ca+2" AttributeType="Scalar" Center="Node">
          <DataItem Format="HDF" Dimensions="22 1">data.h5:/VisualisationVector/0</DataItem>
        </Attribute>
        <Attribute Name="Mg+2" AttributeType="Scalar" Center="Node">
          <DataItem Format="HDF" Dimensions="22 1">data.h5:/VisualisationVector/1</DataItem>
        </Attribute>
      </Grid>
      <Grid Name="Ca+2_2" GridType="Uniform">
        <Topology NumberOfElements="20" TopologyType="Triangle">
          <DataItem Format="HDF" Dimensions="20 3">data.h5:/Mesh/2/topology</DataItem>
        </Topology>
        <Geometry GeometryType="XY">
          <DataItem Format="HDF" Dimensions="22 2">data.h5:/Mesh/2/coordinates</DataItem>
        </Geometry>
        <Attribute Name="Ca+2" AttributeType="Scalar" Center="Node">
          <DataItem Format="HDF" Dimensions="22 1">data.h5:/VisualisationVector/2</DataItem>
        </Attribute>
        <Attribute Name="Mg+2" AttributeType="Scalar" Center="Node">
          <DataItem Format="HDF" Dimensions="22 1">data.h5:/VisualisationVector/3</DataItem>
        </Attribute>
      </Grid>
    </Grid>
  </Domain>
</Xdmf>

The problem is now how to define a good user interface for this intended output.

Your solution could work (I haven't tested). However, it is not practical when thousands of lines would need to be moved to the Grid block.

Maybe you can test it, and if it works, the developers can think about how to implement it better. Until then, you will have to use multiple files, I guess.

What you suggested resulted in the exact behaviour I was expecting from Paraview - well done. Now with a single file I can do a Plot Over Line for all species, for example - a time saver. Well, until this is properly fixed, I will try to implement a python script to parse the xdmf file and modify it according to your suggestion.

2 Answers

+1 vote

What you are doing should work, in theory, and I think the .xdmf file will be valid.
It may be a question of how ParaView is interpreting the file.

I think the problem is that each Function is saved as a separate time series (reasonably enough, as you can specify different times "t" for each Function. What ParaView is probably expecting is several sets of values for the same time.

It's not clear if or how this should be implemented in terms of the user interface. Possibly, we could have a way of accepting a list of Functions as an input to XDMFFile.

You could post a simple .xdmf to demonstrate the problem.

answered Feb 13, 2015 by chris_richardson FEniCS Expert (31,740 points)
+1 vote

Based on Chris suggestions, the following python script will convert a time series xdmf file with different Function instances that can be properly understood by Paraview:

The usage is as follows:

./xdmf-parser.py result.xdmf

assuming that the script is called xdmf-parser.py.

#!/usr/bin/python

import sys
import xml.etree.ElementTree as et

if len(sys.argv) < 2:
    raise RuntimeError('No arguments provided to the parser.')

filename = sys.argv[1]

if len(sys.argv) > 2:
    output = sys.argv[2]
else:
    output = filename[:filename.find('.')] + '-out.xdmf'

tree = et.parse(filename)
root = tree.getroot()
domain = root.find('Domain')
grids = domain.findall('Grid')

timeseries = grids[0]

other_timeseries = grids[1:]

attributes = []

for grid in other_timeseries:
    innergrids = grid.findall('Grid')
    attributes_aux = []
    for x in innergrids:
        attribute = x.find('Attribute')
        attributes_aux.append(attribute)
    attributes.append(attributes_aux)

timeseries_grids = timeseries.findall('Grid')

for grid_attributes in attributes:
    for (grid, attribute) in zip(timeseries_grids, grid_attributes):
        grid.append(attribute)

for x in other_timeseries:
    domain.remove(x)

tree.write(output, encoding="utf-8", xml_declaration=True)
answered Feb 13, 2015 by allanleal FEniCS Novice (280 points)

Great. Thanks for letting me know. Perhaps this script could become a utility method in dolfin that would be invoked at the end of the simulation when the file is no longer being written.

That's a great workaround for now, but I don't think it is a logical long-term solution.

Thanks for posting this! I had the same problem and it worked like a charm.

I'll soon post some changes to the development version that will make this workaround unnecessary.

...