Error in paraview visualization of xdmf-file with multiple solution variables

I try to solve a problem where I end up with multiple solution variables. xdmf output of more than one variable can’t be visualized in paraview (5.9.1), pvd output works. I’m using the latest docker image of dolfinx.
Here is a minimal example code producing a xdmf file with the error:

import dolfinx
import numpy as np
import ufl
from dolfinx import (DirichletBC, Function, FunctionSpace, RectangleMesh, fem)
from dolfinx.cpp.mesh import CellType
from dolfinx.fem import locate_dofs_topological
from import XDMFFile
from dolfinx.mesh import locate_entities_boundary
from mpi4py import MPI
from ufl import ds, dx, grad, inner

# Create mesh and define function space
mesh = RectangleMesh(
    [np.array([0, 0, 0]), np.array([1, 1, 0])], [32, 32],
    CellType.triangle, dolfinx.cpp.mesh.GhostMode.none)

V = FunctionSpace(mesh, ("Lagrange", 1))

# Define boundary condition on x = 0 or x = 1
u0 = Function(V)
with u0.vector.localForm() as u0_loc:
facets = locate_entities_boundary(mesh, 1,
                                  lambda x: np.logical_or(np.isclose(x[0], 0.0),
                                                          np.isclose(x[0], 1.0)))
bc = DirichletBC(u0, locate_dofs_topological(V, 1, facets))

# Define variational problems
x = ufl.SpatialCoordinate(mesh)
f = 10 * ufl.exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
g = ufl.sin(5 * x[0])

# first problem
u0, v0 = ufl.TrialFunction(V), ufl.TestFunction(V)
a0 = inner(grad(u0), grad(v0)) * dx
L0 = inner(f, v0) * dx + inner(g, v0) * ds

problem0 = fem.LinearProblem(a0, L0, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh0 = problem0.solve()

# second problem
u1, v1 = ufl.TrialFunction(V), ufl.TestFunction(V)
a1 = 2 * inner(grad(u1), grad(v1)) * dx
L1 = inner(f, v1) * dx + inner(g, v1) * ds

problem1 = fem.LinearProblem(a1, L1, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh1 = problem1.solve()

# write xdmf file
xdmf = XDMFFile(MPI.COMM_WORLD, "poisson.xdmf", "w")


I’m not sure if that’s a problem with paraview or the xdmf file, the visualization in paraview looks like this:

When I’m running the channel flow example, I end up with the same error.

You need to use the extract block filter, as explained in Using Paraview for visualization — FEniCSx tutorial


Thank you very much.