Case Study 2: Designing step() for the Long Haul

"An interface is a promise. The cost of breaking it is paid by everyone who believed you."

Executive Summary

Where the first case study took apart existing code, this one asks you to design a piece that does not exist yet — and to design it knowing it must survive thirty more chapters of change. The heat solver's update procedure, step, will be reimplemented many times before this book ends: it gets a real finite-difference stencil in Chapter 24, moves into a module in Chapter 8, takes a field_t derived type in Chapter 9, and sprouts OpenMP directives in Chapter 33. Through all of that, its interface — the signature the rest of the code calls — should barely change. This study is about making the decisions now, with the tools of Chapter 6, that keep that promise: the right argument list, the right intents, an optional argument used judiciously, and a pure helper factored out so the compiler and the parallelizer can both do their jobs.

Skills applied: choosing a subroutine over a function (§6.1); assigning intent to every argument including the array (§6.2); adding an optional argument with a present guard (§6.3); factoring a pure kernel helper (§6.4); using an internal procedure with an explicit interface (§6.5); designing around assumed-shape so one routine fits every grid size (§6.6).

Background

You are designing the procedure that advances the temperature field by one time step. You know four things about its future, and each constrains the design today:

1. The real physics is not ready — the genuine stencil, the CFL-stable timestep, and the boundary conditions are Chapter 24's job. So the body must be replaceable without disturbing callers.
2. It must eventually run on grids from $4\times4$ (a unit test) to $10^4\times10^4$ (a production run), so it cannot bake in a size.
3. Some runs want the solver to enforce fixed-temperature edges every step; others manage boundaries elsewhere. The routine must serve both without two separate procedures.
4. It will be parallelized, so the per-cell update must be expressible as an independent, side-effect-free computation the compiler is free to vectorize and spread across cores.

The naive first cut — inline the update in the main program with hard-coded dimensions — fails every one of these. We will design against them instead.

Phase 1 — Fix the Signature Before the Body

The most durable decision is the shape of the call, so make it first and deliberately. The update acts on the field in place and returns nothing to assign, so it is a subroutine, not a function (§6.1). Its arguments and their intents:

Two design commitments are worth defending. First, field is assumed-shape (field(:,:)), so the identical routine serves the $4\times4$ test and the $10^4\times10^4$ run — the size travels with the array and is never passed by hand (§6.6). Second, bc_value is optional (§6.3): a caller who wants the solver to hold the edges passes it; a caller who manages boundaries elsewhere omits it and pays nothing. Making it mandatory would force every call site to supply an edge value even when it is meaningless; leaving it out entirely would rule out the common "hot edge" case. Optional is exactly the seam that serves both.

Phase 2 — Factor Out a pure Kernel

Requirement 4 — future parallelism — is won or lost here. If the per-cell update is a tangle of reads and writes, the compiler cannot know the cells are independent. So we isolate the arithmetic for one cell into a pure function that reads the neighbourhood and returns a number, touching nothing else:

pure function five_point(u, i, j) result(lap)
  real(dp), intent(in) :: u(:,:)      ! assumed-shape; read-only
  integer,  intent(in) :: i, j
  real(dp) :: lap
  lap = u(i-1, j) + u(i+1, j) + u(i, j-1) + u(i, j+1) - 4.0_dp * u(i, j)
end function five_point

Because five_point is pure, the compiler knows a call has no side effects: evaluating it for cell $(i,j)$ cannot disturb the evaluation for any other cell. That is the precise property a do concurrent loop or an OpenMP parallel do needs to spread the update across cores safely — and we have secured it now, in Chapter 6, simply by writing one honest keyword (§6.4). The helper is also where the real stencil will slot in unchanged at Chapter 24: same signature, better arithmetic.

Phase 3 — Build the Procedure

Assemble the pieces. The edges are set first (only if bc_value was supplied), then every interior cell is updated from a single snapshot old, so the sweep order cannot bias the result:

program step_design
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp), allocatable :: field(:,:)
  real(dp), parameter   :: alpha = 0.1_dp, dt = 1.0_dp

  allocate(field(4, 4))
  field = 0.0_dp
  call step(field, alpha, dt, bc_value=100.0_dp)   ! set edges to 100, then update
  call show(field)

contains

  subroutine step(field, alpha, dt, bc_value)
    real(dp), intent(inout)        :: field(:,:)
    real(dp), intent(in)           :: alpha, dt
    real(dp), intent(in), optional :: bc_value
    real(dp), allocatable :: old(:,:)
    integer :: i, j, nx, ny

    nx = size(field, 1)
    ny = size(field, 2)

    if (present(bc_value)) then           ! optional boundary enforcement
       field(1,  :) = bc_value
       field(nx, :) = bc_value
       field(:,  1) = bc_value
       field(:, ny) = bc_value
    end if

    allocate(old(nx, ny))
    old = field                           ! snapshot: interior updates are order-free
    do j = 2, ny - 1
       do i = 2, nx - 1
          field(i, j) = old(i, j) + alpha * dt * five_point(old, i, j)
       end do
    end do
  end subroutine step

  pure function five_point(u, i, j) result(lap)
    real(dp), intent(in) :: u(:,:)
    integer,  intent(in) :: i, j
    real(dp) :: lap
    lap = u(i-1, j) + u(i+1, j) + u(i, j-1) + u(i, j+1) - 4.0_dp * u(i, j)
  end function five_point

  subroutine show(f)
    real(dp), intent(in) :: f(:,:)
    integer :: i
    do i = 1, size(f, 1)
       print '(4f8.2)', f(i, :)
    end do
  end subroutine show

end program step_design

Phase 4 — Exercise It, and Check by Hand

Run the design on a $4\times4$ plate that starts cold, with the edges pinned to 100:

$ gfortran -std=f2018 -Wall -O2 step_design.f90 -o step_design && ./step_design
  100.00  100.00  100.00  100.00
  100.00   20.00   20.00  100.00
  100.00   20.00   20.00  100.00
  100.00  100.00  100.00  100.00

Confirm it by hand, cell by cell, because a design you cannot hand-check is a design you do not understand. After the boundary step, the field is 100 around the entire border and 0 in the four interior cells. The snapshot old therefore holds 100 on the edges and 0 inside. For interior cell $(2,2)$, whose four neighbours in old are the edge cell $(1,2)=100$ above, the interior cell $(3,2)=0$ below, the edge cell $(2,1)=100$ to the left, and the interior cell $(2,3)=0$ to the right:

$$ \text{five\_point} = 100 + 0 + 100 + 0 - 4\times 0 = 200, \qquad \text{field}(2,2) = 0 + \alpha\, dt \times 200 = 0 + 0.1 \times 200 = 20. $$

By the symmetry of the setup all four interior cells see two hot edges and two cold interiors, so every one lands at exactly 20.00 — which is what printed. The physics is a placeholder, but the interface behaves, and that is what we set out to verify.

The reasoning that matters: every requirement from the Background is now met by a Chapter 6 feature. Size-independence came from assumed-shape (§6.6); the dual boundary behavior from an optional argument (§6.3); the parallel-readiness from a pure helper (§6.4); the "replace the body, keep the callers" property from committing to the signature and intents first (§6.1–6.2). None of it required anything from later chapters — the design tools were all in this one.

Phase 5 — Future-Proofing: What the Interface Must Survive

Walk the interface forward through the book and confirm it holds:

• Chapter 24 (real stencil). Only five_point's body changes — a proper Laplacian with dx, dy, and a CFL-limited dt. The signature of step is untouched, so every caller is unaffected.
• Chapter 8 (modules). step and five_point move verbatim into a heat_solver module. They already have explicit interfaces from being internal procedures (§6.5), so the move is mechanical and the calls do not change.
• Chapter 9 (derived type). field(:,:) may become a field_t carrying nx, ny, dx, dy, u(:,:). This is a signature change — the one we tolerate — but the intents and the optional bc_value survive intact.
• Chapter 33 (OpenMP). The interior loop gets a !$omp parallel do (or becomes do concurrent). Because five_point is pure and the update reads only old, the loop is already correct to parallelize — the design paid for this in Phase 2.

The lesson generalizes past this one procedure: design the call, the intents, and the purity first; treat the body as replaceable. An interface you chose deliberately is cheap to keep; one you fell into is expensive to change, because the cost lands on every line that ever called it.

Discussion Questions

1. We made bc_value a single scalar, forcing all four edges to the same temperature. A real plate might want a different value per edge. Sketch two interface options (four optional scalars; one small array of edge values) and argue which better honors "don't break the interface later."
2. step allocates old on every call. Name two ways to avoid that (an argument the caller owns; a module-level workspace) and state the cost each imposes on the interface versus on performance (preview of Chapter 28).
3. five_point takes (u, i, j) and is called once per interior cell. An alternative is a pure function laplacian(u) returning the whole interior at once as an array. Which is friendlier to the Chapter 33 parallelization, and which is friendlier to reading? Do they have to be the same choice?

Your Turn: Extensions

• Option A. Add an optional, intent(out) :: max_change argument that reports the largest $\lvert \text{field}_{\text{new}} - \text{old}\rvert$ over the interior, computed only when the caller asks. This is the convergence check a steady-state solver needs — and a textbook use of an optional output.
• Option B. Write a second driver that calls step without bc_value on a field whose edges you set once before a loop of many steps. Confirm the interior evolves and the edges stay put, demonstrating that the optional argument genuinely serves two workflows.
• Option C. Replace the (i,j) helper with a whole-interior pure function laplacian(u) returning an array, and rewrite step's loop as a single whole-array update. Compare the two designs for readability and for how obviously parallel each is.

Key Takeaways

• Design the interface first — argument list, intents, purity — and treat the procedure body as a replaceable detail. The signature is the promise; the body is an implementation.
• An optional argument is the right tool when one procedure must serve two workflows (enforce boundaries or not) without splitting into two; guard it with present and it costs absent callers nothing.
• Factoring the per-cell arithmetic into a pure helper is what makes future vectorization and parallelism safe by construction — you buy Chapter 33's parallelism with a keyword written in Chapter 6.
• Assumed-shape makes one procedure fit every grid size, so the same step you hand-check on a $4\times4$ is the one that runs on the cluster.