Case Study 26.1: The Inherited Diffusion Code — From a Home-Grown Dump to ParaView

"The first thing you do with an inherited code is not run it. It is read what it writes." — a maxim from more than one national-lab onboarding

Executive Summary

You have inherited a small, working 2-D diffusion code, diffuse.f90, from a departed colleague. It runs, it is numerically fine, and it produces a file called field.out in a format the colleague invented, which a brittle plot.gnu script turns into a single static image. Your advisor wants an animation of the run for a group meeting on Friday, and ParaView cannot read field.out at all. This case study is the analysis and the fix: you will read the home-grown format, verify it against a field you know, diagnose exactly why it neither animates nor opens in ParaView, and port the output to standard legacy VTK so the whole run becomes a scrubbable time series — without touching one line of the (correct) physics.

Skills applied

• Reading an undocumented ASCII output format and recovering its layout (§26.1, §26.2).
• Verifying a writer against a hand-known field before trusting it (§26.2).
• Diagnosing why a single-file, custom dump cannot become a time series (§26.3).
• Porting output to a standard VTK STRUCTURED_POINTS writer with the exact header (§26.2).
• Loading the resulting frame sequence as an animation in ParaView (§26.3).

Background

The inherited writer looks like this — modern enough to compile, but home-grown in format:

subroutine dump(u, nx, ny, fname)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  integer,          intent(in) :: nx, ny
  real(dp),         intent(in) :: u(nx, ny)
  character(len=*), intent(in) :: fname
  integer :: iu, i, j
  open(newunit=iu, file=fname, status='replace', action='write')
  write(iu, '(2i6)') nx, ny                          ! a header of two integers
  do i = 1, nx
    write(iu, '(*(es13.5))') (u(i, j), j = 1, ny)    ! one line per i, all j across
  end do
  close(iu)
end subroutine dump

It writes a two-integer header, then one line per i with all the j values across the row. It is not wrong — it round-trips fine into the colleague's gnuplot script — but it is nobody's standard, and no visualization tool knows how to read it. Your job is to understand it well enough to replace it safely.

Phase 1 — Read the format against a field you know

Never reverse-engineer a format from a big run; use the smallest field whose every value you chose. Drive the inherited dump with a 3×2 field:

real(dp) :: u(3, 2)
u(:, 1) = [0.0_dp, 1.0_dp, 2.0_dp]
u(:, 2) = [3.0_dp, 4.0_dp, 5.0_dp]
call dump(u, 3, 2, 'field.out')

Trace the writer by hand. The header line is nx ny in 2i6, so 3 2. Then the loop runs over i (the outer loop) and writes all j on each line in es13.5:

     3     2
  0.00000E+00  3.00000E+00
  1.00000E+00  4.00000E+00
  2.00000E+00  5.00000E+00

Now you know the format exactly: a line per i, columns are j, values in scientific notation. That is a transposed, header-tagged layout — the opposite of the row-major matrix a plotting tool expects, and the reason the gnuplot script needed a transpose fudge you would have to reproduce. This is the moment the custom format reveals its cost: every consumer must know its private conventions.

Sanity check. The value u(2,1) = 1.0 must appear on the line for i = 2, first column (j = 1): the second data line is 1.00000E+00 4.00000E+00, whose first entry is indeed 1.00000E+00. The format is understood.

Phase 2 — Diagnose: why it can't animate, and why ParaView refuses it

Two independent problems block Friday's animation, and it is worth separating them.

Problem 1 — it is one file, not a series. The run calls dump once, at the end, overwriting field.out each time. There is no record of intermediate states, so there is nothing to animate — the time history was never saved. No format change fixes this alone; the time loop must write a frame per interval.

Problem 2 — the format is not VTK. Even the final snapshot cannot open in ParaView, because field.out has none of the structure a reader needs: no magic string, no declared dataset type, no grid geometry, no POINT_DATA count. ParaView has no way to know it is a 3×2 uniform grid of doubles.

The fix addresses both: write standard VTK, and write it once per interval with an ordered name.

Phase 3 — Recover the geometry the format threw away

VTK needs geometry the custom format never stored: the origin and the spacing. The physics code has them — dx and dy are parameters of the run — they simply were not written to field.out. You confirm the run uses dx = dy = 0.01 on the 3×2 test and a hot top edge. This is the general lesson of §26.1: a custom dump often omits the grid because "everyone here knows it," which is exactly what makes it unportable. VTK forces you to write the geometry down, and that is a feature.

Phase 4 — Port the output to standard legacy VTK

Now replace dump with the canonical write_vtk, giving it the geometry via the field_t type. The physics is untouched; only the output routine changes.

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
  open(newunit=iu, file=filename, status='replace', action='write')
  write(iu, '(a)')             '# vtk DataFile Version 3.0'
  write(iu, '(a, i0)')         'ported from field.out, 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'
  do j = 1, field%ny            ! j outer, i inner: VTK x-fastest order
    do i = 1, field%nx
      write(iu, '(f0.6)') field%u(i, j)
    end do
  end do
  close(iu)
end subroutine write_vtk

Notice the loop nest inverted relative to the inherited dump: the old writer looped i outer and wrote j across; VTK wants j outer, i inner, one value per line. Getting this backwards is exactly the transpose bug of Exercise 10 — the values would all be present but the plate would render sideways. Drive the new writer with the same 3×2 test field (dx = dy = 0.01) and it produces:

# vtk DataFile Version 3.0
ported from field.out, step 0
ASCII
DATASET STRUCTURED_POINTS
DIMENSIONS 3 2 1
ORIGIN 0.000000 0.000000 0.000000
SPACING 0.010000 0.010000 1.000000
POINT_DATA 6
SCALARS temperature double 1
LOOKUP_TABLE default
0.000000
1.000000
2.000000
3.000000
4.000000
5.000000

Sanity check. Six values follow POINT_DATA 6, in VTK order u(1,1),u(2,1),u(3,1),u(1,2),u(2,2), u(3,2) = 0,1,2,3,4,5. Compare with the custom dump from Phase 1: the same six numbers, now in a standard, self-describing container ParaView can read.

Phase 5 — Save a series and open the animation

The last change is in the driver's time loop: call write_vtk every save_every steps with a zero-padded name (using frame_name from Chapter 12).

do step = 0, n_steps
  if (mod(step, save_every) == 0) call write_vtk(field, frame_name(step), step)
  call step_field(field, alpha, dt)
end do

Run it, and the output folder fills with heat_000000.vtk, heat_000100.vtk, … . In ParaView, File ▸ Open collapses them into one heat_..vtk source; Apply, color by temperature, and the play button appears. Friday's animation exists — and the physics code was never touched.

Discussion Questions

1. The inherited dump was not wrong — it worked for its author. What, precisely, made it a liability the moment someone else (or a new tool) needed the data? Relate your answer to the "self-describing" idea from Chapter 25.
2. Porting the writer inverted the loop nest (i outer → j outer). Why is that not just a style choice but a correctness requirement? What would the picture look like if you kept the old order?
3. The geometry (dx, dy) existed in the code but not in field.out. Why does VTK's insistence that you write it down make the data more valuable, not just more verbose?
4. You changed only the output routine, not the physics. Why is that separation (I/O in one module, the solver in another — Chapter 8) what made this a one-afternoon job instead of a rewrite?

Your Turn: Extensions

• Option A (analyze). Take the Phase 1 output and write a Fortran reader that parses field.out back into a field_t (read the 2i6 header, then the transposed values into the right u(i,j)). Confirm it round-trips by writing VTK from the reconstructed field and diffing against the Phase 4 file.
• Option B (port). The old plot.gnu script drew a single heat map. Port its intent to a matplotlib script (§26.5) that reads the ported VTK's values and renders the same map with a perceptually uniform colormap — the modern replacement for the brittle gnuplot script.
• Option C (extend). Add a .pvd collection (§26.3) so the ported series carries physical time (t = step*dt). Present the animation with a real-seconds slider at the group meeting.

Key Takeaways

• Read what a code writes before you trust or replace it — recover the format against a field whose every value you chose, then verify one known entry.
• A custom dump's convenience is a portability debt: no magic string, no geometry, no count means no tool but its author's can read it.
• Porting to VTK is mostly bookkeeping done exactly: the five header parts, the geometry the old format omitted, and the j-outer/i-inner value order.
• Keep I/O separate from physics and a format change touches one routine — the analysis here changed how the run is seen, not what it computes.