Read xdmf file to do pyvista visualizations

[renol-kth]

Dear community,

similar questions on the forum

there have been several questions on how to read xdmf files (that store fenicsx functions and meshes) in order to visualize them in pyvista. With this I mean creating a plot of the mesh or function (see e.g. here); or creating a gif of a time-dependent function (e.g. here).

• Read exported solution back into fenicsx about how to do the above in parallel. It is suggested to use the adios package. Can I achieve my goal without this package in certain cases - e.g. if I am fine without parallel computing?
• there are discussions about this for legacy fenicsx, e.g. Loading xdmf data back in - #4 by dokken or Access function in time dependent XDMFFile - #6 by RemDelaporteMathurin or Save solution and load it in the next iteration
• somebody wanted to read only the last time-step of a time-dependent function from an xdmf file, here: Read solution of last time step from file but this post was created before the publication of the adios package, so it also does not answer my question.

Saving mesh and function data in memory-efficient ways (i.e. by only saving raw function and mesh values) is important for me to be able to reproduce plots and gifs efficiently, without recomputing solutions.

my appraoch fails

My naive approach is to do the following:

from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import fem
from ufl import dx, grad, inner
from basix.ufl import element
from dolfinx.fem import dirichletbc, functionspace, locate_dofs_topological
from dolfinx.mesh import create_unit_square, locate_entities_boundary

# mesh, spaces, BC
msh = create_unit_square(MPI.COMM_WORLD, 10,10)
P1 = element("Lagrange", "triangle", 1) 
XV = functionspace(msh, P1)
facetsV = locate_entities_boundary(msh, dim=1,
                                marker=lambda x: np.logical_or.reduce((
                                    np.isclose(x[0], 0.0),
                                    np.isclose(x[0], 1.0),
                                    np.isclose(x[1], 0.0),
                                    np.isclose(x[1], 1.0))))
dofsV = locate_dofs_topological(XV, entity_dim=1, entities=facetsV)
BCs = [dirichletbc(0.0,dofsV, XV)]

# variational problem
u = ufl.TrialFunction(XV)
v = ufl.TestFunction(XV)
x = ufl.SpatialCoordinate(msh)
f = 1
a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx
problem = LinearProblem(a, L, bcs=BCs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

# xdmf file
from dolfinx import io
xdmf = io.XDMFFile(msh.comm, "test.xdmf", "w")  # Open XDMF file for writing
xdmf.write_mesh(msh)
xdmf.write_function(uh) 
xdmf.close()

# read test.xdmf
import pyvista as pv
reader = pv.get_reader("test.xdmf")
grid = reader.read()

# plot using pyvista
p = pv.Plotter()
p.add_mesh(grid)
p.show()

but the code (in particular the line grid = reader.read() ) returns

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run 
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0

which I cant explain, as .read() seems to be well-implemented, see here. A similar approach for a time-dependent problem (in order to create a snapshot or gif of a function) fails as well.

This means, that at the moment I can only produce visualizations directly inside the same file that also produces the function i want to visualize.

my question:

Is the best practice to use adios, or can it be avoided in simpler cases (serial, Lagrange elements, etc)?

Or do professional fenicsx users - who have a similar workflow consisting of saving solutions as raw data in order to create visualisations later on - stick to paraview?

[dokken]

If the only purpose of the output/input is to read it into pyvista, you could store the data in a pyvista-friendly format, rather than going through DOLFINx.

I.e. make a pyvista unstructured grid from the data and store it with pyvista.

As noted in all Pyvista-visualization, they have no concept of other finite elements than arbitrary order Lagrange.

Note that xdmf is a fairly generic format. This means that even if a mesh has been written to a Paraview compatible xdmf format, it is not always compatible with other formats (an example is trying to use meshio to read data from a xdmf file with two meshes in it).

[renol-kth]

Thank you for your answer.

1. Is it true, that there is no way to produce an xdmf-based workflow (essentially my code above) that

• writes from dolfinx to xdmf
• reads from xdmf to pyvista
• plots with pyvista
  ?

2. Concerning your suggestion to use other file formats - which ones?
   I happen to know that .vtk files were used in Legacy-fenics. Are you recommending to use these ones?
   I remember, that they used something like the following:

# Poisson problem
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import fem
from ufl import dx, grad, inner
from basix.ufl import element
from dolfinx.fem import dirichletbc, functionspace, locate_dofs_topological
from dolfinx.mesh import create_unit_square, locate_entities_boundary

# mesh, spaces, BC
msh = create_unit_square(MPI.COMM_WORLD, 10,10)
P1 = element("Lagrange", "triangle", 1) 
XV = functionspace(msh, P1)
facetsV = locate_entities_boundary(msh, dim=1,
                                marker=lambda x: np.logical_or.reduce((
                                    np.isclose(x[0], 0.0),
                                    np.isclose(x[0], 1.0),
                                    np.isclose(x[1], 0.0),
                                    np.isclose(x[1], 1.0))))
dofsV = locate_dofs_topological(XV, entity_dim=1, entities=facetsV)
BCs = [dirichletbc(0.0,dofsV, XV)]

# variational problem
u = ufl.TrialFunction(XV)
v = ufl.TestFunction(XV)
x = ufl.SpatialCoordinate(msh)
f = 1
a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx
problem = LinearProblem(a, L, bcs=BCs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

# write solution as vtk/pvd file
vtk = io.VTKFile(msh.comm, "out/dolfinx_out.pvd", 'w')
vtk.write_mesh(msh)
vtk.write_function(uh)

# read solution from vtk file with pyvista
import pyvista as pv
reader = pv.get_reader("out/dolfinx_out.pvd")
mesh = reader.read()

# plot
plotter = pv.Plotter()
plotter.add_mesh(mesh)
plotter.show()

However the plot shows an empty domain (or the zero function?), while the function is clearly not zero everywhere.

3. Concerning your comment about xdmf: Are you suggesting, that for more complicated element types, one should use xdmf together with paraview?

[dokken]

@renol-kth

1. Is it true, that there is no way to produce an xdmf-based workflow (essentially my code above) that

• writes from dolfinx to xdmf
• reads from xdmf to pyvista
• plots with pyvista

I don’t know the internal details of the pyvista xdmf reader, so I cannot say if that is true or not.

What I am suggesting is that you go from DOLFINx->Pyvista->whatever format pyvista wants to store data as.

@renol-kth

Concerning your suggestion to use other file formats - which ones?
I happen to know that .vtk files were used in Legacy-fenics. Are you recommending to use these ones?

I’ve written a summary of Paraview-supported formats at: Visualization and outputting formats — FEniCS Workshop

You could alternatively use adios4dolfinx to create a checkpoint (readable function) for DOLFINx, which means that you can read it in at a later stage, and then use the dolfinx.plot.vtk_mesh function to convert the mesh to pyvista data.

[renol-kth]

Thank you for your answer.

I managed to create a VTK-based workflow, which works well for static functions. For time-dependent functions, xdmf is much easier to use and does store the mesh only once, which is nice.

The option to read xdmf in pyvista was requested in this Github discussion, and then later implemented (see docs).

However i have not managed to make it work. Here is an MWE (based on the heat-equation example from the tutorial):

# HEAT EQUATION

import pyvista
import ufl
import numpy as np
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, io, plot
from dolfinx.fem import locate_dofs_topological,dirichletbc
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.mesh import *```

# sim parameters
t = 0  
T = 0.1 
num_steps = 20
dt = T / num_steps

# mesh and functions
domain = create_unit_square(MPI.COMM_WORLD, 10,10)
V = fem.FunctionSpace(domain, ("Lagrange", 1))
u_n, uh = fem.Function(V), fem.Function(V)
def initial_condition(x, a=5): return np.exp(-a * (x[0]**2 + x[1]**2 + x[2]**2))
u_n.interpolate(initial_condition)
uh.interpolate(initial_condition)

# BC
facets = locate_entities_boundary(domain, dim=1,
                                    marker=lambda x: np.logical_or.reduce((
                                        np.isclose(x[0], 0.0),
                                        np.isclose(x[0], 1.0),
                                        np.isclose(x[1], 0.0),
                                        np.isclose(x[1], 1.0))))
dofsV = locate_dofs_topological(V, entity_dim=1, entities=facets)
bc = dirichletbc(0.0,dofsV, V)

# problem
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))
a = u * v * ufl.dx + dt * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = (u_n + dt * f) * v * ufl.dx
bilinear_form = fem.form(a)
linear_form = fem.form(L)
A = assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = create_vector(linear_form)

# solver
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

# xdmf file
from dolfinx import io
xdmf = io.XDMFFile(domain.comm, "test.xdmf", "w")
xdmf.write_mesh(domain)

for i in range(num_steps):
    t += dt

    # update right hand side, apply BC
    with b.localForm() as loc_b:
        loc_b.set(0)
    assemble_vector(b, linear_form)
    apply_lifting(b, [bilinear_form], [[bc]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b, [bc])

    # Solve and update sol
    solver.solve(b, uh.vector)
    uh.x.scatter_forward()
    u_n.x.array[:] = uh.x.array

    # write solution to xdmf
    xdmf.write_function(uh, t) 

xdmf.close()

After running this and producing test.xdmf, I would like to read it in pyvista. However, neither

import pyvista as pv
xdmf_reader = pv.XdmfReader("test.xdmf")
mesh = xdmf_reader.read()
mesh.plot()

nor

import pyvista as pv
reader = pv.get_reader("test.xdmf")
mesh = reader.read()
mesh.plot()

work. They give a segmenation fault.
(These two code snippets are inspired by the github discussion and the pyvista documentation).

Im am very grateful for any advice.