Can't write DG0 time data to VTX

I can’t seem to write the time evolution of a DG0 function with VTXWriter. The values always stay the same for all time steps. This doesn’t occur with other interpolations (I tested CG1 and DG1 so far). Here is a MWE to illustrate what I mean:

from dolfinx import mesh, fem, io
from mpi4py import MPI
import numpy as np

domain = mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)
V = fem.functionspace(domain, ("DG", 0))
v = fem.Function(V)

vtx = io.VTXWriter(domain.comm, "result.bp", v)
vtx.write(0.0)

v.interpolate(lambda x: np.full((x.shape[1],), 1))
vtx.write(1.0)
vtx.close()

I would expect the value of my field to be 0 at time 0 and 1 at time 1, but it is always 0 when I view the .bp in Paraview. Again, using another interpolation the result is as expected. Running v0.8.0 of dolfinx.

At first I thought this was just an issue in DOLFINx, so I modified the writer to add the time step counter into CellData, ref: Comparing main...dokken/time-dependent-celldata · FEniCS/dolfinx · GitHub
Inspecting the results.bp file, with the following MWE i get:

from dolfinx import mesh, fem, io, default_scalar_type
from mpi4py import MPI
import numpy as np

domain = mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)
V = fem.functionspace(domain, ("DG", 0))
v = fem.Function(V)

vtx = io.VTXWriter(domain.comm, "result.bp", v)
vtx.write(0.0)

v.interpolate(lambda x: np.arange(1,1+x.shape[1], dtype=default_scalar_type))
vtx.write(1.0)
vtx.close()
 bpls result.bp/ -a -l
  uint32_t  NumberOfCells        2*{1} = 8 / 8
  uint32_t  NumberOfNodes        2*{1} = 9 / 9
  int64_t   connectivity         2*[1]*{8, 4} = 0 / 8
  double    f                    2*[1]*{8, 1} = 0 / 8
  double    geometry             2*[1]*{9, 3} = 0 / 1
  double    step                 2*scalar = 0 / 1
  uint32_t  types                2*scalar = 69 / 69
  string    vtk.xml              attr   = 
<VTKFile type="UnstructuredGrid" version="0.1">
  <UnstructuredGrid>
    <Piece NumberOfPoints="NumberOfNodes" NumberOfCells="NumberOfCells">
      <Points>
        <DataArray Name="geometry" />
      </Points>
      <Cells>
        <DataArray Name="connectivity" />
        <DataArray Name="types" />
      </Cells>
      <PointData>
        <DataArray Name="vtkOriginalPointIds" />
        <DataArray Name="vtkGhostType" />
      </PointData>
      <CellData>
        <DataArray Name="TIME">step</DataArray>
        <DataArray Name="f" />
      </CellData>
    </Piece>
  </UnstructuredGrid>
</VTKFile>

  uint8_t   vtkGhostType         2*[1]*{9} = 0 / 0
  int64_t   vtkOriginalPointIds  2*[1]*{9} = 0 / 8

which seems correct.
Next, I checked if the data is correctly written to f with

bpls result.bp/ -a -l f -d
  double    f                    2*[1]*{8, 1} = 0 / 8
        step 0: 
          block 0: [0:7, 0:0] = 0 / 0
    (0,0)    0 0 0 0 0 0
    (6,0)    0 0 
        step 1: 
          block 0: [0:7, 0:0] = 1 / 8
    (0,0)    1 2 3 4 5 6
    (6,0)    7 8 

which looks correct. However opening it in paraview, I cannot step between the two time steps:


As you can see in data statistics, it realizes that there are two time steps,l but in data arrays we can observe that f is always 0. This seems to be an issue with Paraview’s VTXWriter.
I’ve made an issue at: Time dependent VTX CellData not rendered correctly in Paraview · Issue #4179 · ornladios/ADIOS2 · GitHub which might be transferred to VTK depending on what the developers says.

1 Like

Thank you for the quick look at the issue. Any apparent reason why this only happens with DG0 interpolation, though?

Its because of how paraview/VTXReader reads cell data. In the original VTX documentation, there was only time dependent point data, no cell data, so I guess this is an oversight on their end.

1 Like

Confirmed as an issue by the developers: Time dependent VTX CellData not rendered correctly in Paraview · Issue #4179 · ornladios/ADIOS2 · GitHub
We cant do much more than wait and see

1 Like

FYI, I tested it out with the nightly binaries that Paraview offers on their website and it was fixed there.

1 Like