How to write binary data with VTK/XDMFFile?

General
1

I would like to write Paraview data in binary to save space. I saw that VTKFile takes a file_mode argument. I tried setting it to wb without success. Looking at the source code I think this argument is unused but I could be wrong. Here is an MVP:

from dolfinx import mesh, io, fem
from mpi4py import MPI

nels_per_side = 2
mesh  = mesh.create_unit_square(MPI.COMM_WORLD, nels_per_side, nels_per_side, mesh.CellType.quadrilateral)

V = fem.functionspace(mesh, ("Lagrange", 1))

u1, u2 = fem.Function(V, name="u1"), fem.Function(V, name="u2")
u1.interpolate( lambda x : x[0]**2 )
u2.interpolate( lambda x : x[1]**2 )

# VTK out
with io.VTKFile(mesh.comm, "out/dolfinx_out.pvd", 'FILE_MODE_CAN_BE_ANYTHING') as iofile:
    iofile.write_function([u1, u2])

This will run and write the data in ASCII format.

2

Use VTXWriter (with engine=“BP4”). As of now VTKFile only supports ascii.

1 Like
3

The issue with VTXWriter is that I can’t mix node and cell data the way VTKFile allows me to. MVP:

from dolfinx import mesh, io, fem
from mpi4py import MPI

nels_per_side = 2
mesh  = mesh.create_unit_square(MPI.COMM_WORLD, nels_per_side, nels_per_side, mesh.CellType.quadrilateral)

V = fem.functionspace(mesh, ("Lagrange", 1))
D = fem.functionspace(mesh, ("DG",0))

u1, u2 = fem.Function(V, name="u1"), fem.Function(D, name="u2")
u1.interpolate( lambda x : x[0]**2 )
u2.interpolate( lambda x : x[1]**2 )

writer = io.VTXWriter(mesh.comm,"test.pvd",[u1,u2])
writer.write(0.0)
writer.close()

This results in:
RuntimeError: All functions in VTXWriter must have the same element type.

What is the best solution, two writers for node and cell data separately? If yes, is there is some way to merge the resulting datasets into one?

4

If you are only using first order Lagrange and DG-0, I would use XDMFFile.

Merging datasets with VTXWriter (.bp) would be trick as of now. I did make a PR for DG-0 fields with VTXWriter: Support DG-0 functions in VTXWriter by jorgensd · Pull Request #3107 · FEniCS/dolfinx · GitHub which in theory (with minor extensions) would make it possible to write Cell-data and function spaces to the same file.

1 Like
5

About this, XDMF is a bit buggy when you write multiple functions. I reported this in this issue. Not only it duplicates the data so that it occupies more memory, the coloring of a dataset can be “hidden” by the other one somehow. See the MWE below:

from dolfinx import mesh, io, fem
from mpi4py import MPI

nels_per_side = 2
mesh  = mesh.create_unit_square(MPI.COMM_WORLD, nels_per_side, nels_per_side, mesh.CellType.quadrilateral)

V = fem.functionspace(mesh, ("Lagrange", 1))
D = fem.functionspace(mesh, ("DG",0))

u1, u2 = fem.Function(V, name="u1"), fem.Function(D, name="u2")
u1.interpolate( lambda x : x[0]**2 )
u2.interpolate( lambda x : x[1]**2 )

writer = io.XDMFFile(mesh.comm,"test.xdmf", 'w')
writer.write_mesh(mesh)
writer.write_function(u1)
writer.write_function(u2)
writer.close()

If I only write u2, the correct number of elements is written and I can see the coloring of u2:

1dset
1dset1243Ă—760 53.4 KB

If I write u1 and u2, a multiblock dataset is written, the elements are duplicated (it takes twice the memory) and I can neither see the coloring of u1 nor u2:

2dsets
2dsets1238Ă—772 62.1 KB

I am using Paraview 5.12. I think there used to be a fix for this in dolfin (see this post). The issue is similar in VTK except for some reason I can still see the coloring of both fields.

Thank you for the update on VTXWriter, it works well. I might give a try to implementing the changes although I am not sure my coding would be up to the standards of your project. From the code in ADIOSWriters.h , it looks like the main tasks are to modify _is_piecewise_constant into an array and modify impl_vtx::extract_function_names to return two vectors of strings for the nodal and DG0 functions repectively.

6

This is not true. If you look at the actual xdmf-file, the data is:

<?xml version="1.0"?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
<Xdmf Version="3.0" xmlns:xi="https://www.w3.org/2001/XInclude">
  <Domain>
    <Grid Name="mesh" GridType="Uniform">
      <Topology TopologyType="Quadrilateral" NumberOfElements="4" NodesPerElement="4">
        <DataItem Dimensions="4 4" NumberType="Int" Format="HDF">test.h5:/Mesh/mesh/topology</DataItem>
      </Topology>
      <Geometry GeometryType="XY">
        <DataItem Dimensions="9 2" Format="HDF">test.h5:/Mesh/mesh/geometry</DataItem>
      </Geometry>
    </Grid>
    <Grid Name="u1" GridType="Collection" CollectionType="Temporal">
      <Grid Name="u1" GridType="Uniform">
        <xi:include xpointer="xpointer(/Xdmf/Domain/Grid[@GridType='Uniform'][1]/*[self::Topology or self::Geometry])" />
        <Time Value="0" />
        <Attribute Name="u1" AttributeType="Scalar" Center="Node">
          <DataItem Dimensions="9 1" Format="HDF">test.h5:/Function/u1/0</DataItem>
        </Attribute>
      </Grid>
    </Grid>
    <Grid Name="u2" GridType="Collection" CollectionType="Temporal">
      <Grid Name="u2" GridType="Uniform">
        <xi:include xpointer="xpointer(/Xdmf/Domain/Grid[@GridType='Uniform'][1]/*[self::Topology or self::Geometry])" />
        <Time Value="0" />
        <Attribute Name="u2" AttributeType="Scalar" Center="Cell">
          <DataItem Dimensions="4 1" Format="HDF">test.h5:/Function/u2/0</DataItem>
        </Attribute>
      </Grid>
    </Grid>
  </Domain>
</Xdmf>

where you ca see that the topology and geometry is only stored once.
Thus it shouldn’t take twice the memory to read in Paraview.

You need to use the Extract Block filter in Paraview to view each individual block in the dataset.

1 Like
7

I’ve added a fix to the VTK duplicate piece issue at: Move piece node outside of function loop in VTKFile by jorgensd · Pull Request #3120 · FEniCS/dolfinx · GitHub

1 Like
8

I tried to follow your steps but I am only shown a single block in the extract block and nothing special happens, see picture below:

extractBlock
extractBlock1132Ă—655 49.3 KB

I am not familiar with XDMF but in VTK, what’s happening is that the VTUs are referenced one more time in the PVTUs for every call to write_function, where they should only be referenced once.

Much appreciated.

9

You need to select the correct block in the top right left corner, see:

image
image3306Ă—1866 241 KB

image
image3306Ă—223 56.3 KB

1 Like
10

Thank you, I managed to get it to work by using a different Paraview reader. Was using Xdmf3_Reader_S and with Xdmf3_Reader_T it works for me even without the extract_block.

I was talking about RAM. In the information window some posts ago, you can see that the dataset occupies 5KBs if I only write a single function and 10KBs if I write 2 functions. I think this can only be the case if the geometry data is read multiple times by Paraview. I think this is confirmed by Paraview displaying 4 and 8 elements in the respective information windows.

11

This is a Paraview issue, not a storage/DOLFINx issue. We provide pointers to the same data in the h5 file for both of them.

As I mentioned earlier, XDMF is no longer maintained by its developers (Kitware) and is not a format that is receiving updates in Paraview. Thus there is very little we will do wrt to this. Therefore I would strongly suggest to use VTXWriter, that uses ADIOS2 (Visualizing Data — ADIOS2 2.10.0-rc1 documentation)

1 Like
12

Thanks for the help and information, much appreciated!