!> @file stepper_kernels.f90 !> @brief Dimension-blind explicit Runge-Kutta stepper kernels (R2, spec §5). !! !! Each kernel advances one time step over a FLAT rank-1 view of the whole !! (halo-padded) solution array, calling back into the dimension's residual !! assembly through `rhs`. Stage arithmetic is copied VERBATIM from the !! historical 1D steppers — the groupings are load-bearing for the project's !! bit-for-bit refactor invariant; do not re-associate them. !! !! Halo-zero invariant (required of every caller): the whole-array stage !! updates span halo cells, so they advance halos as ub_halo + c*dt*resid_halo. !! This is correct only because `resid` halo cells are zero — zeroed once at !! allocation and never written by the residual assembly (which fills interior !! cells only). Halo `u` values left stale by the updates are refreshed by !! halo exchange + BC application inside the next `rhs` call. !! !! Aliasing: `rhs(ctx)` writes the same memory `resid` points at (and reads !! the memory `u` points at). The array dummies are POINTERS so this external !! modification is conforming; do not change them to assumed-shape dummies. module stepper_kernels use precision, only: wp implicit none private public :: rhs_fn public :: kernel_euler, kernel_ssprk22, kernel_tvd_rk3, kernel_rk4, kernel_ssprk54 abstract interface !> Recompute the spatial residual for the current solution state. !! Implementations refresh halos/BCs and fill resid interior cells. !! Implementations must not modify the context's dt mid-step — kernels !! freeze dt at entry. subroutine rhs_fn(ctx) class(*), intent(inout), target :: ctx end subroutine rhs_fn end interface contains !> Explicit Euler: Q^{n+1} = Q^n + dt R(Q^n). subroutine kernel_euler(u, resid, dt, rhs, ctx) real(wp), pointer, intent(in) :: u(:), resid(:) real(wp), intent(in) :: dt procedure(rhs_fn) :: rhs class(*), intent(inout), target :: ctx call rhs(ctx) u = u + dt * resid end subroutine kernel_euler !> SSPRK(2,2) (Shu & Osher 1988). 1D-verbatim stage expressions. subroutine kernel_ssprk22(u, resid, s1, dt, rhs, ctx) real(wp), pointer, intent(in) :: u(:), resid(:), s1(:) real(wp), intent(in) :: dt procedure(rhs_fn) :: rhs class(*), intent(inout), target :: ctx s1 = u call rhs(ctx) u = s1 + dt * resid call rhs(ctx) u = 0.5_wp * s1 + 0.5_wp * u + 0.5_wp * dt * resid end subroutine kernel_ssprk22 !> TVD-RK3 (Shu & Osher 1988). 1D-verbatim stage expressions. subroutine kernel_tvd_rk3(u, resid, s1, dt, rhs, ctx) real(wp), pointer, intent(in) :: u(:), resid(:), s1(:) real(wp), intent(in) :: dt procedure(rhs_fn) :: rhs class(*), intent(inout), target :: ctx s1 = u call rhs(ctx) u = s1 + dt * resid call rhs(ctx) u = 0.75_wp * s1 + 0.25_wp * u + 0.25_wp * dt * resid call rhs(ctx) u = 1.0_wp / 3.0_wp * s1 + 2.0_wp / 3.0_wp * u + 2.0_wp / 3.0_wp * dt * resid end subroutine kernel_tvd_rk3 !> Classic RK4 (not SSP). s1 holds Q^n; s2 accumulates the weighted stages. subroutine kernel_rk4(u, resid, s1, s2, dt, rhs, ctx) real(wp), pointer, intent(in) :: u(:), resid(:), s1(:), s2(:) real(wp), intent(in) :: dt procedure(rhs_fn) :: rhs class(*), intent(inout), target :: ctx s1 = u s2 = 0.0_wp call rhs(ctx) s2 = s2 + (1.0_wp / 6.0_wp) * dt * resid u = s1 + 0.5_wp * dt * resid call rhs(ctx) s2 = s2 + (2.0_wp / 6.0_wp) * dt * resid u = s1 + 0.5_wp * dt * resid call rhs(ctx) s2 = s2 + (2.0_wp / 6.0_wp) * dt * resid u = s1 + dt * resid call rhs(ctx) s2 = s2 + (1.0_wp / 6.0_wp) * dt * resid u = s1 + s2 end subroutine kernel_rk4 !> SSPRK(5,4) — Spiteri & Ruuth (2002), SIAM J. Numer. Anal. 40(2), 469-491, !! Shu-Osher form. Coefficients independently cross-checked against NodePy !! (ketch/nodepy, runge_kutta_method.py, RK['SSP54']/'SSPRK54', fetched !! 2026-07-08): converting this routine's Shu-Osher coefficients to Butcher !! (A, b) form reproduces every one of NodePy's 30-digit A/b entries to !! 8-12 significant digits, consistent with NodePy's own documented note !! that its table is a higher-precision (30-digit) re-solve of the same !! Ruuth-Spiteri method whose originally published coefficients (as coded !! below, ~15 digits) are only accurate to about 10 digits. Five stages, !! order 4, SSP coefficient c ~= 1.508. !! Register use: s1 = u0 (stages 2-4), s2 = u2 then the final stage's !! partial sum, s3 = u3; the final stage's L(u3) piece is accumulated into !! s2 before stage 4 overwrites resid (avoids a fourth scratch array). !! HISTORY: before 2026-07 this routine carried an unpublished 3rd-order !! tableau (initial-commit transcription defect; see the 2026-07-07 !! retrospective). The observed-order tripwire test guards against any !! regression. subroutine kernel_ssprk54(u, resid, s1, s2, s3, dt, rhs, ctx) real(wp), pointer, intent(in) :: u(:), resid(:), s1(:), s2(:), s3(:) real(wp), intent(in) :: dt procedure(rhs_fn) :: rhs class(*), intent(inout), target :: ctx real(wp), parameter :: b10 = 0.391752226571890_wp real(wp), parameter :: a20 = 0.444370493651235_wp real(wp), parameter :: a21 = 0.555629506348765_wp real(wp), parameter :: b21 = 0.368410593050371_wp real(wp), parameter :: a30 = 0.620101851488403_wp real(wp), parameter :: a32 = 0.379898148511597_wp real(wp), parameter :: b32 = 0.251891774271694_wp real(wp), parameter :: a40 = 0.178079954393132_wp real(wp), parameter :: a43 = 0.821920045606868_wp real(wp), parameter :: b43 = 0.544974750228521_wp real(wp), parameter :: a52 = 0.517231671970585_wp real(wp), parameter :: a53 = 0.096059710526147_wp real(wp), parameter :: b53 = 0.063692468666290_wp real(wp), parameter :: a54 = 0.386708617503269_wp real(wp), parameter :: b54 = 0.226007483236906_wp ! Register schedule (three scratch arrays, five rhs calls): ! s1 = u0 (consumed by stages 2-4, never needed after) ! s2 = u2, then reused as the final-stage partial sum ! s3 = u3 ! L(u3) is consumed TWICE (stage 4 and the final stage); the final ! stage's u2/u3/L(u3) piece is therefore accumulated into s2 BEFORE ! stage 4 overwrites resid — this avoids a fourth scratch array. ! The full schedule was verified in 40-digit arithmetic to reproduce the ! plain published recurrence exactly (pin: 1.8096749104456358498...). s1 = u ! u0 call rhs(ctx) u = s1 + b10 * dt * resid ! u1 call rhs(ctx) u = a20 * s1 + a21 * u + b21 * dt * resid ! u2 s2 = u ! save u2 call rhs(ctx) u = a30 * s1 + a32 * u + b32 * dt * resid ! u3 s3 = u ! save u3 call rhs(ctx) ! resid = L(u3) ! accumulate the final stage's u3-residual piece NOW, before stage 4 ! overwrites resid: stash a52*u2 + a53*u3 + b53*dt*L(u3) into s2. s2 = a52 * s2 + a53 * s3 + b53 * dt * resid u = a40 * s1 + a43 * s3 + b43 * dt * resid ! u4 call rhs(ctx) ! resid = L(u4) u = s2 + a54 * u + b54 * dt * resid ! u5 end subroutine kernel_ssprk54 end module stepper_kernels