Time series data in pyvista - XdmfReader

I/O
2

Dear @nav, thanks for your post. I tried to run your code (in serial, i.e. by running python3 filename.py) and obtain the following error:

[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

Do you know what could be the issue?

I have been aware of this xdmf-reader feature, and tried to implement a workflow like you have posted above before, but failed to do so.
One example of what I tried to do is visible in this post. Your code is supposed to do exactly what I was asking for, so it would provide an answer to the question I posted there.

Looking forward to your answer.

EDIT: Here is a minimal example that reproduces the error:

import pyvista as pv
from dolfinx import fem, io, mesh 
from mpi4py import MPI

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = fem.FunctionSpace(domain, ("CG", 1))
phi = fem.Function(V)
phi.name = "phi"  # This determines the name in reader.set_active_scalars("phi")

# Create xdmf file for saving time series
xdmf = io.XDMFFile(MPI.COMM_WORLD, "phi.xdmf", "w")
xdmf.write_mesh(domain)

# Set phi to an arbitrary function at three times and save to xdmf file
t = 0.0
phi.interpolate(lambda x: (t + 0.1) * x[0])
xdmf.write_function(phi, t)
t = 0.5
phi.interpolate(lambda x: (t + 0.1) * x[0])
xdmf.write_function(phi, t)
t = 1.0
phi.interpolate(lambda x: (t + 0.1) * x[0])
xdmf.write_function(phi, t)
xdmf.close()

# read file and plot
xdmf_file = "phi.xdmf"
reader = pv.get_reader(xdmf_file)
reader.set_active_time_value(0.0)
# grid = reader.read()[1]
# grid.set_active_scalars("phi")
# p = pv.Plotter()
# p.add_mesh(grid, cmap="plasma", clim=[0.0,1.0])
# p.show()

The line reader.set_active_time_value(0.0) seems to cause the problem.