Case Study 2: Building an Extensible Solver Framework

"Every select type that decides behavior is a method that has not yet been written."

Executive Summary

Where the first case study read a framework, this one asks you to build one — and to fix a common mistake along the way. You are given a solver driver that chooses between several stepping schemes with a large select type block, and asked to make it extensible: a colleague wants to add a scheme without editing the driver. You will diagnose why the select type design resists that, refactor it into an abstract solver_t with a deferred step, and add a second scheme to prove the new design accepts extensions cleanly. Then you will place the polymorphism deliberately — coarse-grained, one dispatch per timestep — so the abstraction costs nothing where it matters. The result is the optional Project Checkpoint framework, and the reasoning that produced it.

Skills applied: recognizing the select type design smell (§10.3); designing an abstract type and a deferred contract (§10.4); type extension to add an implementation (§10.1); dispatch through a class handle (§10.2); and the coarse-grain performance discipline (§10.5, §10.6).

Background

The inherited driver advances a temperature field and supports two schemes: an explicit diffusion step and a relaxation (averaging) step. The previous author selected between them with the object's type:

! The design we are replacing: behavior chosen by a select type in the driver.
subroutine step_dispatch(scheme, fld)
  class(*),      intent(in)    :: scheme
  type(field_t), intent(inout) :: fld
  select type (scheme)
  type is (heat_scheme_t)
    ! ... explicit diffusion update inline ...
  type is (relax_scheme_t)
    ! ... averaging update inline ...
  class default
    error stop 'unknown scheme'
  end select
end subroutine

It works, and it is exactly wrong for the goal. Every new scheme means a new type is arm inside this subroutine — and inside every other subroutine that switches on the scheme. The variation is scattered across the driver instead of living with each scheme. That is the smell §10.3 warned about: a select type that implements behavior is inheritance turned inside out.

Phase 1 — Design the Contract

The fix is to invert the dependency. Instead of the driver knowing every scheme, each scheme knows how to step itself, and the driver knows only the interface. That interface is an abstract type with one deferred binding:

module solver_base
  use kinds,      only: dp
  implicit none
  private
  public :: solver_t

  type, abstract :: solver_t
    real(dp) :: alpha = 0.0_dp
    real(dp) :: dt    = 0.0_dp
  contains
    procedure(step_i), deferred :: step
  end type solver_t

  abstract interface
    subroutine step_i(self, fld)
      use heat_types, only: field_t
      import :: solver_t
      class(solver_t), intent(in)    :: self
      type(field_t),   intent(inout) :: fld
    end subroutine step_i
  end interface
end module solver_base

The contract is now a single sentence the compiler enforces: a solver is anything that can step a field. Nothing in solver_base mentions heat, or relaxation, or any scheme. That ignorance is the design's strength — the base cannot be broken by a new scheme, because it does not know schemes exist.

Phase 2 — Build Two Implementations

Each scheme becomes an extends(solver_t) type that supplies its own step. The explicit diffusion scheme carries a diffusivity and a timestep; the relaxation scheme carries an over-relaxation factor (with $\omega = 1$ giving plain neighbour-averaging):

module schemes
  use kinds,       only: dp
  use heat_types,  only: field_t
  use solver_base, only: solver_t
  implicit none
  private
  public :: heat_solver_t, relax_solver_t

  type, extends(solver_t) :: heat_solver_t
  contains
    procedure :: step => heat_step
  end type heat_solver_t

  type, extends(solver_t) :: relax_solver_t
    real(dp) :: omega = 1.0_dp                 ! 1.0 = plain Jacobi averaging
  contains
    procedure :: step => relax_step
  end type relax_solver_t

contains

  subroutine heat_step(self, fld)
    class(heat_solver_t), intent(in)    :: self
    type(field_t),        intent(inout) :: fld
    real(dp), allocatable :: u_new(:,:)
    real(dp) :: r
    integer  :: i, j
    r     = self%alpha * self%dt / fld%dx**2
    u_new = fld%u
    do j = 2, fld%ny - 1
      do i = 2, fld%nx - 1
        u_new(i,j) = fld%u(i,j) + r * ( fld%u(i-1,j) + fld%u(i+1,j) &
                   + fld%u(i,j-1) + fld%u(i,j+1) - 4.0_dp*fld%u(i,j) )
      end do
    end do
    fld%u = u_new
  end subroutine heat_step

  subroutine relax_step(self, fld)
    class(relax_solver_t), intent(in)    :: self
    type(field_t),         intent(inout) :: fld
    real(dp), allocatable :: u_new(:,:)
    real(dp) :: avg
    integer  :: i, j
    u_new = fld%u
    do j = 2, fld%ny - 1
      do i = 2, fld%nx - 1
        avg = 0.25_dp * ( fld%u(i-1,j) + fld%u(i+1,j) + fld%u(i,j-1) + fld%u(i,j+1) )
        u_new(i,j) = (1.0_dp - self%omega)*fld%u(i,j) + self%omega*avg
      end do
    end do
    fld%u = u_new
  end subroutine relax_step

end module schemes

Notice what the refactor accomplished. The two update formulas that were tangled inside one select type now live in two separate, self-contained procedures, each bound to its own type. A third scheme would be a third module type with its own step — and, crucially, the driver below would not change a character.

Phase 3 — Drive It Polymorphically

The driver holds a class(solver_t) and steps it. It is written once, against the abstract type, and works for every scheme. Here it runs the diffusion scheme for three steps and the relaxation scheme for one, on identical $3\times3$ plates with the top edge held at $100$:

program cs02
  use kinds,       only: dp
  use heat_types,  only: field_t
  use solver_base, only: solver_t
  use schemes,     only: heat_solver_t, relax_solver_t
  implicit none
  class(solver_t), allocatable :: sim
  type(field_t) :: plate
  integer :: n

  ! --- explicit diffusion, three steps ---
  call setup(plate)
  allocate(heat_solver_t :: sim)
  sim%alpha = 1.0_dp
  sim%dt    = 0.1_dp
  do n = 1, 3
    call sim%step(plate)
  end do
  print '(a, f8.3)', 'heat  solver, interior after 3 steps = ', plate%u(2,2)
  deallocate(sim)

  ! --- relaxation, one step: same interface, one line changed ---
  call setup(plate)
  allocate(relax_solver_t :: sim)              ! the ONLY line that differs
  call sim%step(plate)
  print '(a, f8.3)', 'relax solver, interior after 1 step  = ', plate%u(2,2)

contains

  subroutine setup(f)
    type(field_t), intent(out) :: f
    f%nx = 3; f%ny = 3; f%dx = 1.0_dp; f%dy = 1.0_dp
    if (allocated(f%u)) deallocate(f%u)
    allocate(f%u(f%nx, f%ny))
    f%u       = 0.0_dp
    f%u(1, :) = 100.0_dp
  end subroutine setup
end program cs02

$ gfortran -std=f2018 -Wall kinds.f90 heat_types.f90 solver_base.f90 schemes.f90 cs02.f90 -o cs02 && ./cs02
heat  solver, interior after 3 steps =   19.600
relax solver, interior after 1 step  =   25.000

The two numbers are hand-computed and physically telling. With $r = \alpha\Delta t/\Delta x^2 = 0.1$, the diffusion scheme creeps the interior point up: $0 \to 10 \to 16 \to 19.6$. The relaxation scheme, with $\omega=1$, sets the interior point to the average of its fixed neighbours in a single step: $0.25\,(100+0+0+0) = 25$ — which is exactly the steady state the diffusion scheme is slowly crawling toward. Two solvers, one interface, the same physics answered at different speeds.

Phase 4 — Place the Polymorphism

The design is extensible; now make sure it is also fast, which means placing the dispatch deliberately. Look at where class appears and where it does not. The class(solver_t) handle — the polymorphism — lives in the driver, and its step is dispatched once per timestep. That is coarse-grained: one indirect call amortized over the entire grid update, utterly negligible even for millions of steps.

Now look inside heat_step and relax_step. There is no class anywhere. fld%u is a plain, contiguous, monomorphic 2D array; the inner loop runs over the first index, the fast direction for column-major storage (Chapter 5); nothing is dispatched per cell. The optimizer is free to inline, unroll, and vectorize the kernel exactly as if no object orientation existed — because, at the level that matters for speed, none does.

The reasoning that matters: this is §10.5 made concrete. We abstracted the solver (a decision made once) and kept the cells flat (arithmetic done a billion times). Had we instead made each cell a polymorphic object and dispatched a method per cell, the same physics would run with an un-inlinable indirect call in the hottest loop — extensible and slow. The framework's shape puts the indirection where it is free and forbids it where it would be ruinous.

Phase 5 — What You Built, and What It Costs to Extend

Step back and tally the design's properties against the goal. Extensibility: adding a scheme is one new extends(solver_t) type with a step; the driver and the other schemes are untouched, and the compiler guarantees the newcomer conforms. Correctness: the deferred contract makes "a solver that forgets to step" a compile error, not a run-time surprise. Performance: the one dispatch per step is free, and the kernels are as fast as hand-written non-OOP code. The cost: a little more ceremony than a bare subroutine — an abstract type, an interface block, the import — justified only because there is genuine plural variation (more than one scheme) to manage. If there were only ever one scheme, none of this would be worth it, and §10.5 would tell you to keep the plain module from Chapter 8. The framework earns its complexity precisely because "one of several solvers" is a real requirement.

Discussion Questions

1. The select type version and the abstract-type version compute identical results. Articulate the difference in terms of who must change when a fourth scheme is added. Why is "the driver does not change" the decisive property?
2. relax_solver_t carries an omega component that heat_solver_t does not. How does the framework accommodate schemes with different parameters behind one step interface?
3. The dispatch is once per timestep here. Describe a redesign that would (wrongly) push it to once per cell, and estimate — in words — the damage, referencing §10.5.

Your Turn: Extensions

• Option A. Add a third scheme, source_solver_t, whose step adds a constant source term self%source * self%dt to every interior cell. Confirm you changed only one new module type and one allocate line in the driver. Hand-compute three steps from a zero interior with source=50, dt=0.1.
• Option B. Give solver_t a second deferred binding, name(self), returning a short label, and have the driver print which scheme it is running. What must every existing scheme now provide, and how does the compiler make sure you did not forget one?
• Option C. Compose a second axis of variation: an abstract boundary_t with a deferred apply(self, fld), held as a component of solver_t, so boundary treatment varies independently of the stepping scheme. Sketch the types. Why is composition (a solver has a boundary) better here than trying to express both axes through inheritance alone?

Key Takeaways

• A select type that dispatches behavior scatters variation across the code; an abstract type with a deferred binding gathers each behavior with its type, so adding an implementation touches nothing else.
• The decisive test of an extensible design is who must change to add a case: here, only a new type — never the driver, never the existing schemes.
• Put the polymorphism at the coarse grain (choose a solver once, dispatch per timestep) and keep the kernels monomorphic; you get extensibility and full optimization at the same time.
• OOP earns its ceremony only when there is genuine plural variation to manage. One scheme wants a plain module; "one of several schemes" wants this framework.