Chapter 26: Visualization Output — Writing Data That ParaView, VisIt, and Python Can Read

"The purpose of computing is insight, not numbers." — Richard W. Hamming, Numerical Methods for Scientists and Engineers

Overview

Your solver has been running for twenty-four chapters. Since Chapter 24 it is a real one: a finite-difference heat equation that marches a two-dimensional temperature field forward in time, CFL-safe timestep by CFL-safe timestep, until the plate settles. It is correct. It is even beginning to be fast. And it is completely, frustratingly blind. A real(dp) array of a hundred thousand numbers is the truth about how heat is spreading across your plate, and you cannot see a single thing in it. Print it and you get a wall of decimals. Somewhere in that wall a hot corner is cooling and a front is advancing, but your eyes will never find it.

This chapter gives the solver eyes. The job is narrow and entirely practical: take the field of numbers your simulation already computes and write it to disk in a format that a real scientific visualization tool can open — so that instead of reading 301.4 299.8 298.1 … you watch the heat diffuse, in color, with a slider you can drag through time. The dominant tools for this in computational science are ParaView and VisIt, both free, both built on the VTK library, and both able to render a billion-cell dataset on a laptop or a cluster. Feeding them is not hard — a VTK file is just a text file with a specific header — but it must be done exactly, because these readers are strict, and a header that is off by one keyword or one point count produces either an error dialog or, worse, a plausible picture of the wrong thing.

We will also take the lighter path. Not every plot needs ParaView. When you want a quick heat map or a line through the middle of the plate, a plain ASCII file and three lines of gnuplot, or a matplotlib script reading the array with numpy.loadtxt, will get you there in seconds — and this is where the book's sixth theme, Fortran and Python are better together, does real work: Fortran computes the field at full speed, and Python turns it into the figure that goes in the paper. By the end of the chapter your heat solver writes a VTK file every few steps, you open the whole sequence in ParaView as an animation, and you have a publication-quality figure produced by a post-processing script you can rerun any time the data changes.

In this chapter, you will learn to:

• Read the VTK data model well enough to know which of its formats your data wants, and write a valid legacy .vtk file for a structured field, keyword by keyword.
• Write the modern XML .vti (ImageData) format too, and understand the trade-off between the two.
• Emit one file per timestep, name the files so they sort in time order, and open the whole set as a time series in ParaView or VisIt — with a .pvd collection to attach real physical times.
• Write dead-simple gnuplot and matplotlib-ready ASCII, and know which format each tool wants.
• Turn a saved field into a publication-quality figure with a small, rerunnable matplotlib script.

Learning Paths

How to read this chapter by track. - 🔬 Scientist ("I need to see my results") — this is a core chapter for you. Read §26.2 (write the VTK file), §26.3 (open it in ParaView), and §26.5 (the matplotlib figure). The Project Checkpoint wires it straight into your solver. - 📖 Standard — read straight through; it closes Part VI and hands the capstone (Chapter 38) its figures. - 🔧 Legacy ("I inherited old code") — old scientific codes wrote bespoke ASCII dumps and drove gnuplot; §26.4 is your bridge, and §26.2 shows the small, standard format you should migrate them toward. - ⚡ HPC ("I run at scale") — VTK ASCII does not scale to a billion cells; read §26.2 for the format and §26.3 for the time-series mechanics, but pair them with the binary/parallel formats of Chapter 25. The naming and .pvd machinery here is exactly what a parallel run needs.

26.1 The VTK Format: A Data Model for Visualization

Before writing a single byte, you need a mental model of what a visualization file is. It is tempting to think of it as "a picture," but it is the opposite: a visualization file contains no picture at all. It contains the data — the geometry of your grid and the values on it — and the picture is made later, on demand, by the viewer. This separation is the whole point. You write the temperature field once; then in ParaView you can color it, slice it, contour it, warp it, animate it, and export a figure, all without your program running again. The file is the data; the visualization is a live question you ask of it.

The lingua franca for this data is VTK.

Definition (VTK). VTK is the Visualization Toolkit, an open-source C++ library for scientific visualization, and — the sense we care about — the family of file formats it defines for storing a dataset: a grid plus the field values living on it. ParaView and VisIt are both built on VTK and read its files natively. There are two generations of the format: the legacy format (a simple, human-readable .vtk text file, unchanged for decades) and the newer XML formats (.vti, .vtr, .vts, .vtu, one per grid type, supporting compression and parallel pieces). Both encode the same idea; the legacy format is the one to learn first because you can read it, write it, and debug it by eye.

The single most important choice VTK asks you to make is the type of grid, because it determines how much you must write. VTK's dataset types form a ladder from "the grid is implicit, write almost nothing" to "the grid is arbitrary, write everything":

  STRUCTURED_POINTS   regular box: give origin + spacing + dimensions.
  (a.k.a. ImageData)  Point positions are IMPLICIT — you write only the values.   ← our solver
  RECTILINEAR_GRID    axis-aligned but unequal spacing: give the coordinate arrays x(:), y(:), z(:).
  STRUCTURED_GRID     logically a box (i,j,k) but curved in space: give every point's (x,y,z).
  UNSTRUCTURED_GRID   arbitrary mesh: give every point AND the connectivity of every cell.

Your heat solver lives at the very top of that ladder, and that is a gift. The plate is a uniform grid: the temperature sits at points spaced dx apart in one direction and dy in the other, starting at an origin. Nothing about the geometry varies from run to run except three numbers. So the grid is implicit — VTK can reconstruct every point position from the origin, the spacing, and the dimensions — and the file needs to carry only the temperature values themselves. This is the STRUCTURED_POINTS dataset in the legacy format, called ImageData (extension .vti) in the XML format. It is the simplest thing VTK offers, and it is exactly what a finite-difference code on a regular grid produces.

Definition (structured grid output). Structured grid output writes field values on a logically rectangular grid — one where each point has an implied index (i, j, k) and the neighbors of a point are found simply by adding or subtracting one from an index. The grid's connectivity is therefore implicit and need not be stored; for a uniform grid (like our plate) even the point coordinates are implicit, recoverable from an origin and a spacing. This is the opposite of unstructured output, where points sit at arbitrary locations and every cell's list of corner points must be written out explicitly. Structured output is smaller, simpler, and faster to write — and it is what finite-difference simulations on regular grids naturally produce.

💡 Intuition: think of STRUCTURED_POINTS as a spreadsheet, and UNSTRUCTURED_GRID as a bag of surveyed GPS points with a separate list of which points form each triangle. The spreadsheet needs only its numbers, its row/column counts, and the physical size of a cell; everything else is implied by the grid. The bag of GPS points needs every coordinate and every triangle spelled out. Our plate is a spreadsheet, so we get to write the easy file.

A note on unstructured data. You will meet UNSTRUCTURED_GRID the moment your geometry stops being a box — a finite-element mesh of an engine part, a triangulated aircraft surface, an ocean model that follows a coastline. There you must write two things: a list of every point's (x, y, z) coordinates, and a list of cells, each naming the points that form it and a type code (triangle, tetrahedron, hexahedron). It is more work and a larger file, and VTK's unstructured format (legacy UNSTRUCTURED_GRID, XML .vtu) exists precisely to carry it. We will not write one here, because the heat solver never needs it; but know that the same tools, the same ParaView, read it, and that the leap from our structured file to an unstructured one is a matter of adding the points and cells, not learning a new tool. For a regular grid, resist the temptation — writing your uniform plate as an unstructured mesh would store hundreds of thousands of coordinates the grid already implies, for no gain.

🔗 Connection: VTK is the visualization counterpart to the storage formats of Chapter 25. NetCDF and HDF5 answer "how do I store a terabyte of field data portably and compressed?"; VTK answers "how do I hand a field to a viewer so a human can see it?" They are complementary, and large codes use both: HDF5 (or its VTK cousin) for the archival data, VTK files for the frames you actually look at. In fact ParaView reads NetCDF and HDF5 too, but through readers that expect specific conventions — the VTK file is the path of least resistance.

26.2 Writing Legacy VTK and XML .vti from Fortran

Now the payoff: writing the file. We will do the legacy .vtk format first, in full, because it is short enough to hold in your head and every byte of it is human-readable — which means when ParaView complains, you can open the file in a text editor and see exactly what is wrong.

The legacy VTK file, part by part

A legacy VTK file has a rigid five-part structure, and the reader parses it positionally, so order and spelling matter absolutely. Here are the five parts for a structured-points dataset:

1. # vtk DataFile Version 3.0     <- the FIRST line, verbatim; the magic string
2. <a one-line title/comment>      <- any text up to 256 characters; provenance goes here
3. ASCII                           <- data storage: ASCII or BINARY
4. DATASET STRUCTURED_POINTS       <- the grid type, then its parameters:
     DIMENSIONS  nx ny nz          <-   number of POINTS along each axis
     ORIGIN      x0 y0 z0          <-   physical position of point (1,1,1)
     SPACING     dx dy dz          <-   physical distance between adjacent points
5. POINT_DATA  (nx*ny*nz)          <- values live AT the points; the count MUST match
     SCALARS temperature double 1  <-   an array: name, type, components-per-point
     LOOKUP_TABLE default          <-   the color table (default = viewer's choice)
     <one value per point, in order>

Two details in that skeleton cause almost every VTK bug, so fix them in your mind now.

First: DIMENSIONS counts points, not cells, and POINT_DATA must equal their product. For an nx × ny field with nz = 1, there are nx*ny*1 points, and you must write exactly that many values after POINT_DATA. Write one too few or one too many and the reader either errors or silently misaligns every value — a picture of the wrong thing, which is the worst failure mode in all of visualization.

Second: the point ordering runs with the first axis fastest. VTK stores structured points with the x-index varying fastest, then y, then z. If we map our field's first index i to x and second index j to y, then the required output order is: for each j, for each i, write u(i,j). That is a loop with i on the inside — and, as we are about to see, that is not a coincidence you should ignore.

Here is the writer. It takes a field_t (the derived type from Chapter 9: nx, ny, dx, dy, and u(:,:)), the filename to write, and the timestep number (which it records in the title line for provenance). This is the canonical write_vtk your solver will call from now on:

module vtk_io
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  private
  public :: dp, field_t, write_vtk

  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

contains

  subroutine write_vtk(field, filename, step)
    type(field_t),    intent(in) :: field
    character(len=*), intent(in) :: filename
    integer,          intent(in) :: step
    integer :: iu, i, j, ios

    open(newunit=iu, file=filename, status='replace', action='write', iostat=ios)
    if (ios /= 0) error stop 'write_vtk: cannot open output file'

    ! --- the five mandatory header parts (order and spelling are exact) ---
    write(iu, '(a)')             '# vtk DataFile Version 3.0'
    write(iu, '(a, i0)')         'heat solver output, step ', step
    write(iu, '(a)')             'ASCII'
    write(iu, '(a)')             'DATASET STRUCTURED_POINTS'
    write(iu, '(a, 3(1x, i0))')  'DIMENSIONS', field%nx, field%ny, 1
    write(iu, '(a, 3(1x, f0.6))')'ORIGIN', 0.0_dp, 0.0_dp, 0.0_dp
    write(iu, '(a, 3(1x, f0.6))')'SPACING', field%dx, field%dy, 1.0_dp
    write(iu, '(a, 1x, i0)')     'POINT_DATA', field%nx * field%ny
    write(iu, '(a)')             'SCALARS temperature double 1'
    write(iu, '(a)')             'LOOKUP_TABLE default'

    ! --- the values: x (index i) fastest, matching column-major memory order ---
    do j = 1, field%ny
      do i = 1, field%nx
        write(iu, '(f0.6)') field%u(i, j)
      end do
    end do

    close(iu)
  end subroutine write_vtk

end module vtk_io

Drive it with a tiny 3×2 field whose values you can read at a glance, so the output file is short enough to verify against the format skeleton by eye:

program demo_vtk
  use vtk_io
  implicit none
  type(field_t) :: plate

  plate%nx = 3;  plate%ny = 2
  plate%dx = 0.5_dp;  plate%dy = 0.5_dp
  allocate(plate%u(3, 2))
  plate%u(:, 1) = [0.0_dp, 1.0_dp, 2.0_dp]   ! bottom row  (j = 1)
  plate%u(:, 2) = [3.0_dp, 4.0_dp, 5.0_dp]   ! top row     (j = 2)

  call write_vtk(plate, 'heat_000100.vtk', 100)
  print '(a)', 'wrote heat_000100.vtk'
end program demo_vtk

$ gfortran -std=f2018 -Wall example-01-legacy-vtk.f90 -o demo_vtk && ./demo_vtk
wrote heat_000100.vtk

And here is the exact file it writes — the object of the whole exercise. Compile the program, run it, and open heat_000100.vtk; it should match this character for character:

# vtk DataFile Version 3.0
heat solver output, step 100
ASCII
DATASET STRUCTURED_POINTS
DIMENSIONS 3 2 1
ORIGIN 0.000000 0.000000 0.000000
SPACING 0.500000 0.500000 1.000000
POINT_DATA 6
SCALARS temperature double 1
LOOKUP_TABLE default
0.000000
1.000000
2.000000
3.000000
4.000000
5.000000

Read it against the skeleton. Line 1 is the magic string, verbatim — get one character wrong and the file is not a VTK file. Line 2 is our provenance comment, carrying the step number so a stray file on disk can identify itself. ASCII says the values are text. DATASET STRUCTURED_POINTS with DIMENSIONS 3 2 1 declares a 3-by-2-by-1 grid of points; ORIGIN places point (1,1) at the physical origin; SPACING says adjacent points are 0.5 apart in x and y (the third value, 1.0, is a harmless placeholder for the unused z axis). POINT_DATA 6 promises exactly 3 × 2 × 1 = 6 values, and six values follow. SCALARS temperature double 1 names the array temperature, declares it double precision with one component per point, and LOOKUP_TABLE default hands the coloring to the viewer. Then the six values, in VTK order: u(1,1), u(2,1), u(3,1) — the bottom row — then u(1,2), u(2,2), u(3,2) — the top row — which are exactly 0, 1, 2, 3, 4, 5. That file, tiny as it is, is a complete, valid dataset ParaView will open and shade.

⚠️ Common Pitfall — the POINT_DATA count must match, exactly. The single most common VTK error is a POINT_DATA value that disagrees with the number of values written, usually because someone counted cells instead of points, or looped nx*ny but wrote the header with a stale dimension. If the count is too large, the reader hits end-of-file mid-array and errors; if it is too small, it reads a valid-looking but shifted field and shows you a lie. Always compute the count as field%nx * field%ny from the same variables the loop uses, never as a hard-coded number, so the header and the body cannot drift apart.

⚡ Performance Note — the write order is column-major, and that is on purpose. VTK wants the x-index fastest; Fortran stores u(i,j) with the first index fastest (column-major order, the central idea of Chapter 5). Because we mapped i to x, the required loop — do j; do i; write u(i,j) — walks the array u(1,1), u(2,1), … in exactly the order it sits in memory, so the write streams contiguously through the cache instead of jumping by a row each step. Map your fastest array index to VTK's fastest axis and the file writer is automatically cache-friendly; map them the other way and every write is a stride away from the last. It is a small instance of a large truth in this book — arrays are Fortran's superpower, and performance is not accidental — showing up even in something as mundane as writing a file. (For a plate of a hundred thousand points this hardly matters; for a three-dimensional field of a hundred million, written every few steps, it is the difference between I/O you notice and I/O you do not.)

The modern XML format: .vti

The legacy format is perfect for learning and debugging, but production codes increasingly prefer VTK's XML formats, which support in-file compression, streaming, and parallel pieces (each MPI rank writes its own chunk, and a small master file stitches them together — a natural fit for the distributed solver of Chapter 34). For our uniform grid the XML type is ImageData, written to a .vti file.

Definition (VTI). A .vti file is VTK's XML ImageData format — the XML-era equivalent of the legacy STRUCTURED_POINTS dataset. It stores the same thing (a uniform grid given by whole-extent, origin, and spacing, plus point data) but as well-formed XML, which lets it carry compressed binary payloads, multiple named arrays, and (via a companion collection file) a whole time series. The XML is more verbose than the legacy header, but it is what modern pipelines and parallel writers emit.

The structure mirrors the legacy file, wrapped in XML tags. The one new concept is extent, given as six integers i0 i1 j0 j1 k0 k1 — the inclusive index range of points along each axis. For our 3×2×1 grid the points are indexed 0..2 in x, 0..1 in y, and 0..0 in z, so the whole extent is 0 2 0 1 0 0. (Note VTK's XML uses 0-based point indices in the extent, even though your Fortran array is 1-based — a small translation the writer makes.) Here is a compact .vti writer for the same field:

module vti_io
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use vtk_io, only: field_t
  implicit none
  private
  public :: write_vti

contains

  subroutine write_vti(field, filename, step)
    type(field_t),    intent(in) :: field
    character(len=*), intent(in) :: filename
    integer,          intent(in) :: step
    integer :: iu, i, j
    character(len=64) :: extent, spacing

    open(newunit=iu, file=filename, status='replace', action='write')
    ! build the two attribute strings we need (0-based inclusive extent; spacing) 
    write(extent,  '(5(i0,1x),i0)')        0, field%nx-1, 0, field%ny-1, 0, 0
    write(spacing, '(f0.6,1x,f0.6,1x,f0.6)') field%dx, field%dy, 1.0_dp

    write(iu, '(a)') '<?xml version="1.0"?>'
    write(iu, '(a,i0,a)') '<VTKFile type="ImageData" version="1.0" '// &
                          'byte_order="LittleEndian"> <!-- step ', step, ' -->'
    write(iu, '(a)') '  <ImageData WholeExtent="'//trim(extent)// &
                     '" Origin="0 0 0" Spacing="'//trim(spacing)//'">'
    write(iu, '(a)') '    <Piece Extent="'//trim(extent)//'">'
    write(iu, '(a)') '      <PointData Scalars="temperature">'
    write(iu, '(a)') '        <DataArray type="Float64" Name="temperature" format="ascii">'
    do j = 1, field%ny
      do i = 1, field%nx
        write(iu, '(f0.6)') field%u(i, j)
      end do
    end do
    write(iu, '(a)') '        </DataArray>'
    write(iu, '(a)') '      </PointData>'
    write(iu, '(a)') '    </Piece>'
    write(iu, '(a)') '  </ImageData>'
    write(iu, '(a)') '</VTKFile>'
    close(iu)
  end subroutine write_vti

end module vti_io

For the same 3×2 field, write_vti(plate, 'heat_000100.vti', 100) produces:

<?xml version="1.0"?>
<VTKFile type="ImageData" version="1.0" byte_order="LittleEndian"> <!-- step 100 -->
  <ImageData WholeExtent="0 2 0 1 0 0" Origin="0 0 0" Spacing="0.500000 0.500000 1.000000">
    <Piece Extent="0 2 0 1 0 0">
      <PointData Scalars="temperature">
        <DataArray type="Float64" Name="temperature" format="ascii">
0.000000
1.000000
2.000000
3.000000
4.000000
5.000000
        </DataArray>
      </PointData>
    </Piece>
  </ImageData>
</VTKFile>

The two files carry identical information — same grid, same six values, same ordering — in two dressings. Which should you write? For a first solver, write the legacy .vtk: it is easier to generate, trivial to inspect, and ParaView reads it instantly. Reach for .vti when you need what XML buys — compression for large fields, or the per-rank parallel pieces (.pvti) that a distributed run wants. The book's default, and the Project Checkpoint's, is legacy VTK; the .vti writer is here so that when a colleague's pipeline demands XML, you already have it.

🐍 Python Comparison: In a pure-Python workflow you would not hand-write either format — you would call a library: pyvista or meshio (both real, widely used packages) build a grid object and write .vti/.vtk for you, and evtk writes VTK from NumPy arrays directly. Those libraries are excellent, and if Python owns your I/O, use them. But your solver is Fortran, running in a tight time loop where calling back into Python every few steps would be absurd — so the field gets written where it lives, in Fortran, with the twenty lines above. This is the division of labor the whole chapter rests on: Fortran and Python are better together — Fortran writes the data at speed, and Python (or ParaView) reads it back to make the picture.

🔄 Check Your Understanding. 1. A field is nx = 100, ny = 80. What number must follow POINT_DATA, and how many value lines follow the LOOKUP_TABLE default line? 2. Why does the value-writing loop put i on the inside rather than j? 3. In the .vti file, why is the whole extent 0 99 0 79 0 0 and not 1 100 1 80 1 1?

Answers

1. POINT_DATA must be 100 × 80 × 1 = 8000, and exactly 8000 value lines follow (one scalar per point). The header count and the number of values must agree. 2. VTK stores structured points with the x-index (our i) varying fastest, so the inner loop must be over i. As a bonus, since Fortran is column-major, that inner i loop also walks memory contiguously. 3. VTK's XML extent uses 0-based inclusive point indices: 100 points along x are indexed 0 through 99. Fortran's array is 1-based (1..100), so the writer subtracts one when it emits the extent.

26.3 Opening the Output in ParaView and VisIt — and a Time Series

A file on disk is not yet a picture. This section is the bridge from your .vtk file to the moving, colored plate on your screen, in the two tools that dominate scientific visualization.

ParaView and VisIt are both free, open-source, cross-platform, and built on VTK; both scale from a laptop to the largest supercomputers by running a rendering server alongside your data. They differ in interface and heritage — ParaView (Kitware) and VisIt (Lawrence Livermore) — but for our purposes they are interchangeable, and everything below has a direct equivalent in each. We will describe ParaView, because it is the more common starting point.

Opening a single file is three clicks. File ▸ Open, pick heat_000100.vtk, and press the green Apply button in the Properties panel — nothing renders until you Apply, a surprise for every newcomer. You will see a flat gray rectangle: the grid, uncolored. In the coloring dropdown on the toolbar, change Solid Color to temperature, and the plate lights up — blue where it is cool, red where it is hot, with a color legend you can edit. That is the entire loop: write, open, Apply, color. Everything else ParaView does — contours, slices, warping the surface by temperature into a 3-D relief, computing the gradient — is a filter you add on top of that colored dataset.

But a single frame is not why we are here. We want to watch heat move, and that means a time series.

Definition (time-series output). Time-series output is a sequence of files, one per saved timestep, that a viewer loads together and plays as an animation. Each file is an ordinary snapshot of the field at one instant; what turns the set into a movie is (a) a shared base name with a zero-padded, increasing index — heat_000000.vtk, heat_000100.vtk, heat_000200.vtk, … — that lets the viewer detect the group and order it, and (b) optionally a small collection file that names each snapshot and its physical time. The zero-padding is not cosmetic: it makes the filenames sort in numeric order under plain alphabetical sorting, which is how ParaView orders the frames.

This is where the humble filename helper from Chapter 12 earns its keep. Recall frame_name(step), which turns a step number into a zero-padded name:

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

The i6.6 descriptor zero-pads the step to six digits, so heat_000100.vtk sorts before heat_001000.vtk under the ordinary alphabetical sort ParaView uses — and that is the entire reason we padded. When you open a folder of these files, ParaView collapses the whole numbered set into one entry named heat_..vtk (note the two dots — its wildcard for the varying digits). Open that, Apply, color by temperature, and the toolbar's play button ▶ appears with a time slider. Drag it and you scrub through your simulation; press play and you watch the plate evolve. The picture you have been unable to see for twenty-four chapters is finally moving.

   heat-solver/output/                ParaView sees one time-varying source:
   ├── heat_000000.vtk   ┐
   ├── heat_000100.vtk   │  detected as
   ├── heat_000200.vtk   ├───────────►   heat_..vtk   ◀── open this
   ├── heat_000300.vtk   │               [ |◀  ◀  ▶  ▶| ]  ← time controls
   └── …                 ┘               time slider: ●──────────

VisIt does the same thing through a slightly different door: it groups a numbered sequence into a database automatically, and its Time controls step through the frames. If you point VisIt at heat_*.vtk it offers to open the group as a time-varying database; from there, add a Pseudocolor plot of temperature and use the animation controls. The concept is identical; only the menu names change.

Attaching real physical times with a .pvd file. By default ParaView labels the frames 0, 1, 2, … — the file index, not the physical time. Your step 100 might be t = 0.05 seconds of simulated diffusion, and you would rather the slider said so. The clean way is a PVD collection file: a tiny XML file that lists each snapshot with its physical time.

Definition (PVD). A .pvd file is a small ParaView collection file: XML that maps each data file to a physical timestep value. You open the single .pvd, and ParaView loads the whole series with correct times on the slider — decoupling the physical time from the file index, so you can save every hundredth step yet still see true seconds.

You write it once, after the run, listing the frames you saved:

<?xml version="1.0"?>
<VTKFile type="Collection" version="1.0">
  <Collection>
    <DataSet timestep="0.00" file="heat_000000.vtk"/>
    <DataSet timestep="0.05" file="heat_000100.vtk"/>
    <DataSet timestep="0.10" file="heat_000200.vtk"/>
    <DataSet timestep="0.15" file="heat_000300.vtk"/>
  </Collection>
</VTKFile>

Emitting this from Fortran is a handful of write statements accumulated in the time loop (open the .pvd, write the header, append one <DataSet .../> line each time you save a frame, close with the two end tags); the Project Checkpoint sketches it. Open heat.pvd in ParaView and the slider now reads 0.00, 0.05, 0.10, 0.15 — physical seconds, not frame numbers. The capstone in Chapter 38 uses exactly this to present its results with correct time axes.

🧩 Try It Yourself: Before reading on, write three tiny fields by hand — a plate that is all zeros, one with a hot spot in the center, and one where the hot spot has spread — as frame_000000.vtk, frame_000001.vtk, frame_000002.vtk (3×3 grids are plenty). Open the set in ParaView, color by temperature, and press play. Watching your own three frames animate teaches the pipeline better than any screenshot could, and it takes five minutes.

26.4 Simple Formats for gnuplot and Python

VTK and ParaView are the right tools for a field evolving in time, but they are heavy artillery for a quick look. Half the time what you want is a fast heat map of one snapshot, or a line plot of temperature along the plate's centerline, and for that a plain ASCII file plus gnuplot or a few lines of matplotlib is faster and lighter. The good news is that these tools read almost anything, so the Fortran side is trivial — you already learned every tool you need in Chapter 7.

A matrix for matplotlib. The simplest useful format is the field written as a rectangular block of numbers: ny lines, each with nx values separated by spaces. NumPy's loadtxt reads that straight into a 2-D array with no arguments, and imshow turns the array into a heat map. Writing it is one nested loop:

subroutine write_matrix(field, filename)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use vtk_io, only: field_t
  implicit none
  type(field_t),    intent(in) :: field
  character(len=*), intent(in) :: filename
  integer :: iu, i, j
  open(newunit=iu, file=filename, status='replace', action='write')
  do j = field%ny, 1, -1                       ! top row first, so the file reads
    write(iu, '(*(f8.3, 1x))') (field%u(i, j), i = 1, field%nx)  ! like the plate looks
  end do
  close(iu)
end subroutine write_matrix

The '(*(f8.3, 1x))' format uses the unlimited repeat count *(...), which applies the parenthesized group to as many values as the output list supplies — here the whole row (field%u(i,j), i=1,field%nx) — so one write emits a full line regardless of nx. We loop j from ny down to 1 so the first line of the file is the top row of the plate; then imshow's default orientation shows the plate the right way up without fiddling with origin. For our 3×2 field the file is just:

   3.000    4.000    5.000
   0.000    1.000    2.000

Two lines, three values each: the top row (3 4 5) then the bottom row (0 1 2), exactly the plate seen from above. numpy.loadtxt('heat.dat') reads this as a shape-(2, 3) array, ready to plot.

The three-column format for gnuplot. gnuplot's surface and heat-map modes (splot … with pm3d) like a different shape: three columns x y value, with a blank line separating each scan line (each row of constant y). The blank line is how gnuplot knows the grid's row structure. Writing it adds the coordinates:

subroutine write_gnuplot(field, filename)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use vtk_io, only: field_t
  implicit none
  type(field_t),    intent(in) :: field
  character(len=*), intent(in) :: filename
  integer :: iu, i, j
  open(newunit=iu, file=filename, status='replace', action='write')
  do j = 1, field%ny
    do i = 1, field%nx
      write(iu, '(f0.3, 1x, f0.3, 1x, f0.6)') &
           (i-1)*field%dx, (j-1)*field%dy, field%u(i, j)
    end do
    write(iu, '(a)') ''                       ! blank line ends this scan for gnuplot
  end do
  close(iu)
end subroutine write_gnuplot

For the 3×2 field this writes each point's physical coordinate and value, with a blank line after each row:

0.000 0.000 0.000000
0.500 0.000 1.000000
1.000 0.000 2.000000

0.000 0.500 3.000000
0.500 0.500 4.000000
1.000 0.500 5.000000

Then three lines at the gnuplot prompt draw the heat map:

$ gnuplot
gnuplot> set pm3d map
gnuplot> set palette
gnuplot> splot 'heat.dat' using 1:2:3 with pm3d

⚠️ Common Pitfall — gnuplot's blank line is structural, not decorative. In the three-column grid format, the blank line after each scan (each row of constant y) tells gnuplot where one grid line ends and the next begins. Omit it and gnuplot reads your data as one long unstructured list and either draws nonsense or refuses to make a surface. Emit exactly one blank line between scans and none is missing — a doubled blank line means "the grid is broken here" and also misbehaves. It is the gnuplot analog of the VTK POINT_DATA count: a structural marker the reader depends on.

📜 From History: for decades gnuplot was how scientific Fortran saw its results. A code would write a column of numbers, a gnuplot script would plot them, and that was the visualization pipeline — simple, scriptable, and still perfectly good for a line plot or a quick surface. It is a fine example of the theme that Fortran is not dead: the language never needed a built-in graphics library because the Unix philosophy of small cooperating tools gave it gnuplot, and later ParaView and matplotlib, each reading a plain file. Your solver writes numbers; the visualization is somebody else's specialized job. That separation has aged extremely well.

🔄 Check Your Understanding. 1. Why does write_matrix loop j from ny down to 1 instead of 1 up to ny? 2. What does the * in the format '(*(f8.3, 1x))' do, and why is it convenient here? 3. In the gnuplot format, what breaks if you forget the blank line between scans?

Answers

1. So the file's first line is the plate's top row (largest y). matplotlib's imshow draws the first row at the top by default, so writing top-row-first makes the image match the physical plate without setting origin='lower'. 2. The * is the unlimited repeat count: it applies the group (f8.3,1x) to however many values the I/O list provides, so a single format writes a full row of any length nx without hard-coding a repeat count. 3. gnuplot loses the grid's row structure and can no longer build a surface — it treats the points as one undifferentiated list and draws nonsense (or refuses). The blank line is the row separator it requires.

26.5 Producing Publication-Quality Figures with matplotlib

The final step is the figure that goes in the paper, the thesis, or the slide — and here Python and matplotlib are simply the best tool for the job, reading the ASCII your Fortran solver wrote. This is the mature form of the sixth theme: the fast, correct field comes from Fortran; the polished, labeled, color-mapped figure comes from a post-processing script you can rerun forever. Separating the two means you can refine the figure — a better colormap, a contour overlay, a bigger font for the slide — without ever touching or rerunning the simulation.

Here is a complete, rerunnable post-processing script. It reads the matrix file from §26.4, makes a heat map with a labeled colorbar, and saves it at 300 DPI — the resolution journals ask for:

#!/usr/bin/env python3
"""plot_heat.py -- publication-quality heat map from the solver's matrix output.
   Usage:  python plot_heat.py heat.dat heat.png
   Reads the ASCII matrix written by write_matrix (top row first)."""
import sys
import numpy as np
import matplotlib.pyplot as plt

infile  = sys.argv[1] if len(sys.argv) > 1 else 'heat.dat'
outfile = sys.argv[2] if len(sys.argv) > 2 else 'heat.png'

field = np.loadtxt(infile)                 # shape (ny, nx); top row first
print(f'loaded {field.shape[0]} x {field.shape[1]} grid, '
      f'min={field.min():.3f}, max={field.max():.3f}')

fig, ax = plt.subplots(figsize=(5, 4), constrained_layout=True)
im = ax.imshow(field, cmap='inferno', aspect='equal')   # 'inferno' is perceptually uniform
cbar = fig.colorbar(im, ax=ax)
cbar.set_label('Temperature (K)')
ax.set_title('Heat diffusion on a square plate')
ax.set_xlabel('x node')
ax.set_ylabel('y node')
fig.savefig(outfile, dpi=300)              # 300 DPI: journal-ready raster
print(f'wrote {outfile}')

Run against the tiny 3×2 file from §26.4 and it reports what it read before writing the image, which is the part we can verify by hand:

$ python plot_heat.py heat.dat heat.png
loaded 2 x 3 grid, min=0.000, max=5.000
wrote heat.png

The two numbers are exactly right for our example field — its values run from 0.000 to 5.000, laid out as a 2-row, 3-column grid — which is the check that the file round-tripped correctly from Fortran to NumPy. (The image itself is not text, so this printed line is the honest "expected output"; the heat.png it writes is a small labeled heat map.) One choice in that script deserves a callout, because it is where visualization crosses from cosmetic to ethical.

Definition (colormap). A colormap (or color map) is the function that turns each scalar value into a color for display. The choice is not merely aesthetic: a perceptually uniform colormap — where equal steps in value produce equal-looking steps in color — represents the data honestly, while the old rainbow/"jet" map (still many tools' default) creates false boundaries at its yellow and cyan bands and hides detail in its greens, literally showing structure that is not in the data. matplotlib's viridis, inferno, magma, and plasma are perceptually uniform and safe for readers with color-vision deficiency; reach for one of them and avoid jet.

⚠️ Common Pitfall — the rainbow colormap lies. It is worth stating plainly: a jet/rainbow colormap can invent features. Its sharp yellow-green and cyan transitions read as edges in your data even when the field is perfectly smooth, and its dark-red-to-dark-blue endpoints compress the extremes. Reviewers in many fields now flag rainbow figures for exactly this reason. Defaulting to viridis or inferno — as the script above does — is a one-word change that makes every figure both more honest and more accessible.

🐍 Python Comparison: notice what this script is not: it is not doing any physics. It loads an array and draws it. All the computation — the finite-difference stepping, the CFL check, the boundary conditions — happened in Fortran, at Fortran speed, and left behind a file. matplotlib is doing what it is world-class at (typography, color, layout, export formats) and nothing else. That is the collaboration in its ideal form: had you tried to run the simulation in pure Python it might be fifty times slower (Chapter 15 measures exactly that), and had you tried to make this figure in Fortran you would have written a graphics library from scratch. Each language does the half it is best at. Fortran and Python are better together.

For a time series, the same script becomes a loop: read frame_000000.dat, frame_000100.dat, …, render each to a PNG, and stitch the PNGs into an animated GIF or MP4 with imageio or ffmpeg. That is a common way to put a simulation movie in a talk without carrying ParaView to the podium — and it is a natural extension left to the exercises.

Project Checkpoint

Your solver can finally be seen. This checkpoint gives the heat solver its visualization output: it calls the canonical write_vtk (from §26.2) once every save_every steps, using frame_name (from Chapter 12) to name the frames, producing the time series you open in ParaView.

The write_vtk subroutine belongs in the heat_io module, with the canonical signature you saw in full in §26.2:

subroutine write_vtk(field, filename, step)   ! canonical signature — see §26.2 for the body
  type(field_t),    intent(in) :: field
  character(len=*), intent(in) :: filename
  integer,          intent(in) :: step

The new work is the output cadence in the driver's time loop. You do not want a file every step — that is gigabytes of frames no one will scrub through — so you save every save_every steps, naming each frame by its step number:

integer, parameter :: save_every = 100
integer :: step

do step = 0, n_steps
  if (mod(step, save_every) == 0) then
    call write_vtk(field, frame_name(step), step)   ! frame_name -> 'heat_000100.vtk'
  end if
  call step_field(field, alpha, dt)                  ! the Ch.24 finite-difference update
end do

Trace the naming for the first few saves and you can predict the folder's contents exactly:

   step 0      -> write_vtk(field, 'heat_000000.vtk', 0)
   step 100    -> write_vtk(field, 'heat_000100.vtk', 100)
   step 200    -> write_vtk(field, 'heat_000200.vtk', 200)
   step 300    -> write_vtk(field, 'heat_000300.vtk', 300)

Those zero-padded names sort in time order, so ParaView collapses them into one heat_..vtk source with a play button — exactly the payoff §26.3 described. Open it, color by temperature, press play, and watch the hot edge bleed into the cold interior until the plate reaches steady state. For the first time, the number in the array and the picture in your head are the same thing.

Record three items in your heat-solver/ notes. First, the output cadence: save_every = 100 is a starting point — tune it so a run produces a few hundred frames, not a few hundred thousand. Second, put the frames in their own subfolder (output/) so heat_*.vtk does not clutter your source tree, and add output/ to your .gitignore — generated data does not belong in version control. Third, note the optional next step for the capstone: accumulate a heat.pvd collection (§26.3) as you save, so the Chapter 38 write-up can show a slider in physical seconds. Your solver now computes a field and shows it — the last capability it needed before Part VII makes it fast and Part VIII makes it parallel.

Summary

This chapter gave the solver its eyes: it writes field data in formats that ParaView, VisIt, gnuplot, and matplotlib read, and you learned to open the result as a picture and an animation.

The three things to memorize. First, a legacy VTK structured-points file is five fixed parts, and POINT_DATA must equal nx*ny*nz exactly — the count is the bug that bites everyone. Second, write the values with i innermost (do j; do i; write u(i,j)): it satisfies VTK's x-fastest ordering and Fortran's column-major layout at once. Third, let each tool do its half — Fortran writes the field, ParaView animates it, matplotlib makes the paper figure; that division is the sixth theme, Fortran and Python are better together, at work.

Spaced Review

Retrieval practice on the two chapters this one stands on: I/O (Chapter 7), whose write and edit descriptors produce every byte of the VTK file, and the PDE solver (Chapter 24), whose field is the thing we visualize. Answer before opening the details.

1. (Ch. 7) In write(iu, '(a, i0)') 'POINT_DATA ', 6000, what exact text lands in the record, and why is i0 the right descriptor for a point count?
   

   Answer

   The record is POINT_DATA 6000. The a descriptor writes the string 'POINT_DATA ' (with its trailing space) exactly, and i0 writes the integer in its minimum width — no leading blanks, and it never overflows into asterisks the way a fixed width like i4 would for a large grid. A point count can be anything from a few to hundreds of millions, so a self-sizing i0 is exactly right.

2. (Ch. 7) Why do we open with status='replace' for each output frame rather than status='new'?
   

   Answer

   status='replace' creates the file if it does not exist and overwrites it if it does, so re-running the solver cleanly regenerates its frames. status='new' requires the file not to exist and errors if it does — which would make the second run of a simulation crash on its first output file. For regenerable output, replace is the correct choice.

3. (Ch. 24) The VTK file stores the temperature at grid points, which is why we use POINT_DATA. In the five-point finite-difference stencil, where do the values u(i,j) conceptually live — at points or in cells — and why does that make POINT_DATA (not CELL_DATA) correct?
   

   Answer

   In a finite-difference scheme the unknowns u(i,j) live at the grid points (nodes): the stencil approximates the Laplacian at a point from that point and its four neighbors. Because the values are nodal, they attach to VTK's points, so POINT_DATA is correct. CELL_DATA would be right for a finite-volume scheme that stores one value per cell.

4. (Ch. 24) Your solver's field is nx = 200, ny = 200, and you write a .vtk file every 100 steps for 10,000 steps. How many frames result, and roughly how many scalar values does each frame contain?
   

   Answer

   Saving every 100 steps over steps 0..10000 gives frames at 0, 100, …, 10000 — that is 101 frames (don't forget the frame at step 0). Each frame is a full field of 200 × 200 = 40,000 scalar values, following its POINT_DATA 40000 header.

5. (Ch. 24) Why is a zero-padded step number in the filename (heat_000100.vtk, not heat_100.vtk) essential for a time series, and what would go wrong without it?
   

   Answer

   ParaView orders frames by alphabetical (lexical) sort of the filenames. Zero-padding to a fixed width makes lexical order match numeric order, so the animation plays in time order. Without padding, heat_100.vtk sorts before heat_1000.vtk but also before heat_99.vtk (because '1' < '9'), scrambling the sequence — the plate would appear to jump around in time instead of diffusing smoothly. The i6.6 in frame_name is what prevents this.

What's Next

Your solver is now complete in capability: it reads a configuration, marches a validated finite-difference scheme, and writes output you can watch diffuse in ParaView. What it is not yet is fast — and for the computations Fortran exists to run, fast is the whole point. Part VII opens with Chapter 27, which finally makes rigorous the performance claims this book has been making since Chapter 1: why a change of loop order can make the stencil ten times faster, why Fortran's no-aliasing rule lets the compiler optimize where C cannot, and how the column-major layout you just exploited to write a VTK file efficiently is the same idea that governs how fast your solver runs. You have built a simulation that works and that you can see. Now we make it quick.