Chapter 26 Exercises: Visualization Output

These exercises take you from writing a single valid VTK file to producing an animated time series and a publication figure — the full pipeline from the array in memory to the picture on the page. Work them at a real terminal with gfortran, ParaView (or VisIt), and a Python with NumPy and matplotlib installed; seeing your own field animate is the point.

Difficulty tiers. ⭐ builds a single skill in isolation; ⭐⭐ combines two or three; ⭐⭐⭐ is a small project that integrates the chapter (and often an earlier one). Solutions to problems marked † and to all odd-numbered problems are in appendices/answers-to-selected.md; the rest are yours to check by compiling and by opening the result in ParaView.

A reminder that runs through every problem: never trust "it compiled" — open the file. A VTK writer can compile perfectly and still emit a header that ParaView rejects or, worse, silently misreads. Predict the file's exact contents first, then compile, then open it in a text editor and in ParaView.

Part A — Type, Compile, and Run (predict the output first)

1. † ⭐ Compile and run example-01-legacy-vtk.f90. Before you look at heat_000100.vtk, write down — on paper — exactly what its DIMENSIONS, POINT_DATA, and the six value lines will be. Then open the file and confirm every line. Which single line is the "magic string" that makes it a VTK file at all?

2. ⭐ Change the demo field in example-01 to a 2×2 plate with values u(:,1) = [10, 20], u(:,2) = [30, 40] and dx = dy = 1.0. Predict the complete file, then verify. What does POINT_DATA become, and in what order do the four values appear?

3. † ⭐⭐ Add a second scalar array to the legacy VTK writer. After the temperature block, append a second SCALARS array named error (also double 1) holding field%u - u_exact for some analytic u_exact. Predict what ParaView's coloring dropdown will offer after you open the file. (Hint: multiple SCALARS blocks may follow one POINT_DATA line.)

4. ⭐ Run example-03-simple-formats.f90 and predict both heat.dat and heat_gp.dat before opening them. In heat.dat, why is the first line 3.000 4.000 5.000 and not 0.000 1.000 2.000?

5. † ⭐⭐ Type, compile, and run a program that writes the same 3×2 field as both legacy .vtk and XML .vti (call write_vtk and write_vti). Open both in ParaView. Confirm they render identically, then diff the two files and list three structural differences between the formats.

Part B — Port It (translate to a Fortran output routine)

6. ⭐⭐ Here is a Python routine that dumps a field to legacy VTK. Port it to a Fortran subroutine with the canonical signature write_vtk(field, filename, step), matching its output byte for byte.

def write_vtk(u, dx, dy, filename):
    ny, nx = u.shape
    with open(filename, "w") as f:
        f.write("# vtk DataFile Version 3.0\n")
        f.write("ported from python\n")
        f.write("ASCII\n")
        f.write("DATASET STRUCTURED_POINTS\n")
        f.write(f"DIMENSIONS {nx} {ny} 1\n")
        f.write("ORIGIN 0 0 0\n")
        f.write(f"SPACING {dx} {dy} 1\n")
        f.write(f"POINT_DATA {nx*ny}\n")
        f.write("SCALARS temperature double 1\nLOOKUP_TABLE default\n")
        for j in range(ny):
            for i in range(nx):
                f.write(f"{u[j, i]:.6f}\n")

Watch the index order carefully: NumPy's u[j, i] is row j, column i. What must the Fortran loop be so the same value lands at the same point? (This is the row-major/column-major trap from Chapter 15 in disguise.)

7. † ⭐⭐ Port this MATLAB one-liner, which writes a field as a comma-separated matrix for a spreadsheet, to a Fortran subroutine write_csv(field, filename). Use the unlimited-repeat format from §26.4 and a comma separator.

writematrix(flipud(U), 'heat.csv')   % flipud puts the top row first

8. ⭐⭐⭐ Port a small gnuplot-driving workflow: given a Python script that saves a field and shells out to gnuplot, reproduce the Fortran half (write the 3-column blank-line-separated .dat) and the gnuplot script half (the three commands that draw the pm3d heat map). You are porting the data-writing to Fortran and leaving the plotting to gnuplot — the division of labor §26.4 argues for.

Part C — Find the Bug

9. † ⭐ This writer compiles and runs, but ParaView reports "Error reading ASCII data. Expected 12 points but could only read 6." Find and fix the bug.

write(iu, '(a, 3(1x, i0))') 'DIMENSIONS', field%nx, field%ny, 2   ! nx=3, ny=2
write(iu, '(a, 1x, i0)')    'POINT_DATA', field%nx * field%ny

10. ⭐⭐ A colleague's plate comes out transposed in ParaView — rows and columns swapped — though the values are all present. The header is correct. Here is the value loop:

do i = 1, field%nx
  do j = 1, field%ny
    write(iu, '(f0.6)') field%u(i, j)
  end do
end do

Explain why the picture is transposed and give the one-line fix. What does this cost in memory-access terms (tie your answer to the Performance Note in §26.2)?

11. † ⭐ This first line was meant to be the VTK magic string but ParaView refuses to open the file. What is wrong, and why is VTK so unforgiving here?

write(iu, '(a)') '#vtk DataFile Version 3.0'

12. ⭐⭐ A time series loads into ParaView but plays out of order: step 1000 appears before step 200. The frames are named with this helper. Diagnose the bug and fix it.

write(buf, '(a, i0, a)') 'heat_', step, '.vtk'   ! step 200 -> 'heat_200.vtk'
name = trim(buf)

13. † ⭐⭐ This .vti writer produces a file ParaView opens but shows as empty (no data to color by). The extent and geometry are right. What is missing or wrong in the DataArray/PointData block?

write(iu,'(a)') '      <PointData>'
write(iu,'(a)') '        <DataArray type="Float64" format="ascii">'

Part D — Modernize It

14. ⭐⭐ Here is a genuine-looking FORTRAN 77-era output routine that dumps a field with a bespoke ASCII format and a hard-coded unit number. Modernize it: implicit none, free-form, newunit, intent, assumed-shape argument, kind-parameterized reals — and make it write standard legacy VTK instead of the ad-hoc format, so it opens in ParaView.

      SUBROUTINE DUMP(U, NX, NY)
      DIMENSION U(NX,NY)
      OPEN(7, FILE='dump.txt')
      WRITE(7,*) NX, NY
      DO 10 J=1,NY
      DO 10 I=1,NX
      WRITE(7,*) U(I,J)
   10 CONTINUE
      CLOSE(7)
      RETURN
      END

15. † ⭐⭐ A legacy code writes one giant text file per run with every timestep concatenated (no way to scrub to a frame). Modernize the strategy: describe (and sketch in code) how to split it into a zero-padded per-step VTK time series, and explain what capability in ParaView this unlocks that the single concatenated file cannot offer.

Part E — Design It (extend the heat solver)

16. ⭐⭐ Extend the Project Checkpoint driver to also accumulate a heat.pvd collection file (§26.3) as it saves frames, so ParaView shows the slider in physical seconds. Compute each frame's time as step * dt. Sketch the three writing moments: the header (once, before the loop), one <DataSet .../> line per saved frame, and the two closing tags (once, after the loop).

17. † ⭐⭐⭐ Design a write_output dispatcher for the solver: a single subroutine write_output(field, step, fmt) where fmt is a character(len=*) selecting 'vtk', 'vti', 'matrix', or 'gnuplot', dispatching to the right writer with a select case. Give it a sensible default and an error stop on an unknown format. Which chapter's construct is select case (Chapter 4)?

18. ⭐⭐ Add a write_centerline(field, filename) routine that writes a two-column x u file of the temperature along the horizontal centerline (j = ny/2), for a quick line plot in gnuplot or matplotlib. This is the 1-D cross-section scientists reach for constantly. Predict its length for nx = 101.

19. † ⭐⭐⭐ Design the 3-D generalization of write_vtk: a field u(:,:,:) on an nx × ny × nz grid. What changes in DIMENSIONS, POINT_DATA, the loop nest, and the SPACING line? Write the routine and state the VTK point ordering for three dimensions.

Part F — Back of the Envelope

20. ⭐ A run saves a 256 × 256 field as ASCII legacy VTK every 50 steps for 50,000 steps. About how many frames is that, and — taking each value as roughly 9 ASCII bytes (0.123456\n) — about how large is the whole time series? Would you keep it as ASCII VTK, or reach for the binary formats of Chapter 25?

21. † ⭐⭐ The same 256 × 256 field in double-precision binary is 256*256*8 bytes per frame. Compare the per-frame size of ASCII VTK (from problem 20) with binary, and state the two costs ASCII pays for its human-readability. Which cost matters more at a billion cells?

22. ⭐⭐ You want a smooth 10-second animation at 30 frames per second. Your run is 90,000 timesteps. What save_every gives you enough frames, and how many files is that? (Round to a convenient value.)

23. † ⭐ Estimate: writing one 1000 × 1000 ASCII VTK frame is a million write statements of ~9 bytes. If your disk sustains 500 MB/s, is the write more likely limited by disk bandwidth or by the per-write formatting overhead? (You need not be exact — reason about the order of magnitude and name the likely bottleneck.)

Part G — Interleaved (earlier chapters resurface)

24. ⭐⭐ (Ch. 12) The frame_name helper uses i6.6, which holds up to 999,999 steps. Your run needs 2 million steps. Rewrite the format so filenames stay sortable, and state the new fixed width. What happens if you don't widen it — what does i6.6 write for step 1,000,000?

25. † ⭐⭐ (Ch. 7) Explain, using the edit descriptors from Chapter 7, exactly what write(iu, '(a, 3(1x, i0))') 'DIMENSIONS', 100, 80, 1 places in the record, character group by character group. Why i0 and not i4?

26. ⭐⭐ (Ch. 9) The writers all take a field_t carrying nx, ny, dx, dy, u(:,:). Why is passing the single derived type better than passing five separate arguments write_vtk(u, nx, ny, dx, dy, ...)? Tie your answer to the derived-type motivation of Chapter 9.

27. † ⭐⭐⭐ (Ch. 5, Ch. 24) The value-writing loop do j; do i; write u(i,j) walks the array in column-major order. (a) Explain why that is the cache-friendly order for a Fortran array. (b) The Chapter 24 five-point stencil update u(i,j) = f(u(i±1,j), u(i,j±1)) — does it also want i innermost for the same reason? (c) State the general rule connecting VTK's required write order, Fortran's memory layout, and loop nesting.

28. ⭐⭐ (Ch. 24) Your solver visualizes beautifully for a while, then every frame goes uniformly to a huge value and ParaView shows a single flat color. The output code is fine. What numerical failure from Chapter 24 does the visualization just make visible, and what does the picture look like at the instant it goes wrong?