Time series data in pyvista - XdmfReader

I/O
1

Version 0.39.0 of pyvista, released three days ago, now features an XdmfReader.

I’ve been using it quite a lot for post-processing time-series data created with fenicsx, so I thought I’d provide a simple example in case anyone else was interested.

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

comm = MPI.COMM_WORLD
rank = comm.rank 

ncells = 100
domain = mesh.create_unit_square(comm, ncells, ncells)

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(comm, "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()

# Plot with pyvista on one rank - most likely you would put the following in a 
# separate script and only run it with one core
if rank == 0:
    def plot_sol(func, min, max):
        """ Analyse 2d solution over time """
        cmap = "plasma"
        def load_time(value):
            """ Load solution at value specified using the slider """
            reader.set_active_time_value(value)
            grid = reader.read()[1]
            grid.set_active_scalars(func)
            p.add_mesh(grid, cmap=cmap, clim=[min, max], scalar_bar_args=sargs)

        # Load in data
        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(func)
        # Create plot
        width = 0.6
        sargs = dict(width=width, height=0.1, position_x=(1-width)/2, position_y=0.01)
        p = pv.Plotter()
        p.add_mesh(grid, cmap=cmap, clim=[min, max], scalar_bar_args=sargs)
        p.show_bounds(xtitle="", ytitle="", ztitle="")
        p.view_xy()
        p.add_slider_widget(load_time, [0.0, 1.0], value=0.0, title="Time",
                            interaction_event="always", pointa=(0.25, 0.93), 
                            pointb=(0.75, 0.93))
        p.show()

    plot_sol("phi", 0.0, 1.1)

Some example screenshots:

p1
p11024×768 11.1 KB

p2
p21024×768 11.9 KB

p3
p31024×768 11.2 KB

3 Likes
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.

3

Hi @renol-kth,

Just to check, does the xdmf file save with no errors when you run the script with python3?

Also, what version of pyvista are you using?