Parallelizing the SWAT+ engine across CPU cores

Can the SWAT+ engine use the cores a server already has?

SWAT+ runs one daily time step at a time, walking every object in the watershed — HRUs, routing units, channels, reservoirs, aquifers — in a single sequential order. For a small basin that is fine. For regional, hyper-resolution models with tens of thousands of HRUs it is the wall: a one-year run of the 57,998-HRU Peace River benchmark takes ~3 minutes on a fast modern core (and 5–6 on an older one), and a multi-decade calibration multiplies that by thousands of evaluations.

Modern servers have many idle cores. The question this page answers: can the SWAT+ engine itself use them — on one machine, shared memory, without changing the science — by running independent parts of the watershed at the same time?

A full-DAG wavefront over the routing network

Within one day, two objects can run concurrently only if neither is downstream of the other. The routing network is a directed acyclic graph (DAG): headwater HRUs feed routing units, which feed channels, which feed the next channel down to the outlet. We compute, for every object, its longest path from a headwater leaf — its level. Objects that share a level are mutually independent and run together; level L+1 waits for level L.

We then drive the daily step level-by-level under an OpenMP parallel loop. The land phase (tens of thousands of HRUs) is one enormous wide level; the channel network narrows as it converges on the main stem. Making this safe required auditing every shared 'current-object' scratch variable in the engine and giving each thread its own copy (threadprivate), so concurrent objects never clobber one another.

How we proved it correct

We validate by thread-count invariance: independent runs at different thread counts (and 1-thread vs the stock serial engine) must produce identical output to the last digit. Any run-to-run difference is, by definition, a real race — and points straight at the unprotected shared state. This turns correctness into a precise, automatable test, and it is what later exposed the one remaining race (in the particulate constituents) that a streamflow-only check would have missed.

Table 2. Engine fork, benchmark model, and dedicated hosts used for the measurements on this page.

The optimization journey: from I/O to CPU time

Making a very large SWAT+ model fast was not one fix but a sequence — each step removed the bottleneck the previous step exposed, and the dominant cost kept moving: from disk, to startup, to the land phase, to channel routing, and now to thread synchronization. The colour shows the class of each fix; none of them change the model's science.

Why it unlocks scaling: PSO calibration is capped at ~40 vCPU per model — with ≤40 parameters and good initial simulations, a larger population doesn't help (and can hurt). A parallel engine is the only way to put more than 40 cores on a single model, so one model's calibration can scale out to the whole cluster.

All four upstream pull requests (#213, #214, #219, #220) are open against swat-model/swatplus; the engine runs them today in our fork.

Two products: a serial gain for everyone, and thread scaling for one model

Peace River HUC-8, 57,998 HRUs, 1 simulated year, dedicated AWS c8a.8xlarge (AMD EPYC "Turin", 32 cores, SMT off) · thread scaling of the final engine, self-relative to its own 1-thread wall; ifx -O3 -ipo; daily channel_sd output; thread-pinned, best of two on a dedicated box.

Figure 1. Measured thread scaling of the final engine on the dedicated c8a box (AMD EPYC "Turin", 32 cores, SMT off), 57,998-HRU benchmark, one simulated year, full wavefront, self-relative to its own 1-thread wall. Bars are measured speedup; the dashed line is ideal linear scaling. The full wavefront reaches 5.33× at 24 threads; with the 1.67× serial gain on top, the end-to-end acceleration versus the stock single-core engine is ~7.1×.

Table 1. Wall-clock time and self-relative speedup by thread count (final engine, dedicated AWS c8a.8xlarge, one simulated year, ifx -O3 -ipo, full wavefront, daily channel_sd output, best of two).

The campaign yields two independent products. Group A — serial overhead removal (the same fixes that drive the I/O→CPU journey above) — needs no threads and changes no model behavior, so it accelerates every existing run on a single core: a cumulative 1.67× on this c8a node (the stock engine's one-year wall drops 191.9 s → 115.2 s), and larger on more memory-starved cores (up to 2.47× on the Intel "Granite Rapids" part). Group B — the routing-DAG wavefront — then scales a single model across cores.

On the dedicated 32-core c8a box the full wavefront scales to 5.33× at 24 threads (table and chart, self-relative to its own 1-thread wall). Multiplying the two layers, the end-to-end acceleration versus the stock single-core engine is ~7.1× — a one-year run drops from 191.9 s to 26.9 s. Two channel-side fixes got the scaling there: removing the dominant ch_temp serial hotspot (which shrank the Amdahl serial fraction and so lifted the achievable ceiling), then fusing the width-1 main-stem wave levels.

Scaling is bounded by physical cores and memory bandwidth, not thread count — this is a memory-bound model. The 32-core Turin keeps climbing to 24 threads; the 16-core Intel and AMD parts peak at their physical-core count and regress when threads spill onto SMT siblings. The full cross-hardware study on three dedicated cloud CPUs is below.

Cross-hardware scaling on dedicated cloud CPUs

The identical engine and the same one-simulated-year benchmark workload, thread-pinned (OMP_PROC_BIND=close, OMP_PLACES=cores), best of two, on three dedicated AWS instances. Speedup is each machine's own 1→N ratio on the parallel binary; the stock single-core wall and the Group-A serial gain are listed per machine so the two layers can be composed.

• AMD EPYC — Zen 5 "Turin" — 32 cores · SMT off · stock 191.9 s · serial gain 1.67× · AWS c8a.8xlarge
• Intel Xeon — "Granite Rapids" — 16 cores · SMT on · stock 272.2 s · serial gain 2.47× · AWS c8i.8xlarge
• AMD EPYC — Zen 2 "Rome" — 16 cores · SMT on · stock 350.6 s · serial gain 2.16× · AWS c5a.8xlarge

Values are speedup vs each machine's own 1 worker. ◁ marks each CPU's peak (beyond it, hyperthreads only contend).

AMD Turin (32 real cores, SMT off) is the scaling champion — still climbing at 24 threads (5.33×) and fastest in absolute wall time (a one-year run in 26.9 s). With the 1.67× serial gain on top, that is ~7.1× end-to-end versus the stock single-core engine.

Intel Granite Rapids has the fastest single core, so it wins the Group-A serial round (2.47×) — but it peaks at exactly 16 workers (its physical-core count, 4.01×) and REGRESSES at 20/24 as threads land on SMT siblings. The Rome box hits the same wall at its 16 cores (3.81×).

Speedup is bounded by physical cores and memory bandwidth, not thread count: this is a memory-bound model, so 32 real cores beat 16-cores-plus-SMT. Pushing past the physical-core count never helps and usually hurts.

The two layers are orthogonal and multiply: the serial gain is largest on the most memory-starved cores (2.47× on Intel), while the thread scaling is largest where there are the most real cores (5.33× on the 32-core AMD). The interactive calculator below composes them for any CPU/core/mode choice.

Estimate your speedup

Pick a CPU, a core count, and a correctness mode to estimate how much faster one engine run (e.g. one calibration evaluation) becomes versus the stock single-core engine. Every curve is directly measured on a dedicated AWS instance (one model per box, threads pinned, best of two runs) with the 57,998-HRU benchmark and daily channel_sd output. Speedup vs stock = stock single-thread wall / final-engine N-thread wall; serial and parallel gains compose multiplicatively.

Total = serial × parallel = stock single-core wall (192 s) ÷ final-engine 8-thread wall (36.2 s), both at the same daily channel_sd output scope. All three curves are directly measured on this CPU.

Total speedup (serial × parallel) versus the stock single-core engine, AWS c8a — AMD EPYC "Turin", 32 cores, full wavefront mode. Dashed line = ideal linear.

Where the time goes: an Intel VTune profile

An Intel VTune hotspots profile of the final parallel engine at 8 threads (180-day window) shows exactly why the ceiling is where it is — and that high CPU usage is not the same as useful work.

Peace River HUC-8 · final engine · 8 threads · 180-day window · Intel VTune hotspots

Figure 2. Intel VTune CPU-time breakdown of the final engine at 8 threads (Peace River, 180-day window): of ~1,267 CPU-seconds, 73% is effective work, 25% is spin (idle threads waiting at barriers), 2% is overhead. The effective work is dominated by serial channel routing (sd_channel_control3), ~7.7× the parallel land phase.

Of the ~1,267 CPU-seconds the eight threads burned, only ~73% was real work — ~25% was threads spinning at barriers, idle, waiting for the serial main stem to finish. That spin is why a CPU monitor shows 80–90% utilization even though the per-thread payoff is bounded: busy cores are not the same as productive cores.

And the real work is lopsided: channel routing (sd_channel_control3, 698 s) costs about 7.7× the HRU land phase (hru_control, 90 s) — and it is the part that runs serially down the main stem. The land phase we parallelized is only ~12% of the compute, so even perfect HRU scaling can't move the total much. The next real lever is the channel network itself (overlapping its levels, or speeding the per-channel water-quality and sediment math), not more threads.

Acting on this, source-line profiling pinned the single hottest line in the entire engine: the stream-temperature routine (ch_temp) was re-zeroing the ENTIRE landscape-unit array — every routing unit in the watershed — once per channel, every day, even though each channel only reads its own few units. Resetting only the units a channel actually uses (an output-identical fix) barely changes the serial time but markedly raised the multi-thread ceiling, because the cut was on the serial main stem and shrank the Amdahl serial fraction. The methodology working: profile, remove overhead, re-measure.

Re-profiling the fixed build confirms the shift: ch_temp has dropped out of the hotspots entirely (its hottest line went from ~731 s to ~0.8 s), and no channel routine is dominant anymore. What now dominates is not computation but synchronization — with the heavy serial channel work removed, the worker threads finish the wide land-phase level quickly and then sit idle at the per-level barriers while the single-width main stem is walked. That barrier spin (~25% of CPU time at 8 threads) is the new frontier — the next gains come from smarter scheduling (task-dataflow over the routing network), not from shaving more arithmetic.

Correctness: byte-identical at one thread; model-equivalent above, with one localized race

How we check: thread-count invariance — exact byte-level diffs of the daily output files between independent runs at different thread counts, and 1-thread vs the stock serial engine. 'Bit-identical' means every value matched to the last digit; a record is counted as differing if it is off by more than 1e-6, a deliberately strict tripwire, not a measure of magnitude.

• At 1 thread the parallel engine is byte-identical to the unmodified stock serial engine — verified bit-for-bit on a one-year run — so an exact reference run is always available. The wavefront activates only at 2+ threads, so single-thread keeps the original execution order exactly.
• At higher thread counts the hydrology and dissolved constituents are model-equivalent: flow, sediment, nitrate (NO₃) and soluble-phosphorus basin loads agree to <2% (mostly <0.5%) — the residual is floating-point reassociation from the order in which upstream hydrographs are summed, not a change in the modeled physics.
• The HRU land phase is bit-identical across thread counts — water balance (hru_wb), nutrient balance (hru_nb) and sediment losses (hru_ls) all match run-to-run.
• A non-ipo build is bit-identical to the serial source; the ~0.008% residual under the production -ipo build is compiler reassociation, not a difference in physics.

The open item: a localized particulate race, and a byte-identical mode that sidesteps it

With in-stream water quality active over a multi-year run, thread-count invariance exposes a 15–20% divergence in exactly three PARTICULATE / sediment-associated species — organic-N, ammonia and sediment-P. Run-to-run variation confirms it is a genuine data race, not reassociation. The dissolved species (NO₃, soluble-P) and the entire HRU land phase are spared, which localizes the race to the channel in-stream constituent-transport path.

The root cause is the compiler, not the algorithm. ThreadSanitizer traces it not to a missing threadprivate — every per-object buffer is privatized — but to ifx allocating the temporaries of the 152-byte hyd_output derived type (operator results, function returns, argument copy-ins) in shared static module storage rather than on the per-thread stack, so concurrent channels clobber the same static temporary. It is invariant to every compiler control we tried, and a field-by-field source rewrite is byte-identical at one thread but does not remove it.

Because the race lives entirely in the routing phase, the engine offers a byte-identical mode: run the land phase in parallel and the routing phase serially (SWATPLUS_ROUTING_SERIAL=1). That is race-free and bit-reproducible — including in-stream water quality — to the floating-point-reassociation floor. The cost is the routing phase's parallelism: byte-identical mode scales ~2.5× at 8 threads versus the full wavefront's ~4× (about a 37% forfeit at 8 threads). So practitioners pick per run: byte-identical for in-stream water quality, the full wavefront (bit-exact for flow/sediment/dissolved species) otherwise.

The races we found and fixed

Thread-count-invariance testing surfaced several shared 'current-object' scratch variables that raced under the parallel wave. Each was given a per-thread copy:

A footnote: a pre-existing upstream bug, surfaced not introduced

Our bounds-checked debug build kept crashing in the stream-temperature routine. The cause turned out to be upstream, not ours: a channel whose first inflow is another channel reads array element hin_d(0), but the array is allocated from index 1. The production binary, built without bounds checking, silently reads the adjacent memory and continues — which is exactly why the water-temperature column is non-deterministic even in the stock engine.

Because stream temperature is output-only in this configuration (it feeds nothing in the water/sediment/nutrient balance), we left the behavior unchanged to match production, and flagged the bounds read for a separate upstream fix.

What remains

• Fix the particulate-constituent race at its compiler root — move the hyd_output derived-type temporaries off shared static storage (a source-level or compiler change) so the full wavefront is byte-identical for in-stream water quality too, removing the ~37% byte-identical-mode forfeit.
• Cut the barrier spin — now the dominant cost. With the channel-compute bottleneck removed, the wave's per-level barriers leave threads idle on the narrow deep levels. Schedule the routing network as a task-dataflow so a downstream channel starts the moment its own upstreams finish.
• Mop up the smaller remaining compute: per-step string name-matching (replace obtyp comparisons with integer codes) and threadprivate-access overhead.
• Confirm over a multi-decade run that the dissolved/flow reassociation residual stays bounded rather than slowly accumulating.

• The SWAT+ engine can be parallelized on a single shared-memory node without changing the science — and at 1 thread it stays byte-identical to the stock production engine.
• The campaign yields two orthogonal products: a machine-portable serial gain (1.67× on c8a, up to 2.47× on memory-starved cores) that accelerates every user on a single core with no behavior change, and thread parallelism that scales a single model to 5.33× at 24 cores — together ~7.1× end-to-end versus the stock single-core engine.
• Scaling is bounded by physical cores and memory bandwidth, not thread count; the VTune profile shows serial channel routing costs ~7.7× the parallel land phase, so the main-stem routing — not more threads — is the remaining lever.
• Correctness is validated by thread-count invariance: flow, sediment, nitrate and soluble-P are model-equivalent (<2%, mostly <0.5%); a localized, compiler-induced race in three particulate species (organic-N, ammonia, sediment-P) is handled by a byte-identical routing-serial mode, so in-stream water-quality runs stay bit-reproducible at a measured ~37% scaling cost.