Case Study 1: The Simulation That Spent Its Life Formatting

"The fastest I/O is the I/O you do not convert."

Executive Summary

A colleague hands you a working 2D diffusion code. The physics is fine and the results are correct, but the program is slow — a 5,000-step run on a $1000 \times 1000$ grid takes far longer than the arithmetic can explain. You suspect, correctly, that the culprit is not the solver but the output: every step, the code writes the entire temperature field to a text file. In this case study you read the offending routine, diagnose the two costs of text I/O (it is slow, and it is lossy), estimate the damage with a back-of-envelope calculation, and port the output to unformatted binary — turning an I/O-bound program back into a compute-bound one. By the end you will recognize a text-I/O bottleneck on sight, which is one of the most common and most fixable performance bugs in scientific code.

Skills applied: reading and formatting output with edit descriptors (§7.1); the precision and speed cost of text versus unformatted/stream binary (§7.5); robust, guarded I/O with iostat (§7.6); the back-of-envelope habit from the exercises.

Background

The code discretizes a square plate into a $1000 \times 1000$ grid of real(dp) temperatures and marches it forward for 5,000 steps. Each step performs a modest number of floating-point operations per cell — the five-point stencil you will build in Chapter 24 — and then calls a routine dump_text to save the field for later plotting. On the colleague's workstation the whole run takes about nine minutes, and profiling (a skill from Chapter 28) points the finger squarely at dump_text. Our job is to understand why, and fix it, without touching the physics.

Phase 1 — Read the Output Routine

Here is the routine, lightly cleaned up. Read it before you read our commentary.

subroutine dump_text(f, step)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  real(dp),     intent(in) :: f(:,:)
  integer,      intent(in) :: step
  integer :: u, i, j
  character(len=32) :: name

  write(name, '(a, i0, a)') 'field_', step, '.txt'    ! e.g. field_137.txt
  open(newunit=u, file=trim(name), status='replace', action='write')
  do i = 1, size(f, 1)
     do j = 1, size(f, 2)
        write(u, '(es24.16)') f(i, j)                  ! one value per line
     end do
  end do
  close(u)
end subroutine dump_text

Two things jump out. First, it writes one value per write — a million write statements per step, each formatting one double into 24 characters of decimal text. The formatting work alone is enormous. Second, it is called every step, producing 5,000 files totaling gigabytes, almost none of which the colleague will ever look at. The routine is correct, and it is a performance disaster.

The clue that matters: es24.16 is a lossy as well as slow choice — but here it is actually the least lossy text option, keeping all 16 significant digits. The routine's author was being careful about precision, and paid for it with size. That tension — text is either lossy or bulky, never both cheap and exact — is the heart of the problem.

Phase 2 — Estimate the Cost

Put numbers on it before changing anything. A $1000 \times 1000$ field is $10^6$ values. As es24.16 text, each value is about 24 bytes, so the file is roughly 24 MB; as raw binary at 8 bytes per double, it is 8 MB. And the time: a decimal conversion costs on the order of 100 ns per value, a bulk byte copy about 1 ns. The following calculator makes the comparison concrete (round, illustrative figures):

program io_cost
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: n, bin_mb, txt_mb, txt_step, bin_step
  n        = 1000.0_dp * 1000.0_dp        ! 1.0e6 values
  bin_mb   = n * 8.0_dp  / 1.0e6_dp       ! MB as binary
  txt_mb   = n * 24.0_dp / 1.0e6_dp       ! MB as es24.16 text
  txt_step = n * 100.0e-9_dp              ! ~100 ns / value as text
  bin_step = n * 1.0e-9_dp                ! ~1 ns / value as binary

  print '(a, f6.1)', 'binary size (MB)     : ', bin_mb
  print '(a, f6.1)', 'text size (MB)       : ', txt_mb
  print '(a, f6.3)', 'text write / step (s): ', txt_step
  print '(a, f6.3)', 'bin  write / step (s): ', bin_step
  print '(a, f6.1)', 'text 5000 steps (s)  : ', txt_step * 5000.0_dp
  print '(a, f6.1)', 'bin  5000 steps (s)  : ', bin_step * 5000.0_dp
end program io_cost

$ gfortran -std=f2018 -Wall -O2 io_cost.f90 -o cost && ./cost
binary size (MB)     :    8.0
text size (MB)       :   24.0
text write / step (s):  0.100
bin  write / step (s):  0.001
text 5000 steps (s)  :  500.0
bin  5000 steps (s)  :    5.0

There is the nine minutes. Roughly 500 seconds — over eight minutes — of the run is spent converting numbers to text, versus about 5 seconds to write the same data as binary. The output is nearly two orders of magnitude more expensive as text, and it is three times larger on disk. The physics was never the problem.

Phase 3 — The Precision Trap Hiding Underneath

Before we fix the speed, notice the second, quieter defect. Text output is a lossy conversion: the value you read back is only as precise as the digits you wrote. es24.16 happens to keep enough digits to reconstruct a double almost exactly, but the moment a well-meaning author "tidies up" the format to f8.2 for readability, the saved field is rounded to two decimals — and a restart from that file cannot reproduce the original run. A simulation you cannot restart bit-for-bit is a reproducibility problem (Chapter 37). Binary output sidesteps the trap entirely: it stores the exact bits, so a restart continues exactly where it left off.

Phase 4 — Port to Unformatted Binary

The fix is mechanical and safe. Replace the per-value text writes with a single unformatted write of the whole array, and — because the colleague will hand the files to a Python plotting script — use stream access so the bytes are a clean, marker-free sequence NumPy can read.

program port_demo
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: f(3) = [1.5_dp, 2.5_dp, 3.5_dp]   ! stand-in for the field
  real(dp) :: back(3)
  integer  :: u

  ! NEW dump: whole array, one statement, stream binary, exact bytes.
  open(newunit=u, file='field.raw', access='stream', form='unformatted', &
       status='replace', action='write')
  write(u) f
  close(u)

  ! Prove the round trip is exact by reading it straight back.
  open(newunit=u, file='field.raw', access='stream', form='unformatted', &
       status='old', action='read')
  read(u) back
  close(u)

  print '(a, 3f6.2)', 'recovered:', back
end program port_demo

$ gfortran -std=f2018 -Wall -O2 port_demo.f90 -o port && ./port
recovered:  1.50  2.50  3.50

The million per-value writes collapse to one write(u) f; the conversion cost vanishes; the values return exactly. On the Python side the file is read in a line — field = numpy.fromfile('field.raw', dtype='float64').reshape(1000, 1000, order='F') — where the three things the two languages must agree on (dtype float64, the shape, and the column-major order='F') are the caveats §7.5 warned about, and the subject of Chapter 15.

Phase 5 — Write Less, and Write Safely

Speed is not the only lever; frequency is the other. The colleague did not need output every step — every 100th step is plenty for an animation. Make the cadence a namelist parameter (§7.4) so it is tuned without recompiling, and guard the open with iostat (§7.6) so a full disk fails with a message instead of corrupting the run:

! In the driver, read once at startup:
!   namelist /output/ write_every
!   write_every = 100
!   read(cfg_unit, nml=output)
! Then, in the time loop:
if (mod(step, write_every) == 0) then
   call dump_binary(field, step)     ! stream-binary version of the routine
end if

Writing every 100th step instead of every step is a further 100× reduction in output volume and time, on top of the ~100× per-write speedup — the two fixes compound. Sanity check: with write_every = 100, the 5,000-step run now writes 50 binary files of 8 MB each (400 MB total) in a few seconds, versus 5,000 text files totaling ~120 GB over eight minutes. The nine-minute run becomes a few seconds of arithmetic plus a few seconds of I/O — compute-bound again, as it should be.

Discussion Questions

1. The original author chose es24.16 specifically to preserve precision. Explain why that was the right instinct but the wrong tool, and what binary I/O gives them that no text format can.
2. Writing every 100th step reduces output by 100×, but loses the intermediate states. When is that acceptable, and when would a scientist genuinely need every step? (Consider animation versus restart versus debugging a blow-up.)
3. Stream binary is fast and exact but not self-describing — the file records neither its shape nor its dtype. What are the risks of that on a shared cluster, and what does Chapter 25 offer instead?

Your Turn: Extensions

• Option A. Take the dump_text routine and measure it for real: write a $500 \times 500$ field once as es24.16 text and once as stream binary, time both with system_clock (Chapter 28), and compare file sizes with ls -l. Do your measured numbers match the order-of-magnitude estimate here?
• Option B. Add a text header to the binary file (write nx and ny as the first two records) so the file becomes minimally self-describing, and write the matching reader that uses the header to size its array before reading the data.
• Option C. Read one of your binary fields back in Python with numpy.fromfile, reshape it with order='F', and plot it with matplotlib. Confirm the plate looks right — and deliberately use order='C' once to see the transposed garbage that a row/column-major mismatch produces.

Key Takeaways

• A program that writes text every step is almost always I/O-bound for a reason you can eliminate: text conversion is ~100× slower than a raw byte copy and several times larger on disk.
• Text I/O forces a choice between lossy (f8.2) and bulky (es24.16); binary is exact and compact, which is why checkpoints and restart files are binary.
• The two fixes compound: write each field as one unformatted/stream write, and write fewer fields via a namelist-configured cadence. Together they turn minutes back into seconds.
• Diagnosing an I/O bottleneck is a reading skill: spot the per-value write in a hot loop, estimate its cost on the back of an envelope, and you have found the bug before you profile.