Case Study 1: From Six Arrays to One Type

"The parallel arrays compiled fine for ten years. Then someone inserted a particle in the middle of five of the six, and the physics quietly broke."

Executive Summary

You have inherited a small N-body kernel — the kind of code that simulates gravitating bodies, or atoms in a molecular-dynamics run — and it stores its particles the way pre-1990 Fortran forced everyone to: as a fistful of parallel arrays, one array per physical quantity. It works, but it is a standing invitation to a whole family of bugs, and it obscures the physics. In this case study you read that code, diagnose precisely what makes it fragile, and refactor it into a single body derived type with type-bound-ready structure. Then — because a computational scientist never refactors without asking about performance — you confront the one honest cost of the change: the Array-of-Structures versus Structure-of-Arrays trade, which decides whether your new type helps or hurts the hot loop. By the end you can look at any bag-of-parallel-arrays code and refactor it into types with your eyes open to the layout consequences.

Skills applied: defining derived types and nested types (§9.1); passing a type as one argument instead of many (§9.1); the array-of-derived-types pattern and component sections (§9.1, §9.5); reasoning about memory layout and cache behaviour (§9.5, foreshadowing Part VII).

Background

The code models $N$ point masses. For each body it must know a position $(x, y, z)$, a velocity, and a mass. The inherited version — call it the "before" — declares six separate arrays:

real(dp) :: px(N), py(N), pz(N)     ! positions
real(dp) :: vx(N), vy(N), vz(N)     ! velocities
real(dp) :: mass(N)                 ! masses

and every routine that touches a particle takes seven arguments, in an order everyone is expected to remember. It has run correctly for a decade. Our job is not to prove it wrong but to make it safe to change and clear to read — and to do so without accidentally making it slower.

Phase 1 — Read the Code and Name the Hazards

Three concrete hazards live in the parallel-array design, and naming them is half the work.

Hazard 1 — the arrays can fall out of step. Nothing ties px(i), vy(i), and mass(i) together except the shared index i and the programmer's discipline. Insert a body at position 500 by shifting px, py, pz, vx, vy but forget vz (the epigraph's bug), and every body past 500 now has the wrong $z$-velocity — a corruption the compiler cannot see, because as far as it knows these are seven unrelated arrays.

Hazard 2 — the argument lists are long and orderable-wrong. A routine subroutine kick(px, py, pz, vx, vy, vz, mass, dt) has eight arguments of nearly identical type. Transpose two real(dp) arrays at a call site and it still compiles and runs — and computes nonsense. intent cannot save you here; the types match even when the meanings do not.

Hazard 3 — the physics is buried. 0.5_dp * mass(i) * (vx(i)**2 + vy(i)**2 + vz(i)**2) is kinetic energy, but you have to decode it. There is no body in the code, only a swarm of subscripts.

The diagnosis: every one of these hazards is the same root cause — data that belongs together is not together. The seven arrays describe one concept, "a body," that the language has never been told about. A derived type is exactly the act of telling it.

Phase 2 — Refactor to a body Type

We introduce a vec3 for the two 3-vectors and a body that nests them with the mass. Then the seven arrays collapse into one array of body, and the seven-argument routines become one-argument routines that take a whole body (or an array of them). Here is the refactored kernel, complete and compilable, computing two quantities every N-body code needs — the total mass, the center of mass, and the total momentum:

module nbody_types
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  private
  public :: vec3, body, center_of_mass, total_momentum

  type :: vec3
    real(dp) :: x = 0.0_dp, y = 0.0_dp, z = 0.0_dp
  end type vec3

  type :: body
    type(vec3) :: pos          ! position  — nested type
    type(vec3) :: vel          ! velocity  — nested type
    real(dp)   :: mass = 0.0_dp
  end type body

contains

  pure function center_of_mass(bodies) result(c)
    type(body), intent(in) :: bodies(:)
    type(vec3) :: c
    real(dp) :: m
    integer  :: i
    m = sum(bodies%mass)                       ! component section: all masses
    c = vec3(0.0_dp, 0.0_dp, 0.0_dp)
    do i = 1, size(bodies)
      c%x = c%x + bodies(i)%mass * bodies(i)%pos%x
      c%y = c%y + bodies(i)%mass * bodies(i)%pos%y
      c%z = c%z + bodies(i)%mass * bodies(i)%pos%z
    end do
    c%x = c%x / m;  c%y = c%y / m;  c%z = c%z / m
  end function center_of_mass

  pure function total_momentum(bodies) result(p)
    type(body), intent(in) :: bodies(:)
    type(vec3) :: p
    integer :: i
    p = vec3(0.0_dp, 0.0_dp, 0.0_dp)
    do i = 1, size(bodies)
      p%x = p%x + bodies(i)%mass * bodies(i)%vel%x
      p%y = p%y + bodies(i)%mass * bodies(i)%vel%y
      p%z = p%z + bodies(i)%mass * bodies(i)%vel%z
    end do
  end function total_momentum

end module nbody_types

program nbody_demo
  use nbody_types
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  type(body) :: swarm(3)
  type(vec3) :: com, mom

  swarm(1) = body(pos=vec3(0.0_dp,0.0_dp,0.0_dp), vel=vec3(1.0_dp,0.0_dp,0.0_dp), mass=1.0_dp)
  swarm(2) = body(pos=vec3(1.0_dp,0.0_dp,0.0_dp), vel=vec3(0.0_dp,1.0_dp,0.0_dp), mass=2.0_dp)
  swarm(3) = body(pos=vec3(0.0_dp,1.0_dp,0.0_dp), vel=vec3(0.0_dp,0.0_dp,1.0_dp), mass=3.0_dp)

  com = center_of_mass(swarm)
  mom = total_momentum(swarm)

  print '(a, f6.2)',  'total mass    = ', sum(swarm%mass)
  print '(a, 3f8.3)', 'center of mass= ', com%x, com%y, com%z
  print '(a, 3f8.3)', 'momentum      = ', mom%x, mom%y, mom%z
end program nbody_demo

$ gfortran -std=f2018 -Wall -O2 nbody.f90 -o nbody && ./nbody
total mass    =   6.00
center of mass=    0.333   0.500   0.000
momentum      =    1.000   2.000   3.000

Sanity check by hand. Total mass is $1 + 2 + 3 = 6$. The center of mass is $\frac{\sum_i m_i \mathbf{r}_i}{\sum_i m_i} = \frac{1(0,0,0) + 2(1,0,0) + 3(0,1,0)}{6} = \frac{(2,3,0)}{6} = (0.333, 0.5, 0)$. The total momentum is $\sum_i m_i \mathbf{v}_i = 1(1,0,0) + 2(0,1,0) + 3(0,0,1) = (1,2,3)$. Every printed number matches. Note what the refactor bought: center_of_mass(swarm) takes one argument, a body cannot be assembled with its velocity where its position belongs, and bodies(i)%vel%z reads like the physics it is.

Phase 3 — The Honest Cost: Array-of-Structures vs Structure-of-Arrays

Here is where a careless refactor can betray you. Our body type lays each particle's fields out contiguously — pos.x, pos.y, pos.z, vel.x, vel.y, vel.z, mass, seven real(dp)s, then the next body. An array of body is therefore an Array-of-Structures (AoS): the data for one particle is together, but the same field across particles is scattered. Consecutive masses are $7 \times 8 = 56$ bytes apart.

The original parallel arrays were the opposite — a Structure-of-Arrays (SoA) — where all the masses sit in one contiguous mass(:) array, 8 bytes apart.

Why does this matter? A CPU fetches memory in cache lines, typically 64 bytes at a time. Consider a loop that needs only the masses (a common case — computing total mass, or scaling forces):

In the SoA layout the mass loop streams through memory using every byte it fetches; in the AoS layout it drags an entire particle into cache to use one number of it, wasting roughly seven-eighths of the memory bandwidth. For a bandwidth-bound loop over millions of particles, that can be several times slower — a real penalty you will measure in Part VII.

⚡ The rule of thumb: AoS is clearer and is faster when you touch most fields of one particle at a time (as a force calculation does — it needs every body's full state). SoA is faster when you sweep one field across all particles (as a whole-array update does). Neither wins universally; the access pattern decides.

Phase 4 — Resolve the Trade Deliberately

You are not forced to choose blind. Three defensible resolutions, in order of how often they are right:

1. Keep the AoS body type if the dominant loop is the force calculation (it usually is in N-body and MD), because that loop reads a body's whole state and AoS keeps it contiguous. Clarity and the common case agree.
2. Adopt an SoA type — one particles type whose components are arrays (real(dp), allocatable :: x(:), y(:), …) — if profiling shows the hot loops are single-field sweeps. You still get one named object and safe argument passing; the arrays just live inside it. This is Exercise 9.23(b).
3. Only refactor the interface, measure, then decide the layout — the professional path. The type is a contract; you can change its internal layout later without touching the call sites, precisely because the data is now behind one name.

The key insight is that the derived type did not cause the layout question — it surfaced it. The parallel-array code had already committed to SoA silently; wrapping the data in a type just made the choice visible and changeable. That visibility is a feature.

Phase 5 — What We Ported, and What We Gained

Line for line, the refactor was modest: seven arrays became two types, and seven-argument routines became one-argument routines. But the three hazards of Phase 1 are gone. The arrays can no longer fall out of step, because there are no longer seven arrays — there is one array of one type, and inserting a body moves all of its fields together or none. The argument-order bugs are gone, because there is one argument. And the physics is legible. We paid for this with a layout decision we then made deliberately rather than by default. That is the whole transaction: a derived type trades a silent, fragile data layout for an explicit, safe, and changeable one — at a memory-layout cost you now know how to reason about.

Discussion Questions

1. The epigraph's bug — shifting five of six arrays — is impossible in the AoS body design. Explain precisely why, in terms of what "one array of body" guarantees that "six parallel arrays" does not.
2. A colleague says "derived types are slower, so real HPC codes avoid them." Using the AoS/SoA analysis, explain why this is half-true and half-myth, and what actually determines the speed.
3. You must add a charge field to every particle. Compare the edit required in the parallel-array version with the edit required in the body-type version. Which is safer, and why?

Your Turn: Extensions

• Option A. Take the SoA path: write a particles type whose components are allocatable arrays (x(:), y(:), z(:), vx(:), …, mass(:)), give it a type-bound n() returning the count, and rewrite total_momentum to sweep the component arrays. Which version reads more clearly? Which would you expect to vectorize better?
• Option B. Add a type-bound kinetic_energy(self) to body returning $\tfrac{1}{2}m|\mathbf{v}|^2$, and a module function total_kinetic_energy(bodies). Verify by hand for the three-body swarm above (you should get $\tfrac{1}{2}(1\cdot 1 + 2\cdot 1 + 3\cdot 1) = 3.0$).
• Option C. Instrument both layouts (AoS and SoA) with system_clock (Chapter 28) on a ten-million-particle mass-sum loop and measure the ratio. Does it match the seven-eighths-wasted-bandwidth prediction?

Key Takeaways

• Parallel arrays are silent AoS-or-SoA decisions dressed up as "no data structure at all"; they are fragile because the language cannot see that the arrays belong together.
• A derived type makes the concept explicit, collapses long argument lists to one, and makes the physics legible — and it makes bugs of the "arrays fell out of step" family unwriteable.
• The one honest cost is memory layout: an array of a type is Array-of-Structures, which is fast for whole-particle loops and slow for single-field sweeps. Choose the layout from the access pattern, and hide it behind the type so you can change your mind after profiling.