Case Study 26.2: Building a Time-Series Visualization Module for the Heat Solver

"A simulation you cannot see is a simulation you cannot debug, cannot present, and cannot trust." — folklore of computational science, and true

Executive Summary

Case Study 26.1 ported an inherited writer; here you build the real thing: a small, reusable visualization module, heat_viz, that turns the heat solver's time loop into a proper ParaView time series with physical time on the slider. The design centerpiece is a stateful series object that opens a .pvd collection, writes a VTK frame and appends a collection entry each time you save, and closes the collection cleanly at the end — so a single object manages the whole animation's bookkeeping. You will then add the Python side: a post-processing script that renders the frames into a movie and a publication figure. The result is the visualization half of the Chapter 38 capstone, built once and reused.

Skills applied

• Designing a stateful writer as a derived type with begin/add/end lifecycle (§26.2, §26.3; Chapter 9).
• Emitting a VTK time series with zero-padded names (§26.3; Chapter 12).
• Accumulating a .pvd collection so ParaView shows physical seconds (§26.3).
• Driving a matplotlib movie/figure pipeline from the saved frames (§26.5).
• Keeping I/O in its own module, decoupled from the solver (Chapter 8).

Background

The Project Checkpoint writes a VTK file per interval — enough to animate, but the frames carry only their index, not their physical time, and each run leaves an unmanaged pile of files. A production code wants three things the checkpoint lacks: physical time on the slider (a .pvd), a single object that owns the output session so the driver stays clean, and a clear seam between "compute" and "visualize." We design heat_viz to provide all three, and we hold the public interface small.

Phase 1 — Design the interface before the implementation

Good I/O modules expose a lifecycle, not a pile of subroutines. A time series has three moments — begin, add a frame, end — so the interface is three procedures around one state-carrying type:

type :: vtk_series_t
  integer                   :: pvd_unit = -1      ! open unit for the .pvd file
  integer                   :: nframes  = 0       ! how many frames written so far
  logical                   :: active   = .false. ! guards against misuse
end type vtk_series_t

! series_begin(series, pvd_name)          -- open the .pvd, write its header
! series_add(series, field, step, time)   -- write heat_<step>.vtk, append a .pvd entry
! series_end(series)                      -- write the .pvd footer, close it

The vtk_series_t object is the whole design: it hides the open .pvd unit and the frame count behind a three-call lifecycle. The driver never touches a filename or an XML tag — it says begin, add each save, end, and the module does the bookkeeping. This is the derived-type-as-owned-resource pattern from Chapter 9, applied to an output session.

Design rule. The driver's time loop should read like physics with one visualization call in it, not like a file-format tutorial. If the loop is cluttered with write statements and tag strings, the seam between compute and visualize is in the wrong place.

Phase 2 — Implement the module

The implementation reuses the canonical write_vtk (from §26.2) verbatim and adds the .pvd session around it. Here is the heart of heat_viz:

module heat_viz
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  private
  public :: field_t, vtk_series_t, series_begin, series_add, series_end

  type :: field_t
    integer  :: nx = 0, ny = 0
    real(dp) :: dx = 0.0_dp, dy = 0.0_dp
    real(dp), allocatable :: u(:,:)
  end type field_t

  type :: vtk_series_t
    integer :: pvd_unit = -1
    integer :: nframes  = 0
    logical :: active   = .false.
  end type vtk_series_t

contains

  subroutine series_begin(series, pvd_name)
    type(vtk_series_t), intent(out) :: series
    character(len=*),   intent(in)  :: pvd_name
    open(newunit=series%pvd_unit, file=pvd_name, status='replace', action='write')
    write(series%pvd_unit, '(a)') '<?xml version="1.0"?>'
    write(series%pvd_unit, '(a)') '<VTKFile type="Collection" version="1.0">'
    write(series%pvd_unit, '(a)') '  <Collection>'
    series%active  = .true.
    series%nframes = 0
  end subroutine series_begin

  subroutine series_add(series, field, step, time)
    type(vtk_series_t), intent(inout) :: series
    type(field_t),      intent(in)    :: field
    integer,            intent(in)    :: step
    real(dp),           intent(in)    :: time
    character(:), allocatable :: fname
    if (.not. series%active) error stop 'series_add: call series_begin first'
    fname = frame_name(step)
    call write_vtk(field, fname, step)                      ! the §26.2 writer
    write(series%pvd_unit, '(a, f0.4, a)') &
         '    <DataSet timestep="', time, '" file="' // fname // '"/>'
    series%nframes = series%nframes + 1
  end subroutine series_add

  subroutine series_end(series)
    type(vtk_series_t), intent(inout) :: series
    write(series%pvd_unit, '(a)') '  </Collection>'
    write(series%pvd_unit, '(a)') '</VTKFile>'
    close(series%pvd_unit)
    series%active = .false.
  end subroutine series_end

  function frame_name(step) result(name)          ! Chapter 12
    integer, intent(in)       :: step
    character(:), allocatable :: name
    character(len=32)         :: buf
    write(buf, '(a, i6.6, a)') 'heat_', step, '.vtk'
    name = trim(buf)
  end function frame_name

  ! write_vtk is the canonical §26.2 subroutine; omitted here for brevity.

end module heat_viz

The active flag turns a misuse (adding a frame before beginning the session) into a loud error stop rather than a write to a closed unit — the defensive habit of Chapter 13. Everything about the file format lives in this module; the solver knows none of it.

Phase 3 — Drive it, and read the collection it builds

The payoff is the driver's time loop, which now reads like physics with three visualization calls bracketing it:

type(vtk_series_t) :: series
real(dp), parameter :: dt = 0.0005_dp
integer, parameter  :: n_steps = 200, save_every = 100
integer :: step

call series_begin(series, 'heat.pvd')
do step = 0, n_steps
  if (mod(step, save_every) == 0) &
       call series_add(series, field, step, step * dt)
  call step_field(field, alpha, dt)
end do
call series_end(series)

Saving at steps 0, 100, 200 with dt = 0.0005 gives frame times 0.0, 0.05, 0.10. Trace the .pvd the session accumulates and you can predict it line for line:

<?xml version="1.0"?>
<VTKFile type="Collection" version="1.0">
  <Collection>
    <DataSet timestep="0.0000" file="heat_000000.vtk"/>
    <DataSet timestep="0.0500" file="heat_000100.vtk"/>
    <DataSet timestep="0.1000" file="heat_000200.vtk"/>
  </Collection>
</VTKFile>

Sanity check. Three <DataSet> lines for three saves; the timesteps are step*dt = 0*0.0005, 100*0.0005, 200*0.0005 = 0.0, 0.05, 0.10, formatted by f0.4 as 0.0000, 0.0500, 0.1000; and each file= matches the frame_name for that step. Open heat.pvd (not the individual .vtk files) in ParaView and the slider reads physical seconds, 0.00 to 0.10.

Phase 4 — The Python half: a movie from the frames

ParaView animates interactively, but for a slide you often want a self-contained MP4 or GIF. This is pure post-processing — the frames already exist — so it belongs in Python (§26.5). The script renders each saved matrix to a PNG and stitches them:

#!/usr/bin/env python3
"""make_movie.py -- render the solver's saved frames into an animated GIF."""
import glob
import numpy as np
import matplotlib.pyplot as plt
import imageio.v2 as imageio

frames = sorted(glob.glob("heat_*.dat"))     # zero-padded names sort in time order
vmin, vmax = 0.0, 100.0                        # fixed color scale across ALL frames
images = []
for fn in frames:
    field = np.loadtxt(fn)
    fig, ax = plt.subplots(figsize=(4, 4), constrained_layout=True)
    im = ax.imshow(field, cmap="inferno", vmin=vmin, vmax=vmax)
    fig.colorbar(im, ax=ax, label="T (K)")
    png = fn.replace(".dat", ".png")
    fig.savefig(png, dpi=120)
    plt.close(fig)
    images.append(imageio.imread(png))
imageio.mimsave("heat.gif", images, fps=10)
print(f"stitched {len(frames)} frames into heat.gif")

Design note — fix the color scale. The single most common mistake in simulation movies is letting each frame autoscale its colormap: the colors then mean different temperatures in every frame, and a cooling plate can look like it is not changing at all. Pin vmin/vmax to the run's global range (here 0–100 K), as above, so a color means the same temperature throughout the animation. It is the moving-picture cousin of the perceptually-uniform-colormap rule from §26.5: the visualization must not lie about magnitude.

Phase 5 — The publication figure

The final steady-state frame becomes the paper figure: a heat map with a contour overlay marking isotherms. The plot_heat.py of §26.5, extended with two lines, does it:

cs = ax.contour(field, levels=[20, 40, 60, 80], colors="white", linewidths=0.6)
ax.clabel(cs, inline=True, fontsize=7, fmt="%d K")

Reading the same saved data, this adds white isotherm contours labeled in kelvin over the inferno map, at 300 DPI. The figure and the animation are now two views of one dataset the solver wrote once — the mature form of Fortran and Python are better together: the field is computed at Fortran speed and rendered, in as many ways as the paper needs, by Python that never re-runs the physics.

Discussion Questions

1. The vtk_series_t object hides the open .pvd unit and the frame count. What bugs does bundling that state into one object prevent, compared with passing a bare unit number and an integer counter through the driver?
2. series_add guards with if (.not. series%active) error stop. Name two concrete misuse scenarios this catches, and relate it to the defensive-programming discipline of Chapter 13.
3. Phase 4 fixes the colormap range across all frames. Construct a scenario where autoscaling each frame would make a correct simulation look wrong to a viewer.
4. The whole module changes how the run is seen, not what it computes. Why is that a sign the seam between the heat_viz module and the solver modules is in the right place (Chapter 8)?

Your Turn: Extensions

• Option A (design). Give series_add an optional second scalar (say, the analytic error field) so each frame carries both temperature and error arrays; ParaView's dropdown then offers both. What changes in the VTK writer, and why can multiple SCALARS blocks share one POINT_DATA?
• Option B (optimize). For a large 3-D run, ASCII VTK is too big and too slow. Redesign series_add to write XML .vti with a binary DataArray (base64 or raw appended) instead of ASCII, and estimate the size reduction versus the ASCII series (tie to Chapter 25).
• Option C (build). Add a series_summary(series) that, at series_end, prints the frame count, the time span, and the total bytes written — a one-line provenance report for the run's log.

Key Takeaways

• Model an output session as an object with a lifecycle — begin/add/end around a state-carrying type keeps the driver clean and the bookkeeping in one place.
• A .pvd collection buys physical time on the slider for a handful of write statements — index time is rarely what you want to show.
• Fix the colormap range across an animation, or the colors lie about magnitude frame to frame.
• Compute in Fortran, render in Python — one saved dataset feeds an interactive ParaView animation, a stitched movie, and a contoured publication figure, none of which re-runs the simulation.