Support list of DG-0 + other space in VTXWriter

[ordinary-slim]

Continuation of #3107 so that we're able to output both DG0 + another space in the same bp file. The previous PR by @jorgensd says:

I've decided to keep it such that only single function-space functions should be saved in a single bp file. This could of course be changed now, to store a mixture of DG-0 and a single other element function.

This is up for discussion.

Since VTX is the recommended writer, I think it's necessary to support this feature. This code doesn't modify previous behavior regarding the written post files.

MWE

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

cells_per_side = 4

domain  = mesh.create_unit_square(MPI.COMM_WORLD,
                                  cells_per_side,
                                  cells_per_side,
                                  )

cg1 = fem.functionspace(domain, ("Lagrange", 1),)
dg0 = fem.functionspace(domain, ("Discontinuous Lagrange", 0),)

u_cg1 = fem.Function(cg1, name="u_cg1")
u_dg0 = fem.Function(dg0, name="u_dg0")

expression =  lambda x : x[0]*x[1] 
u_cg1.interpolate(expression)
u_dg0.interpolate(expression)

writer = io.VTXWriter(domain.comm,
                      "result_folder/post.bp",
                      output=[u_cg1,u_dg0])
writer.write(0.0)
writer.close()

[jorgensd]

Currently DG-0 support in VTX is limited by: https://github.com/ornladios/ADIOS2/issues/4179

[ordinary-slim]

Currently DG-0 support in VTX is limited by: ornladios/ADIOS2#4179

Thank you for the information, a bit of a bummer. This PR could still be useful for steady data.

[garth-wells]

We're directing efforts towards VTKHDF5 - it is actively developed/maintained.

[ordinary-slim]

What is the recommended way to output both CG1 and DG0 time dependent data?

[jorgensd]

What is the recommended way to output both CG1 and DG0 time dependent data?

I would interpolate the CG1 and DG0 data into a compatible function space (DG1).

As Paraview 6.0 was just released (https://www.paraview.org/download/) it might be worthwhile getting this PR merged (as VTKHDF currently has very limited visualization capabilities of time dependent data in Paraview). @ordinary-slim do you mind testing if time-dependent data is rendering correctly in paraview 6.0?

[ordinary-slim]

@ordinary-slim do you mind testing if time-dependent data is rendering correctly in paraview 6.0?

I've been using this branch without issues with Paraview 5.13. I can confirm it also works with 6.0 @jorgensd . I am using this example:

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

cells_per_side = 4

domain  = mesh.create_unit_square(MPI.COMM_WORLD,
                                  cells_per_side,
                                  cells_per_side,
                                  )

cg1 = fem.functionspace(domain, ("Lagrange", 1),)
dg0 = fem.functionspace(domain, ("Discontinuous Lagrange", 0),)

u_cg1 = fem.Function(cg1, name="u_cg1")
u_dg0 = fem.Function(dg0, name="u_dg0")

expression =  lambda x : x[0]*x[1] 
u_cg1.interpolate(expression)
u_dg0.interpolate(expression)

writer = io.VTXWriter(domain.comm,
                      "result_folder/post.bp",
                      output=[u_cg1,u_dg0])

for i in range(10):
    time = np.float64(i)
    u_cg1.x.array[:] *= (time + 1.0)
    u_dg0.x.array[:] *= (time + 1.0)
    writer.write(time)

writer.close()
``

[jorgensd]

This doesn't seem to give the right result for me if using a vector-function for the non-constant space.

Minimal test that I added to python/test/unit/io/test_adios2.py:

    def test_dg_0_data(self, tempdir):
        """Test reusage of mesh by VTXWriter."""
        from dolfinx.io import VTXWriter

        adios2 = pytest.importorskip("adios2", minversion="2.10.0")
        if not adios2.is_built_with_mpi:
            pytest.skip("Require adios2 built with MPI support")

        mesh = generate_mesh(2, False)
        v = Function(functionspace(mesh, ("Lagrange", 2, (2, ))))
        filename = Path(tempdir, "v.bp")
        v.name = "v"
        v.interpolate(lambda x: (x[0], x[1]))
        z = Function(v.function_space)
        z.name = "z"
        z.interpolate(lambda x : (np.sin(x[0]), np.cos(x[1])))


        q = Function(functionspace(mesh, ("DG", 0)))
        q.name = "q"
        q.x.array[:] = np.arange(q.x.array.size, dtype=q.x.array.dtype)
 
        # Save three steps
        writer = VTXWriter(mesh.comm, filename, [v,q])
        writer.write(0)
        writer.close()

When not storing q, I get but when storing q it becomes

[jorgensd]

This doesn't seem to give the right result for me if using a vector-function for the non-constant space.

Resolved in https://github.com/FEniCS/dolfinx/pull/3375/commits/09af8dbf38bcff36d165575c2e9a61808f6e7370