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: ornladios/ADIOS2#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]

@ordinary-slim could you elaborate on this logic?
Here you look through all input functions, and replace the first function space with any (the last) function space if it is not a "constant" space. Is this so that you can mix the non-constant and constant spaces?

[ordinary-slim]

This V0 is used to check that all functions have the same dofmap. In the vtk scheme, we can provide both point data and cell data. Since the DG0 functions have this separate path, we only have to verify that the remaining non-DG0 functions use the same dofmap.

[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]

I've made a slight change to ensure that this is working as expected when mixing different function spaces.
Other than that, I think this is good to go. @jhale, shall we merge?:)

[jorgensd]

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

Resolved in 09af8db