Appendix I: The Complete Heat-Solver Code

This appendix reproduces, in full, the program you built across the book: the 2D heat-equation solver, assembled as one coherent, compilable codebase. Each Project Checkpoint added one piece — the dp kind in Chapter 3, the field_t type in Chapter 9, the finite-difference core in Chapter 24, VTK output in Chapter 26, timing in Chapter 28, the src/+app/+test/ layout in Chapter 36, the tests in Chapter 37, and the final assembly in the Chapter 38 capstone. Here they are gathered in one place.

Everything below compiles with gfortran -std=f2018 -Wall (the OpenMP variant in §I.6 needs -fopenmp); see Appendix C for the toolchain and the flag reference. As everywhere in this book, no code was executed to produce the outputs — every ! Expected output: was computed by hand. The solver's signature computation, the 5×5 plate marched two steps, is reproduced digit-for-digit by every version here: 28 / 32 / 28, 4 / 4 / 4, maximum 100.

The one idea that held it together. Through every rewrite — serial array update, OpenMP parallel do, MPI halo exchange — the public signatures never moved: step(field, alpha, dt), laplacian(u, dx, dy), write_vtk(field, filename, step). The interface is the contract; the body is free to change. That is why forty chapters could keep bolting capability onto one program without ever disturbing the driver, the I/O, or the tests.

I.1 The directory tree

The solver is an fpm package: a thin app/ driver on top, the library in src/, tests in test/, and a manifest that lets fpm derive the compile order from the use statements. Read it bottom-up and it is a topological sort you can see — kinds.f90 depends on nothing, main.f90 depends on everything.

heat-solver/
├── fpm.toml                     the build manifest (name, version, dependencies)
├── README.md                    what it solves; how to build, run, reproduce
├── heat.nml                     the input deck (run parameters, read by read_config)
├── src/                         THE LIBRARY
│   ├── kinds.f90                    UTILITY: the dp precision kind             (Ch. 3)
│   ├── heat_types.f90               CORE:    the field_t derived type          (Ch. 9)
│   ├── heat_solver.f90              SOLVER:  laplacian + step + stable_dt       (Ch. 24)
│   ├── heat_io.f90                  I/O:     read_config, write_field, write_vtk, frame_name (Ch. 7, 26)
│   └── timers.f90                   UTILITY: tic/toc around system_clock        (Ch. 28)
├── app/
│   └── main.f90                     DRIVER:  program heat — setup, time loop, output
└── test/
    └── test_solver.f90              regression + verification vs the analytical solution (Ch. 37)

I.2 The modules, in full

Presented in dependency order — kinds → heat_types → heat_solver → heat_io → timers → main. Each is one .f90 file under src/ (the driver under app/). A module on its own produces no output, so only the driver (§I.2.6) carries an ! Expected output:.

I.2.1 src/kinds.f90

The precision foundation. One line of real content — dp = selected_real_kind(15, 307) requests a real kind with at least 15 significant decimal digits and an exponent range to $10^{307}$, which on every mainstream machine is IEEE 754 double precision. Defining it once, here, and writing real(dp) and 1.0_dp everywhere else is what keeps the whole codebase at one, deliberate, portable precision. This file has survived unchanged since Chapter 3.

! src/kinds.f90 -- the single precision kind for the whole project (Ch. 3).
module kinds
  implicit none
  private
  public :: dp
  integer, parameter :: dp = selected_real_kind(15, 307)   ! IEEE double: ~15-16 digits
end module kinds

I.2.2 src/heat_types.f90

The central data structure. A field_t bundles the grid geometry (nx, ny, dx, dy) with the temperature field u(:,:) as an allocatable component, so one object carries everything a solver routine needs instead of five loose arguments. The type-bound init sizes and zeroes the field; because its self is class(field_t), intent(out), calling init on an already-used field cleanly deallocates and resets it. Every routine downstream takes a type(field_t).

! src/heat_types.f90 -- the field_t derived type: grid + temperature field (Ch. 9).
module heat_types
  use kinds, only: dp
  implicit none
  private
  public :: field_t

  type :: field_t
    integer  :: nx = 0, ny = 0                  ! grid points in x (first index) and y
    real(dp) :: dx = 0.0_dp, dy = 0.0_dp        ! grid spacing
    real(dp), allocatable :: u(:,:)             ! the temperature field, sized at run time
  contains
    procedure :: init => field_init
  end type field_t

contains

  subroutine field_init(self, nx, ny, dx, dy)
    class(field_t), intent(out) :: self         ! class(...), intent(out) resets/reallocates
    integer,  intent(in) :: nx, ny
    real(dp), intent(in) :: dx, dy
    self%nx = nx;  self%ny = ny
    self%dx = dx;  self%dy = dy
    allocate(self%u(nx, ny))
    self%u = 0.0_dp
  end subroutine field_init
end module heat_types

I.2.3 src/heat_solver.f90

The numerical heart, from Chapter 24. Three public procedures, all with frozen signatures:

• laplacian(u, dx, dy) — the scaled five-point stencil. It approximates $\nabla^2 u$ at each interior point as the sum of two central second differences, each divided by the spacing squared: $\nabla^2 u|_{i,j} \approx (u_{i-1,j} - 2u_{i,j} + u_{i+1,j})/\Delta x^2 + (u_{i,j-1} - 2u_{i,j} + u_{i,j+1})/\Delta y^2$. It is second-order accurate ($O(h^2)$) and exact for quadratics — the basis of the unit test in §I.7. Boundary rows and columns are returned as zero; they are held fixed and never enter the update.
• step(field, alpha, dt) — one explicit forward-Euler (FTCS) step. It evaluates the Laplacian on the old field, then advances the interior by $\alpha\,\Delta t\,\nabla^2 u$. The Dirichlet edges are held.
• stable_dt(alpha, dx, dy, safety) — a CFL-safe timestep. The explicit scheme is only conditionally stable: the diffusion number $r = \alpha\,\Delta t / h^2$ must satisfy $r \le \tfrac14$ in 2D, i.e. $\Delta t \le h^2/(4\alpha)$. This returns safety times that limit (default $0.9$).

! src/heat_solver.f90 -- the finite-difference core: five-point stencil, FTCS, CFL (Ch. 24).
module heat_solver
  use kinds,      only: dp
  use heat_types, only: field_t
  implicit none
  private
  public :: laplacian, step, stable_dt

contains

  ! Scaled five-point Laplacian; interior filled, boundary rows/cols left 0. Second-order
  ! accurate, and EXACT for fields polynomial of degree <= 3.
  pure function laplacian(u, dx, dy) result(lap)
    real(dp), intent(in) :: u(:,:), dx, dy
    real(dp) :: lap(size(u,1), size(u,2))
    integer  :: nx, ny
    nx = size(u,1);  ny = size(u,2)
    lap = 0.0_dp
    lap(2:nx-1,2:ny-1) =                                                             &
          ( u(1:nx-2,2:ny-1) - 2.0_dp*u(2:nx-1,2:ny-1) + u(3:nx,2:ny-1) ) / dx**2   &
        + ( u(2:nx-1,1:ny-2) - 2.0_dp*u(2:nx-1,2:ny-1) + u(2:nx-1,3:ny) ) / dy**2
  end function laplacian

  ! One explicit forward-Euler (FTCS) step. Frozen signature step(field, alpha, dt).
  ! Evaluate the Laplacian on the OLD field (a clean snapshot), then advance the interior;
  ! the Dirichlet boundary rows/cols are never written.
  subroutine step(field, alpha, dt)
    type(field_t), intent(inout) :: field
    real(dp),      intent(in)    :: alpha, dt
    real(dp), allocatable :: lap(:,:)
    integer :: nx, ny
    nx = field%nx;  ny = field%ny
    lap = laplacian(field%u, field%dx, field%dy)
    field%u(2:nx-1,2:ny-1) = field%u(2:nx-1,2:ny-1)  &
                           + alpha*dt * lap(2:nx-1,2:ny-1)
  end subroutine step

  ! A CFL-safe timestep. 2D limit: alpha*dt*(1/dx^2 + 1/dy^2) <= 1/2, i.e. dt <= h^2/(4a)
  ! on a square grid. 'safety' is the fraction of that limit to take (default 0.9).
  pure function stable_dt(alpha, dx, dy, safety) result(dt)
    real(dp), intent(in)           :: alpha, dx, dy
    real(dp), intent(in), optional :: safety
    real(dp) :: dt, s
    s = 0.9_dp
    if (present(safety)) s = safety
    dt = s / (2.0_dp*alpha*(1.0_dp/dx**2 + 1.0_dp/dy**2))
  end function stable_dt
end module heat_solver

⚡ Performance / correctness note. The interior update is written as array sections, not nested loops. The compiler vectorizes the section assignment and — because Fortran forbids argument aliasing — is free to do so aggressively. The serial step here and the OpenMP step in §I.6 are algebraically identical (same stencil, same $r$) and produce the identical field, digit for digit, as Chapter 33 proves.

I.2.4 src/heat_io.f90

All contact with the outside world, gathered from Chapters 7 and 26. read_config parses the run parameters from a Fortran namelist on an already-open unit (the mechanism is Chapter 7's; the driver opens heat.nml and passes the unit). write_field writes the field as text, one grid row per line. write_vtk writes a legacy VTK STRUCTURED_POINTS frame that ParaView or VisIt opens directly — note the value loop runs the first index i fastest, matching column-major memory order. frame_name builds a zero-padded filename with an internal-file write: frame_name(123) returns "heat_000123.vtk".

! src/heat_io.f90 -- configuration in, results out: namelist, text, and VTK (Ch. 7, 26).
module heat_io
  use kinds,      only: dp
  use heat_types, only: field_t
  implicit none
  private
  public :: read_config, write_field, write_vtk, frame_name

contains

  ! Read run parameters from a namelist on an already-open unit (Ch. 7-8). The caller
  ! opens the deck and supplies its unit; the group name in the file must be &config.
  subroutine read_config(unit, nx, ny, alpha, dt, n_steps)
    integer,  intent(in)  :: unit
    integer,  intent(out) :: nx, ny, n_steps
    real(dp), intent(out) :: alpha, dt
    namelist /config/ nx, ny, alpha, dt, n_steps
    read(unit, nml=config)
  end subroutine read_config

  ! Write the field as text, one grid row per line. The unlimited format item *(f8.2)
  ! repeats f8.2 across however many columns the row has.
  subroutine write_field(field, filename)
    type(field_t),    intent(in) :: field
    character(len=*), intent(in) :: filename
    integer :: unit, i, ios
    open(newunit=unit, file=filename, status='replace', action='write', iostat=ios)
    if (ios /= 0) then
       print '(a)', 'write_field: cannot open ' // trim(filename)
       return
    end if
    do i = 1, field%nx
       write(unit, '(*(f8.2))') field%u(i, :)
    end do
    close(unit)
  end subroutine write_field

  ! Write a legacy VTK STRUCTURED_POINTS frame for ParaView/VisIt (Ch. 26). The nine
  ! header lines are mandatory and their order and spelling are exact; the values follow
  ! with the first index i varying fastest (column-major order).
  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'
    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'
    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

  ! A zero-padded frame filename: frame_name(123) -> "heat_000123.vtk" (Ch. 12, 26).
  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'
    name = trim(buf)
  end function frame_name
end module heat_io

🔗 Config as a derived type. read_config here returns five loose intent(out) scalars, exactly as the reader built it in Chapter 8. A natural next refactor — foreshadowed there and used in this chapter's case studies — is to bundle them into a config_t derived type and read them into it (read_config(unit, cfg)), so the run parameters travel as one object the way the field does. The namelist body is unchanged; only the argument list tightens.

I.2.5 src/timers.f90

Wall-clock timing, from Chapter 28. tic starts a stopwatch and toc returns the elapsed seconds, built on the intrinsic system_clock. The tick rate is read once and cached in module state; the counts are integer(int64) so a fast, high-resolution clock does not overflow mid-run. Wall time, not CPU time, is what you want for a parallel solver — it is the time the user actually waits.

! src/timers.f90 -- tic/toc wall-clock timing around system_clock (Ch. 28).
module timers
  use kinds, only: dp
  use, intrinsic :: iso_fortran_env, only: int64
  implicit none
  private
  public :: tic, toc

  integer(int64) :: start_count = 0_int64      ! module state (saved by design)
  integer(int64) :: rate        = 1_int64
  logical        :: have_rate   = .false.

contains

  subroutine tic()                             ! start (or restart) the stopwatch
    if (.not. have_rate) then
      call system_clock(count_rate=rate)       ! read the tick rate exactly once
      have_rate = .true.
    end if
    call system_clock(count=start_count)
  end subroutine tic

  function toc() result(seconds)               ! elapsed wall-clock seconds since tic()
    real(dp) :: seconds
    integer(int64) :: now
    call system_clock(count=now)
    seconds = real(now - start_count, dp) / real(rate, dp)
  end function toc
end module timers

I.2.6 app/main.f90 — the driver

The program heat driver orchestrates the library and contains no physics of its own: it reads the deck, builds the plate, marches the time loop while writing a VTK frame every save_every steps, times the run, and reports. Its printed output is fully determinate — the configuration it read, the diffusion number, the frame log, and the invariant maximum temperature (the hot edge is held at 100, and $r \le \tfrac14$ keeps the interior in $[0,100]$ by the discrete maximum principle) — with only the elapsed time being machine-dependent.

! app/main.f90 -- program heat: setup, time loop, periodic VTK output, timing.
! Build & run in the package:  fpm run   (compile order for a plain build is in section I.4).
! NEVER run during authoring; the output below is hand-computed.
program heat
  use kinds,       only: dp
  use heat_types,  only: field_t
  use heat_solver, only: step, stable_dt
  use heat_io,     only: read_config, write_field, write_vtk, frame_name
  use timers,      only: tic, toc
  implicit none

  type(field_t)      :: f
  integer            :: nx, ny, n_steps, n, unit, ios
  integer, parameter :: save_every = 100          ! write a VTK frame this often
  real(dp)           :: alpha, dt, r, secs
  character(len=200) :: msg

  ! [1] Read the run parameters from the namelist deck heat.nml (Ch. 7-8).
  open(newunit=unit, file='heat.nml', status='old', action='read', iostat=ios, iomsg=msg)
  if (ios /= 0) error stop 'heat: cannot open heat.nml: ' // trim(msg)
  call read_config(unit, nx, ny, alpha, dt, n_steps)
  close(unit)

  ! [2] Build the plate. Unit grid spacing keeps r equal to the book's worked value 0.2
  !     (a physical run would set dx = L/(nx-1)). Hot top edge held at 100 (Dirichlet).
  call f%init(nx, ny, dx=1.0_dp, dy=1.0_dp)
  f%u(1,:) = 100.0_dp

  ! [3] Guard the CFL limit: if the deck's dt is unstable, fall back to a safe one.
  r = alpha*dt / f%dx**2
  if (r > 0.25_dp) then
     dt = stable_dt(alpha, f%dx, f%dy, safety=0.9_dp)
     r  = alpha*dt / f%dx**2
  end if

  print '(a)',            'heat-solver: 2D explicit heat equation (FTCS)'
  print '(a, i0, a, i0)', 'grid        : ', nx, ' x ', ny
  print '(a, es9.2)',     'alpha       : ', alpha
  print '(a, f6.3)',      'dt          : ', dt
  print '(a, i0)',        'n_steps     : ', n_steps
  print '(a, f6.3, a)',   'r = a*dt/h^2: ', r, '  (<= 0.25, stable)'

  ! [4] March in time; write a VTK frame every save_every steps (Ch. 26).
  call tic()
  do n = 0, n_steps
     if (mod(n, save_every) == 0) then
        call write_vtk(f, frame_name(n), n)
        print '(a)', 'wrote ' // frame_name(n)
     end if
     if (n < n_steps) call step(f, alpha, dt)     ! exactly n_steps updates
  end do
  secs = toc()

  ! [5] Final text field + a determinate summary.
  call write_field(f, 'heat_final.txt')
  print '(a, i0, a)',     'done: ', n_steps, ' steps'
  print '(a, f8.2)',      'max temperature : ', maxval(f%u)   ! = 100.00 (held hot edge)
  print '(a, es12.4, a)', 'elapsed         : ', secs, ' s   (machine-dependent)'
  print '(a)',            'wrote heat_final.txt'
end program heat

! Hand computation (heat.nml sets nx=ny=101, alpha=1.0, dt=0.2, n_steps=500; dx=dy=1.0):
!   r = alpha*dt/dx^2 = 1.0*0.2/1.0 = 0.2, not > 0.25 -> no override.
!   es9.2 of 1.0 -> " 1.00E+00";  f6.3 of 0.2 -> " 0.200".
!   Frames at n = 0,100,200,300,400,500 (save_every=100) -> six "wrote ..." lines.
!   Row 1 is held at 100 for all steps, so maxval(u) = 100.00 exactly; f8.2 -> "  100.00".
!
! Expected output (all lines determinate except the elapsed time):
! heat-solver: 2D explicit heat equation (FTCS)
! grid        : 101 x 101
! alpha       :  1.00E+00
! dt          :  0.200
! n_steps     : 500
! r = a*dt/h^2:  0.200  (<= 0.25, stable)
! wrote heat_000000.vtk
! wrote heat_000100.vtk
! wrote heat_000200.vtk
! wrote heat_000300.vtk
! wrote heat_000400.vtk
! wrote heat_000500.vtk
! done: 500 steps
! max temperature :   100.00
! elapsed         :   4.5000E-02 s   (machine-dependent)
! wrote heat_final.txt

I.3 The input deck heat.nml

read_config parses a Fortran namelist whose group is named config. The deck is order-free, case-insensitive, and any parameter it omits keeps the default the driver set. This is the exact deck the Expected output in §I.2.6 assumes.

&config
  nx      = 101,
  ny      = 101,
  alpha   = 1.0,
  dt      = 0.2,
  n_steps = 500,
/

I.4 fpm.toml and building

The manifest is short because fpm uses convention over configuration: it discovers every .f90 under src/, app/, and test/, reads the use statements, and compiles in dependency order for you.

name = "heat-solver"
version = "1.0.0"
license = "MIT"
author = "Your Name"
maintainer = "you@example.com"

[build]
auto-executables = true
auto-tests = true

With that in place, three commands do everything:

$ fpm build      # compile the library, the app, and the tests, in dependency order
$ fpm run        # build if needed, then run app/main.f90 (reads ./heat.nml)
$ fpm test       # build if needed, then run test/test_solver.f90

Without fpm, compile the sources by hand in dependency order — a module must appear before anything that uses it:

$ gfortran -std=f2018 -Wall -O2 \
    src/kinds.f90 src/heat_types.f90 src/heat_solver.f90 \
    src/heat_io.f90 src/timers.f90 app/main.f90 -o heat && ./heat

For a debug build add -fcheck=all -g -O0; for production use -O2 (or -O3 -march=native -flto) and drop -fcheck=all. See Appendix C for the full flag table.

I.5 What a VTK frame looks like

write_vtk produces a plain-text VTK file ParaView opens directly. Here is a tiny illustration — one frame of a 3×2 plate — so the format is concrete. In the package this program compiles alongside the library (fpm build), because it uses the same write_vtk shown in §I.2.4.

! A one-frame illustration; uses write_vtk from src/heat_io.f90.
program vtk_demo
  use kinds,      only: dp
  use heat_types, only: field_t
  use heat_io,    only: write_vtk
  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 vtk_demo

! Expected output (stdout):
! wrote heat_000100.vtk

The file it writes, byte for byte (the f0.6 leading-zero form is gfortran 10+ / Fortran-2018 behaviour for a zero field width):

# 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

I.6 The OpenMP-parallel step

The frozen step(field, alpha, dt) interface never changes; its body gains one !$omp parallel do over the interior grid. Scoping is the whole lesson: the fields are shared (every thread reads the same old field%u and writes distinct cells of u_new, so there is no race), and the loop indices i, j are private — the inner index i is not auto-private and must be listed, or it races. Static scheduling suits the uniform per-cell work. Compile with -fopenmp; without it the !$omp lines are plain comments and you get correct serial code. The result is deterministic and identical to §I.2 on any thread count.

In the package you do not add a second module — you replace the body of step inside src/heat_solver.f90 with the one below and rebuild with -fopenmp (this is exactly what the Chapter 38 capstone ships). It is shown here as a complete, self-contained program so you can compile and run it directly.

! An OpenMP build of the solver's update. Compile:
!   gfortran -std=f2018 -Wall -fopenmp appendix-i-openmp.f90 -o heat_omp
!   OMP_NUM_THREADS=4 ./heat_omp
! NEVER run during authoring; the output is hand-computed (identical to section I.2).

module kinds
  implicit none
  private
  public :: dp
  integer, parameter :: dp = selected_real_kind(15, 307)
end module kinds

module heat_types
  use kinds, only: dp
  implicit none
  private
  public :: field_t
  type :: field_t
    integer  :: nx = 0, ny = 0
    real(dp) :: dx = 0.0_dp, dy = 0.0_dp
    real(dp), allocatable :: u(:,:)
  contains
    procedure :: init => field_init
  end type field_t
contains
  subroutine field_init(self, nx, ny, dx, dy)
    class(field_t), intent(out) :: self
    integer,  intent(in) :: nx, ny
    real(dp), intent(in) :: dx, dy
    self%nx = nx;  self%ny = ny
    self%dx = dx;  self%dy = dy
    allocate(self%u(nx, ny))
    self%u = 0.0_dp
  end subroutine field_init
end module heat_types

module heat_solver                             ! the drop-in body of step (Ch. 33)
  use kinds,      only: dp
  use heat_types, only: field_t
  implicit none
  private
  public :: step
contains
  subroutine step(field, alpha, dt)            ! frozen signature (Ch. 6)
    type(field_t), intent(inout) :: field
    real(dp),      intent(in)    :: alpha, dt
    real(dp) :: u_new(field%nx, field%ny)
    real(dp) :: rx, ry
    integer  :: i, j, nx, ny
    nx = field%nx;  ny = field%ny
    rx = alpha*dt / field%dx**2                ! diffusion numbers (Ch. 24)
    ry = alpha*dt / field%dy**2
    u_new = field%u                            ! copy keeps the Dirichlet edges
    !$omp parallel do default(none)                        &
    !$omp   shared(field, u_new, nx, ny, rx, ry)           &
    !$omp   private(i, j) schedule(static)
    do j = 2, ny-1                             ! outer index: auto-private
      do i = 2, nx-1                           ! inner index: MUST be private
        u_new(i,j) = field%u(i,j)                                              &
                   + rx*(field%u(i-1,j) - 2.0_dp*field%u(i,j) + field%u(i+1,j)) &
                   + ry*(field%u(i,j-1) - 2.0_dp*field%u(i,j) + field%u(i,j+1))
      end do
    end do
    !$omp end parallel do
    field%u = u_new                            ! commit the whole new field at once
  end subroutine step
end module heat_solver

program heat_omp
  use kinds,       only: dp
  use heat_types,  only: field_t
  use heat_solver, only: step
  implicit none
  type(field_t) :: plate
  real(dp), parameter :: alpha = 1.0_dp, dt = 0.2_dp   ! r = alpha*dt/dx^2 = 0.2 <= 1/4
  integer :: s

  call plate%init(nx=5, ny=5, dx=1.0_dp, dy=1.0_dp)
  plate%u(1,:) = 100.0_dp                              ! top edge (row 1) held at 100

  do s = 1, 2
    call step(plate, alpha, dt)
    print '(a,i0)', 'after step ', s
    call show(plate%u)
  end do

contains
  subroutine show(a)
    real(dp), intent(in) :: a(:,:)
    integer :: i
    do i = 1, size(a,1)
      print '(5f8.2)', a(i,:)
    end do
    print '(a)', ''
  end subroutine show
end program heat_omp

! Hand check (rx = ry = 0.2). Step 1, cell (2,2): 0 + 0.2*(100 - 0 + 0) + 0.2*(0) = 20.
! Step 2, (2,2): 20 + 0.2*(100 - 40 + 0) + 0.2*(0 - 40 + 20) = 20 + 12 - 4 = 28.
! Step 2, (2,3): 20 + 0.2*(100 - 40 + 0) + 0.2*(20 - 40 + 20) = 20 + 12 + 0 = 32.
! Step 2, (3,2): 0  + 0.2*(20 - 0 + 0)  + 0.2*(0)            = 4. Matches section I.2, any thread count.
!
! Expected output:
! after step 1
!   100.00  100.00  100.00  100.00  100.00
!     0.00   20.00   20.00   20.00    0.00
!     0.00    0.00    0.00    0.00    0.00
!     0.00    0.00    0.00    0.00    0.00
!     0.00    0.00    0.00    0.00    0.00
!
! after step 2
!   100.00  100.00  100.00  100.00  100.00
!     0.00   28.00   32.00   28.00    0.00
!     0.00    4.00    4.00    4.00    0.00
!     0.00    0.00    0.00    0.00    0.00
!     0.00    0.00    0.00    0.00    0.00

I.7 The test: test/test_solver.f90

The suite uses three kinds of oracle, each an independent check on the solver (Chapter 37). The unit tests exploit facts about the stencil that hold to machine precision: the Laplacian of $x^2+y^2$ is exactly $4$ everywhere (the stencil is exact for quadratics), the Laplacian of a linear field is exactly $0$, and stable_dt returns $r \le \tfrac14$. The regression test reproduces the golden 5×5 two-step field — 28 / 32 / 28, 4 / 4 / 4 — and compares within a tolerance (portable across compilers and flags). The verification tests use analytical truths: a linear temperature ramp is an exact fixed point of step, and the discrete maximum principle bounds every value in $[0,100]$. Finally the suite reports the one-step decay of an analytical sine mode against its exact amplitude $e^{-2\alpha\pi^2\Delta t}$; the small residual is the coarse-grid $O(h^2)$ discretization error. On any failure it calls error stop 1, so fpm test and CI see a nonzero exit code.

! test/test_solver.f90 -- regression + verification for the heat solver (Ch. 37-38).
! Uses the src/ modules; run with `fpm test`, or compile:
!   gfortran -std=f2018 -Wall -O2 src/kinds.f90 src/heat_types.f90 \
!            src/heat_solver.f90 test/test_solver.f90 -o test_solver && ./test_solver
! NEVER run during authoring; every expected value is hand-computed from Chapter 24.
program test_solver
  use kinds,       only: dp
  use heat_types,  only: field_t
  use heat_solver, only: laplacian, step, stable_dt
  implicit none

  real(dp), parameter :: pi = 3.141592653589793_dp

  ! The golden reference: the exact 5x5 field after two steps (Chapter 24), row-major.
  real(dp), parameter :: golden(5,5) = reshape([ &
      100.0_dp, 100.0_dp, 100.0_dp, 100.0_dp, 100.0_dp, &
        0.0_dp,  28.0_dp,  32.0_dp,  28.0_dp,   0.0_dp, &
        0.0_dp,   4.0_dp,   4.0_dp,   4.0_dp,   0.0_dp, &
        0.0_dp,   0.0_dp,   0.0_dp,   0.0_dp,   0.0_dp, &
        0.0_dp,   0.0_dp,   0.0_dp,   0.0_dp,   0.0_dp], [5,5], order=[2,1])

  type(field_t) :: f, s, mode
  real(dp) :: u(5,5), lap(5,5), dt, maxdev, residual, numeric, exact, x, y
  integer  :: i, j, n_fail

  n_fail = 0

  ! --- UNIT 1: laplacian exact on the quadratic x^2 + y^2 (== 4 everywhere).
  do j = 1, 5
    do i = 1, 5
      u(i,j) = (real(i-1,dp)*0.5_dp)**2 + (real(j-1,dp)*0.5_dp)**2
    end do
  end do
  lap = laplacian(u, 0.5_dp, 0.5_dp)
  call check('unit: laplacian(x^2+y^2) == 4', all(abs(lap(2:4,2:4) - 4.0_dp) < 1.0e-12_dp), n_fail)

  ! --- UNIT 2: laplacian exact on a linear field (== 0).
  do j = 1, 5
    do i = 1, 5
      u(i,j) = 25.0_dp * real(j-1, dp)
    end do
  end do
  lap = laplacian(u, 1.0_dp, 1.0_dp)
  call check('unit: laplacian(linear)  == 0', all(abs(lap(2:4,2:4)) < 1.0e-12_dp), n_fail)

  ! --- UNIT 3: stable_dt returns a CFL-safe step, r = alpha*dt/h^2 <= 1/4.
  dt = stable_dt(1.0_dp, 1.0_dp, 1.0_dp, safety=0.8_dp)            ! = 0.2
  call check('unit: stable_dt gives r <= 1/4', (1.0_dp*dt/1.0_dp**2) <= 0.25_dp, n_fail)

  ! --- REGRESSION: the 5x5 two-step field vs the golden reference (tolerance).
  call f%init(nx=5, ny=5, dx=1.0_dp, dy=1.0_dp)
  f%u(1,:) = 100.0_dp
  call step(f, 1.0_dp, 0.2_dp)
  call step(f, 1.0_dp, 0.2_dp)
  maxdev = maxval(abs(f%u - golden))
  call check('regression: 5x5 two-step matches golden (28/32/28, 4/4/4)', maxdev < 1.0e-9_dp, n_fail)

  ! --- VERIFY 1: the linear steady state is an exact fixed point of step.
  call s%init(nx=5, ny=5, dx=1.0_dp, dy=1.0_dp)
  do j = 1, 5
    do i = 1, 5
      s%u(i,j) = 25.0_dp * real(j-1, dp)          ! exact analytical steady state
    end do
  end do
  residual = s%u(3,3)
  call step(s, 1.0_dp, 0.2_dp)                     ! stepping must not change it
  residual = abs(s%u(3,3) - residual)
  call check('verify: linear steady state is a fixed point', residual < 1.0e-12_dp, n_fail)

  ! --- VERIFY 2: the maximum principle -- 0 <= u <= 100 (reuse the hot-top run f).
  call check('verify: maximum principle 0 <= u <= 100', &
             all(f%u >= 0.0_dp) .and. all(f%u <= 100.0_dp), n_fail)

  ! --- ANALYTICAL: one step of a sine mode vs its exact amplitude decay (reported).
  call mode%init(nx=5, ny=5, dx=0.25_dp, dy=0.25_dp)
  do j = 1, mode%ny
    do i = 1, mode%nx
      x = real(i-1, dp) * mode%dx
      y = real(j-1, dp) * mode%dy
      mode%u(i,j) = sin(pi*x) * sin(pi*y)         ! discrete fundamental mode
    end do
  end do
  dt = stable_dt(1.0_dp, mode%dx, mode%dy, safety=0.8_dp)         ! = 0.0125, r = 0.2
  call step(mode, 1.0_dp, dt)
  numeric = mode%u(3,3)                            ! center, x = y = 0.5
  exact   = exp(-2.0_dp*pi**2 * dt)                ! amplitude decay of the mode
  print '(a)',       'analytical: sine-mode center after 1 step'
  print '(a,f10.6)', '  numerical = ', numeric
  print '(a,f10.6)', '  exact     = ', exact
  print '(a,f10.6)', '  |error|   = ', abs(numeric - exact)

  print '(a,i0,a)', '--- ', 6 - n_fail, ' / 6 checks passed'
  if (n_fail > 0) error stop 1                     ! nonzero exit -> fpm test / CI go red

contains

  subroutine check(name, cond, nf)
    character(*), intent(in)    :: name
    logical,      intent(in)    :: cond
    integer,      intent(inout) :: nf
    if (cond) then
      print '(a,a)', 'PASS  ', name
    else
      nf = nf + 1
      print '(a,a)', 'FAIL  ', name
    end if
  end subroutine check
end program test_solver

! Hand-computed expectations (all from Chapter 24 physics; never run):
!   laplacian(x^2+y^2) = 4 at every interior point (exact for quadratics);
!   laplacian(linear ramp) = 0 (second difference of a line is 0);
!   stable_dt = 0.8/4 = 0.2 -> r = 0.2 <= 0.25;
!   two steps of the hot-top plate give the golden field, maxdev = 0;
!   the linear ramp is unchanged by step (residual 0); every value lies in [0,100];
!   sine mode: numeric = 0.765685, exact = exp(-0.2467401) = 0.781344, |error| = 0.015658.
!
! Expected output:
! PASS  unit: laplacian(x^2+y^2) == 4
! PASS  unit: laplacian(linear)  == 0
! PASS  unit: stable_dt gives r <= 1/4
! PASS  regression: 5x5 two-step matches golden (28/32/28, 4/4/4)
! PASS  verify: linear steady state is a fixed point
! PASS  verify: maximum principle 0 <= u <= 100
! analytical: sine-mode center after 1 step
!   numerical =   0.765685
!   exact     =   0.781344
!   |error|   =   0.015658
! --- 6 / 6 checks passed

I.8 The distributed and GPU variants

Three chapters take the same step interface onto larger machines by rewriting only its body and the data it exchanges. They are not reproduced here in full — each is a chapter-length treatment — but the complete, hand-verified programs live in those chapters' Project Checkpoints, and every one reproduces the identical field on any process, image, or thread count:

• Coarrays — Chapter 32. The plate is cut into vertical strips, one per image; each image runs the Chapter 24 stencil on its strip and exchanges its shared edge columns as halos with coindexed access u(:,k)[q]. Splitting along columns makes each halo a contiguous slice in column-major memory. Build with -fcoarray=single (one image) or OpenCoarrays caf/cafrun (many).
• MPI — Chapter 34. One-dimensional domain decomposition into horizontal strips with ghost rows; each step exchanges halos with MPI_Sendrecv (deadlock-free), then runs the ordinary stencil over the owned rows. Build with mpif90 and run with mpirun -np N.
• GPU — Chapter 35 (optional, advanced). A single OpenACC !$acc data region wraps the entire time loop, so the field is copied to the device once and swept there every step with !$acc parallel loop collapse(2) — no per-step host↔device transfer. Build with nvfortran -acc (or gfortran -fopenacc).

That is the whole solver: one small program, grown across forty chapters from a hard-coded array update into a modular, validated, timed, visualized, and parallel simulation — and readable, at every stage, by someone who has never seen it before. That readability is not decoration; it is what makes a computational result something another scientist can trust.