!> @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
