!> @file spatial_discretization_2d.f90
!> @brief Unsplit method-of-lines spatial residual for 2D Euler.
!! R(Q) = -( (Fx_{i+1/2,j}-Fx_{i-1/2,j})/dx + (Gy_{i,j+1/2}-Gy_{i,j-1/2})/dy ).
!! Per-line component-wise reconstruction (reused via state % reconstruct) + 4-component
!! face flux (euler_riemann_2d). Characteristic projection is a later refinement.
!!
!! Bounds note: state%ub has declared bounds (neq2d, 1-h:nx+h, 1-h:ny+h).  Stencil
!! gathers are inlined (not delegated to assumed-shape subroutine dummies) so that
!! Fortran component-access retains those explicit lower bounds and halo indices like
!! -2 map to the correct ghost cells.
module spatial_discretization_2d
  use precision, only: wp
  use solver_state_2d, only: solver_state_2d_t, neq2d
  use boundary_2d, only: apply_halos_2d, wall_edges_2d
  use euler_riemann_2d, only: face_flux_2d, face_flux_2d_normal, wall_flux_2d_normal
  use euler_physics_2d, only: split_flux_2d, split_contravariant_2d
  use option_registry, only: method_fdm
  use mpi_runtime, only: parallel_fatal
  implicit none
  private
  public :: compute_resid_2d

contains

  !> Allocate (or re-allocate on a dimension change) the reconstruction scratch
  !! arrays shared by the uniform and curvilinear residual paths. Consistent
  !! stat= handling: the two paths previously diverged (one omitted stat=).
  subroutine ensure_resid_scratch_2d(state, nx, ny, sw)
    type(solver_state_2d_t), intent(inout) :: state
    integer, intent(in) :: nx, ny, sw
    integer :: astat
    if (allocated(state % rs_stencil)) then
      if (size(state % rs_stencil, 2) /= sw) then
        deallocate (state % rs_stencil, stat=astat)
        if (astat /= 0) error stop 'spatial_discretization_2d: rs_stencil deallocation failed'
      end if
    end if
    if (.not. allocated(state % rs_stencil)) then
      allocate (state % rs_stencil(neq2d, sw), stat=astat)
      if (astat /= 0) error stop 'spatial_discretization_2d: rs_stencil allocation failed'
    end if
    if (allocated(state % rs_fx)) then
      if (size(state % rs_fx, 2) /= nx + 1) then
        deallocate (state % rs_fx, stat=astat)
        if (astat /= 0) error stop 'spatial_discretization_2d: rs_fx deallocation failed'
      end if
    end if
    if (.not. allocated(state % rs_fx)) then
      allocate (state % rs_fx(neq2d, nx + 1), stat=astat)
      if (astat /= 0) error stop 'spatial_discretization_2d: rs_fx allocation failed'
    end if
    if (allocated(state % rs_gy)) then
      if (size(state % rs_gy, 2) /= ny + 1) then
        deallocate (state % rs_gy, stat=astat)
        if (astat /= 0) error stop 'spatial_discretization_2d: rs_gy deallocation failed'
      end if
    end if
    if (.not. allocated(state % rs_gy)) then
      allocate (state % rs_gy(neq2d, ny + 1), stat=astat)
      if (astat /= 0) error stop 'spatial_discretization_2d: rs_gy allocation failed'
    end if
  end subroutine ensure_resid_scratch_2d

  !> Zhang-Shu positivity limiter for a 2D Euler face state (neq2d=4).
  !! Blends the reconstructed face state toward the cell average until
  !! density rho >= eps and pressure p >= eps.  A no-op in smooth regions.
  pure subroutine limit_pos_2d(qf, qa, gam)
    real(wp), intent(inout) :: qf(neq2d)   !< reconstructed face state
    real(wp), intent(in) :: qa(neq2d)   !< cell-average conserved state
    real(wp), intent(in) :: gam
    real(wp), parameter :: eps = 1.0e-13_wp
    real(wp) :: rho_f, rho_a, p_f, p_a, ke_f, ke_a, theta, denom

    rho_f = qf(1); rho_a = qa(1)

    ! Step 1: enforce rho >= eps
    if (rho_f < eps) then
      ! Guard denominator: when rho_a == rho_f the blend is degenerate; θ → 1.
      denom = rho_a - rho_f
      ! if/else, not merge — merge evaluates both args, so the division would
      ! still trap under -ffpe-trap when denom == 0.
      if (abs(denom) > tiny(1.0_wp)) then
        theta = (rho_a - eps) / denom
      else
        theta = 1.0_wp
      end if
      theta = min(1.0_wp, max(0.0_wp, theta))
      qf = qa + theta * (qf - qa)
    end if

    ! Step 2: enforce p >= eps  (p = (E - 0.5*(m_x^2+m_y^2)/rho)*(gam-1))
    ! Floor rho_f to eps before dividing (covers degenerate rho_a==rho_f<eps).
    rho_f = max(qf(1), eps)
    qf(1) = rho_f
    ke_f = 0.5_wp * (qf(2)**2 + qf(3)**2) / rho_f
    p_f = (qf(4) - ke_f) * (gam - 1.0_wp)
    if (p_f < eps) then
      ke_a = 0.5_wp * (qa(2)**2 + qa(3)**2) / rho_a
      p_a = (qa(4) - ke_a) * (gam - 1.0_wp)
      denom = p_a - p_f
      if (abs(denom) > tiny(1.0_wp)) then
        theta = (p_a - eps) / denom
      else
        theta = 1.0_wp
      end if
      theta = min(1.0_wp, max(0.0_wp, theta))
      qf = qa + theta * (qf - qa)
    end if
  end subroutine limit_pos_2d

  subroutine compute_resid_2d(state)
    type(solver_state_2d_t), intent(inout) :: state
    integer :: i, j, k, nx, ny, sw, so, col, row
    real(wp) :: qL(neq2d), qR(neq2d)
    logical :: limit

    if (.not. associated(state % reconstruct)) &
      error stop 'compute_resid_2d: reconstruction scheme not initialised (call init_recon_scheme_2d first)'

    nx = state % nx_local; ny = state % ny_local
    ! ---- method-aware dispatch ----
    ! FDM is uniform-Cartesian only. FDM+FVS routes to its own operator
    ! (Phase 6, stubbed below); FDM+FDS reuses the shared uniform body that
    ! follows -- the nodal FDS operator IS the FVM uniform operator. FVM
    ! (method never == fdm) keeps the existing uniform/curvilinear split.
    if (trim(state % blocks(1) % method) == method_fdm) then
      if (.not. state % mesh % uniform) then
        ! Curvilinear FDM (nodal transformation metrics).  FDS routes to the
        ! contravariant nodal residual; FVS routes to the contravariant
        ! flux-vector-splitting nodal residual.
        if (state % use_fds) then
          call compute_resid_fdm_curv_2d(state)
        else
          call compute_resid_fdm_curv_fvs_2d(state)
        end if
        return
      end if
      if (.not. state % use_fds) then
        call compute_resid_fdm_fvs_2d(state)   ! Phase 6
        return
      end if
      ! fdm + fds (uniform): fall through to the shared uniform body below.
    else if (.not. state % mesh % uniform) then
      call compute_resid_2d_curvilinear(state)
      return
    end if
    sw = state % stencil_width; so = state % stencil_start_offset
    limit = state % cfg % use_positivity_limiter
    call ensure_resid_scratch_2d(state, nx, ny, sw)

    associate (fx => state % rs_fx, gy => state % rs_gy, stencil => state % rs_stencil)
      call apply_halos_2d(state)
      state % resid(:, 1:nx, 1:ny) = 0.0_wp

      ! ---- x sweep: faces i = 1..nx+1 (face i is between cell i-1 and cell i) ----
      ! Left-biased stencil for qL at face i:  cells i+so, i+so+1, ..., i+so+sw-1
      !   (so < 0 for standard schemes, e.g. weno5: so=-3 => i-3..i+1)
      ! Right-biased stencil for qR at face i: cells i-so-1, i-so-2, ..., i-so-sw
      !   (mirror, e.g. weno5: so=-3 => i+2, i+1, ..., i-2)
      ! Inlined to preserve state%ub's explicit bounds (1-h:nx_local+h).
      do j = 1, ny
        do i = 1, nx + 1
          ! Left-biased stencil: stencil(:,k) = ub(:, i+so + (k-1), j)
          do k = 1, sw
            col = i + so + (k - 1)
            stencil(:, k) = state % ub(:, col, j)
          end do
          call state % reconstruct(stencil, qL)

          ! Right-biased stencil: stencil(:,k) = ub(:, i-so-1 - (k-1), j)
          do k = 1, sw
            col = i - so - 1 - (k - 1)
            stencil(:, k) = state % ub(:, col, j)
          end do
          call state % reconstruct(stencil, qR)

          ! Zhang-Shu positivity: blend toward the left/right cell averages.
          ! Guards removed: at np>1 face 1 and face nx_local+1 are interior
          ! global faces; halos fill ub(:,0,j) and ub(:,nx_local+1,j) so both
          ! cell-average references are always valid.
          if (limit) call limit_pos_2d(qL, state % ub(:, i - 1, j), state % cfg % gam)
          if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam)

          ! Non-positive state guard (mirrors 1D FDS path and 2D FVS path).
          ! Check both qL and qR for non-positive density or pressure before
          ! calling face_flux_2d, so a reconstruction overshoot aborts cleanly
          ! instead of producing a nonphysical flux silently.
          block
            real(wp) :: rho_l, p_l, rho_r, p_r
            ! Check density BEFORE computing pressure — the pressure formula
            ! divides by rho, so a non-positive rho must abort first (else the
            ! division itself would trap under -ffpe-trap).
            rho_l = qL(1)
            rho_r = qR(1)
            if (rho_l <= 0.0_wp .or. rho_r <= 0.0_wp) &
              call parallel_fatal('compute_resid_2d: non-positive density in FDS ' &
                                  //'reconstructed face (x-sweep); enable use_positivity_limiter')
            p_l = (qL(4) - 0.5_wp * (qL(2)**2 + qL(3)**2) / rho_l) * (state % cfg % gam - 1.0_wp)
            p_r = (qR(4) - 0.5_wp * (qR(2)**2 + qR(3)**2) / rho_r) * (state % cfg % gam - 1.0_wp)
            if (p_l <= 0.0_wp .or. p_r <= 0.0_wp) &
              call parallel_fatal('compute_resid_2d: non-positive pressure in FDS ' &
                                  //'reconstructed face (x-sweep); enable use_positivity_limiter')
          end block

          call face_flux_2d(qL, qR, 1, state % cfg % flux_scheme, state % cfg % gam, fx(:, i))
        end do
        do i = 1, nx
          state % resid(:, i, j) = state % resid(:, i, j) &
              & - (fx(:, i + 1) - fx(:, i)) / state % dx
        end do
      end do

      ! ---- y sweep: faces j = 1..ny+1 (face j is between cell j-1 and cell j) ----
      do i = 1, nx
        do j = 1, ny + 1
          ! Left-biased stencil along y
          do k = 1, sw
            row = j + so + (k - 1)
            stencil(:, k) = state % ub(:, i, row)
          end do
          call state % reconstruct(stencil, qL)

          ! Right-biased stencil along y
          do k = 1, sw
            row = j - so - 1 - (k - 1)
            stencil(:, k) = state % ub(:, i, row)
          end do
          call state % reconstruct(stencil, qR)

          ! Zhang-Shu positivity: blend toward the bottom/top cell averages.
          ! Guards removed: at np>1 face 1 and face ny_local+1 are interior
          ! global faces; halos fill ub(:,i,0) and ub(:,i,ny_local+1) so both
          ! cell-average references are always valid.
          if (limit) call limit_pos_2d(qL, state % ub(:, i, j - 1), state % cfg % gam)
          if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam)

          ! Non-positive state guard (mirrors 1D FDS path and 2D FVS path).
          block
            real(wp) :: rho_l, p_l, rho_r, p_r
            ! Check density BEFORE computing pressure (see x-sweep note).
            rho_l = qL(1)
            rho_r = qR(1)
            if (rho_l <= 0.0_wp .or. rho_r <= 0.0_wp) &
              call parallel_fatal('compute_resid_2d: non-positive density in FDS ' &
                                  //'reconstructed face (y-sweep); enable use_positivity_limiter')
            p_l = (qL(4) - 0.5_wp * (qL(2)**2 + qL(3)**2) / rho_l) * (state % cfg % gam - 1.0_wp)
            p_r = (qR(4) - 0.5_wp * (qR(2)**2 + qR(3)**2) / rho_r) * (state % cfg % gam - 1.0_wp)
            if (p_l <= 0.0_wp .or. p_r <= 0.0_wp) &
              call parallel_fatal('compute_resid_2d: non-positive pressure in FDS ' &
                                  //'reconstructed face (y-sweep); enable use_positivity_limiter')
          end block

          call face_flux_2d(qL, qR, 2, state % cfg % flux_scheme, state % cfg % gam, gy(:, j))
        end do
        do j = 1, ny
          state % resid(:, i, j) = state % resid(:, i, j) &
              & - (gy(:, j + 1) - gy(:, j)) / state % dy
        end do
      end do
    end associate
  end subroutine compute_resid_2d

  !> Curvilinear residual: R = -(1/V)[ (Fhat_{i+1}-Fhat_i) + (Ghat_{j+1}-Ghat_j) ]
  !! with contravariant face fluxes Fhat = |S_xi| * f_n(qL,qR; n_xi),
  !! Ghat = |S_eta| * f_n(qL,qR; n_eta). Reconstruction is identical to the
  !! uniform path (uniform computational stencil); geometry enters via s_xi/s_eta/vol.
  subroutine compute_resid_2d_curvilinear(state)
    type(solver_state_2d_t), intent(inout) :: state
    integer :: i, j, k, nx, ny, sw, so, col, row
    real(wp) :: qL(neq2d), qR(neq2d), nrm(2), smag, fn(neq2d)
    logical :: limit, w_xlo, w_xhi, w_ylo, w_yhi

    nx = state % nx_local; ny = state % ny_local
    sw = state % stencil_width; so = state % stencil_start_offset
    limit = state % cfg % use_positivity_limiter
    call ensure_resid_scratch_2d(state, nx, ny, sw)

    associate (fx => state % rs_fx, gy => state % rs_gy, stencil => state % rs_stencil, &
               mesh => state % mesh)
      call apply_halos_2d(state)
      ! Physical reflecting-wall edges get a weak pressure-only wall flux at the
      ! boundary face (zero mass/energy flux, momentum = p*n) instead of a Riemann
      ! solve on the reflected ghost; see euler_riemann_2d % wall_flux_2d_normal.
      call wall_edges_2d(state, w_xlo, w_xhi, w_ylo, w_yhi)
      state % resid(:, 1:nx, 1:ny) = 0.0_wp

      ! ---- xi sweep: faces i = 1..nx+1 (ξ-face at node column i) ----
      do j = 1, ny
        do i = 1, nx + 1
          do k = 1, sw
            col = i + so + (k - 1)
            stencil(:, k) = state % ub(:, col, j)
          end do
          call state % reconstruct(stencil, qL)
          do k = 1, sw
            col = i - so - 1 - (k - 1)
            stencil(:, k) = state % ub(:, col, j)
          end do
          call state % reconstruct(stencil, qR)
          if (limit) call limit_pos_2d(qL, state % ub(:, i - 1, j), state % cfg % gam)
          if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam)
          smag = sqrt(mesh % s_xi(1, i, j)**2 + mesh % s_xi(2, i, j)**2)
          ! A degenerate face metric is rank-local data; error stop here would
          ! strand the peers in the next collective (halo exchange / par_sum) —
          ! tear down collectively like the rest of this module (audit 2026-07-06 N4).
          if (smag <= 0.0_wp) call parallel_fatal('compute_resid_2d_curvilinear: degenerate xi-face')
          nrm = mesh % s_xi(:, i, j) / smag
          if (w_xlo .and. i == 1) then
            call wall_flux_2d_normal(qR, nrm, state % cfg % gam, fn)        ! interior side = qR
          else if (w_xhi .and. i == nx + 1) then
            call wall_flux_2d_normal(qL, nrm, state % cfg % gam, fn)        ! interior side = qL
          else
            call face_flux_2d_normal(qL, qR, nrm, state % cfg % flux_scheme, state % cfg % gam, fn)
          end if
          fx(:, i) = smag * fn
        end do
        do i = 1, nx
          state % resid(:, i, j) = state % resid(:, i, j) &
              & - (fx(:, i + 1) - fx(:, i)) / mesh % vol(i, j)
        end do
      end do

      ! ---- eta sweep: faces j = 1..ny+1 (η-face at node row j) ----
      do i = 1, nx
        do j = 1, ny + 1
          do k = 1, sw
            row = j + so + (k - 1)
            stencil(:, k) = state % ub(:, i, row)
          end do
          call state % reconstruct(stencil, qL)
          do k = 1, sw
            row = j - so - 1 - (k - 1)
            stencil(:, k) = state % ub(:, i, row)
          end do
          call state % reconstruct(stencil, qR)
          if (limit) call limit_pos_2d(qL, state % ub(:, i, j - 1), state % cfg % gam)
          if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam)
          smag = sqrt(mesh % s_eta(1, i, j)**2 + mesh % s_eta(2, i, j)**2)
          ! A degenerate face metric is rank-local data; error stop here would
          ! strand the peers in the next collective (halo exchange / par_sum) —
          ! tear down collectively like the rest of this module (audit 2026-07-06 N4).
          if (smag <= 0.0_wp) call parallel_fatal('compute_resid_2d_curvilinear: degenerate eta-face')
          nrm = mesh % s_eta(:, i, j) / smag
          if (w_ylo .and. j == 1) then
            call wall_flux_2d_normal(qR, nrm, state % cfg % gam, fn)        ! interior side = qR
          else if (w_yhi .and. j == ny + 1) then
            call wall_flux_2d_normal(qL, nrm, state % cfg % gam, fn)        ! interior side = qL
          else
            call face_flux_2d_normal(qL, qR, nrm, state % cfg % flux_scheme, state % cfg % gam, fn)
          end if
          gy(:, j) = smag * fn
        end do
        do j = 1, ny
          state % resid(:, i, j) = state % resid(:, i, j) &
              & - (gy(:, j + 1) - gy(:, j)) / mesh % vol(i, j)
        end do
      end do
    end associate
  end subroutine compute_resid_2d_curvilinear

  !> Curvilinear-FDM (FDS) residual on NODAL transformation metrics:
  !! R = -(1/J)[ (Fhat_{i+1/2}-Fhat_{i-1/2}) + (Ghat_{j+1/2}-Ghat_{j-1/2}) ]
  !! with contravariant face fluxes Fhat = |S_xi| * f_n(qL,qR; n_xi),
  !! Ghat = |S_eta| * f_n(qL,qR; n_eta) and S_xi=mesh%sx_xi, S_eta=mesh%sx_eta,
  !! J=mesh%jac.  Structurally identical to compute_resid_2d_curvilinear (the FVM
  !! curvilinear residual); only the metric SOURCE differs — FVM face-area vectors
  !! (s_xi/s_eta/vol) become the nodal transformation metrics here.  Reconstruction
  !! and the contravariant flux machinery are reused verbatim.  Freestream/GCL
  !! preservation (the Phase 1 gate) relies on the consistent metric evaluation in
  !! compute_metrics_fdm_2d: a constant flow telescopes the face fluxes to ~0.
  !!
  !! BCs: apply_halos_2d routes the FDM method to the nodal kernel apply_bcs_fdm_2d
  !! (duplicate-endpoint-aware periodic wrap, node-mirror reflecting).  Reflecting
  !! walls reuse the curvilinear weak pressure-wall flux with the nodal metric
  !! normal (Phase 4 refines the curvilinear nodal wall); periodic BCs (the
  !! Phase 1 gate) need no wall handling.
  subroutine compute_resid_fdm_curv_2d(state)
    type(solver_state_2d_t), intent(inout) :: state
    integer :: i, j, k, nx, ny, sw, so, col, row
    real(wp) :: qL(neq2d), qR(neq2d), nrm(2), smag, fn(neq2d)
    logical :: limit, w_xlo, w_xhi, w_ylo, w_yhi

    nx = state % nx_local; ny = state % ny_local
    sw = state % stencil_width; so = state % stencil_start_offset
    limit = state % cfg % use_positivity_limiter
    call ensure_resid_scratch_2d(state, nx, ny, sw)

    associate (fx => state % rs_fx, gy => state % rs_gy, stencil => state % rs_stencil, &
               mesh => state % mesh)
      call apply_halos_2d(state)
      call wall_edges_2d(state, w_xlo, w_xhi, w_ylo, w_yhi)
      state % resid(:, 1:nx, 1:ny) = 0.0_wp

      ! ---- xi sweep: faces i = 1..nx+1 (xi-face between node i-1 and node i) ----
      do j = 1, ny
        do i = 1, nx + 1
          do k = 1, sw
            col = i + so + (k - 1)
            stencil(:, k) = state % ub(:, col, j)
          end do
          call state % reconstruct(stencil, qL)
          do k = 1, sw
            col = i - so - 1 - (k - 1)
            stencil(:, k) = state % ub(:, col, j)
          end do
          call state % reconstruct(stencil, qR)
          if (limit) call limit_pos_2d(qL, state % ub(:, i - 1, j), state % cfg % gam)
          if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam)
          smag = sqrt(mesh % sx_xi(1, i, j)**2 + mesh % sx_xi(2, i, j)**2)
          ! A degenerate face metric is rank-local data; error stop here would
          ! strand the peers in the next collective (halo exchange / par_sum) —
          ! tear down collectively like the rest of this module (audit 2026-07-06 N4).
          if (smag <= 0.0_wp) call parallel_fatal('compute_resid_fdm_curv_2d: degenerate xi-face')
          nrm = mesh % sx_xi(:, i, j) / smag
          if (w_xlo .and. i == 1) then
            call wall_flux_2d_normal(qR, nrm, state % cfg % gam, fn)        ! interior side = qR
          else if (w_xhi .and. i == nx + 1) then
            call wall_flux_2d_normal(qL, nrm, state % cfg % gam, fn)        ! interior side = qL
          else
            call face_flux_2d_normal(qL, qR, nrm, state % cfg % flux_scheme, state % cfg % gam, fn)
          end if
          fx(:, i) = smag * fn
        end do
        do i = 1, nx
          state % resid(:, i, j) = state % resid(:, i, j) &
              & - (fx(:, i + 1) - fx(:, i)) / mesh % jac(i, j)
        end do
      end do

      ! ---- eta sweep: faces j = 1..ny+1 (eta-face between node j-1 and node j) ----
      do i = 1, nx
        do j = 1, ny + 1
          do k = 1, sw
            row = j + so + (k - 1)
            stencil(:, k) = state % ub(:, i, row)
          end do
          call state % reconstruct(stencil, qL)
          do k = 1, sw
            row = j - so - 1 - (k - 1)
            stencil(:, k) = state % ub(:, i, row)
          end do
          call state % reconstruct(stencil, qR)
          if (limit) call limit_pos_2d(qL, state % ub(:, i, j - 1), state % cfg % gam)
          if (limit) call limit_pos_2d(qR, state % ub(:, i, j), state % cfg % gam)
          smag = sqrt(mesh % sx_eta(1, i, j)**2 + mesh % sx_eta(2, i, j)**2)
          ! A degenerate face metric is rank-local data; error stop here would
          ! strand the peers in the next collective (halo exchange / par_sum) —
          ! tear down collectively like the rest of this module (audit 2026-07-06 N4).
          if (smag <= 0.0_wp) call parallel_fatal('compute_resid_fdm_curv_2d: degenerate eta-face')
          nrm = mesh % sx_eta(:, i, j) / smag
          if (w_ylo .and. j == 1) then
            call wall_flux_2d_normal(qR, nrm, state % cfg % gam, fn)        ! interior side = qR
          else if (w_yhi .and. j == ny + 1) then
            call wall_flux_2d_normal(qL, nrm, state % cfg % gam, fn)        ! interior side = qL
          else
            call face_flux_2d_normal(qL, qR, nrm, state % cfg % flux_scheme, state % cfg % gam, fn)
          end if
          gy(:, j) = smag * fn
        end do
        do j = 1, ny
          state % resid(:, i, j) = state % resid(:, i, j) &
              & - (gy(:, j + 1) - gy(:, j)) / mesh % jac(i, j)
        end do
      end do
    end associate
  end subroutine compute_resid_fdm_curv_2d

  !> Curvilinear-FDM FLUX-VECTOR-SPLITTING residual on NODAL transformation
  !! metrics (the FVS counterpart of compute_resid_fdm_curv_2d):
  !!   R = -(1/J)[ (Fhat_{i+1/2}-Fhat_{i-1/2}) + (Ghat_{j+1/2}-Ghat_{j-1/2}) ]
  !! with the contravariant split-flux reconstruction
  !!   Fhat_{i+1/2} = recon_L[ Fhat^+ ] + recon_R[ Fhat^- ],
  !! mirroring the uniform FVS residual compute_resid_fdm_fvs_2d but with
  !! CONTRAVARIANT split fluxes and the 1/J divisor.
  !!
  !! FREESTREAM (the headline gate): the metric used to SPLIT each stencil node is
  !! FROZEN at the FACE metric (mesh%sx_xi(:,i,j) for face i / mesh%sx_eta(:,i,j)
  !! for face j) — the same GCL-telescoping face metric the FDS branch uses — NOT
  !! a per-node metric.  Then for a constant free-stream the split is identical at
  !! every stencil node, the (constant-reproducing) reconstruction returns it
  !! unchanged, and the face flux collapses to Fhat_face = S_face . (F,G).  The
  !! conservative difference therefore telescopes through the discrete GCL identity
  !! (Delta_xi S_xi + Delta_eta S_eta = 0) to machine zero, exactly as the FDS
  !! branch does.  A naive per-node-metric split would instead leak the
  !! reconstruction's truncation error into the free-stream residual.
  !!
  !! Cost note: because the metric is per-face, the split is recomputed for each
  !! stencil node at each face (2*sw split_contravariant_2d calls per face) rather
  !! than precomputed once per node as in the uniform FVS path.  This keeps the
  !! formulation freestream-consistent and high-order; the LF split (no LAPACK)
  !! used by the accuracy/positivity gates keeps the cost modest.
  !!
  !! BCs: the full nodal BC set is applied via apply_halos_2d's nodal kernel
  !! (apply_bcs_fdm_2d): periodic, outflow, supersonic_inlet, characteristic
  !! farfield, and the curvilinear slip wall (reflection about the boundary-node
  !! metric normal).  This is the production FVS residual for the curvilinear
  !! showcases (cylinder/airfoil/blunt body) and the bump/ramp validation.  Unlike
  !! the FDS path it applies no weak wall-flux at the boundary face (the weak
  !! pressure wall is FDS-only by design); the nodal ghost reflection is its wall
  !! treatment.  No positivity limiter is applied (FVS reconstructs fluxes, not
  !! states), mirroring compute_resid_fdm_fvs_2d.
  subroutine compute_resid_fdm_curv_fvs_2d(state)
    type(solver_state_2d_t), intent(inout) :: state
    integer :: i, j, k, nx, ny, sw, so, col, row
    real(wp) :: fpL(neq2d), fmR(neq2d), gam, sxi(2), seta(2), fp_k(neq2d), fm_k(neq2d)

    nx = state % nx_local; ny = state % ny_local
    sw = state % stencil_width; so = state % stencil_start_offset
    gam = state % cfg % gam
    call ensure_resid_scratch_2d(state, nx, ny, sw)

    associate (fx => state % rs_fx, gy => state % rs_gy, stencil => state % rs_stencil, &
               ub => state % ub, mesh => state % mesh)
      call apply_halos_2d(state)
      state % resid(:, 1:nx, 1:ny) = 0.0_wp

      ! ---- xi sweep: faces i = 1..nx+1, face metric S_xi = mesh%sx_xi(:,i,j) ----
      do j = 1, ny
        do i = 1, nx + 1
          sxi = mesh % sx_xi(:, i, j)
          ! Left-biased stencil on Fhat^+ (split each stencil node with the FACE metric).
          do k = 1, sw
            col = i + so + (k - 1)
            call fvs_vacuum_guard_2d(ub(:, col, j), gam, 'fdm-curv-fvs xi-sweep')
            call split_contravariant_2d(ub(:, col, j), sxi, state % cfg % flux_scheme, gam, fp_k, fm_k)
            stencil(:, k) = fp_k
          end do
          call state % reconstruct(stencil, fpL)
          ! Right-biased (mirror) stencil on Fhat^-.
          do k = 1, sw
            col = i - so - 1 - (k - 1)
            call fvs_vacuum_guard_2d(ub(:, col, j), gam, 'fdm-curv-fvs xi-sweep')
            call split_contravariant_2d(ub(:, col, j), sxi, state % cfg % flux_scheme, gam, fp_k, fm_k)
            stencil(:, k) = fm_k
          end do
          call state % reconstruct(stencil, fmR)
          fx(:, i) = fpL + fmR
        end do
        do i = 1, nx
          state % resid(:, i, j) = state % resid(:, i, j) &
              & - (fx(:, i + 1) - fx(:, i)) / mesh % jac(i, j)
        end do
      end do

      ! ---- eta sweep: faces j = 1..ny+1, face metric S_eta = mesh%sx_eta(:,i,j) ----
      do i = 1, nx
        do j = 1, ny + 1
          seta = mesh % sx_eta(:, i, j)
          do k = 1, sw
            row = j + so + (k - 1)
            call fvs_vacuum_guard_2d(ub(:, i, row), gam, 'fdm-curv-fvs eta-sweep')
            call split_contravariant_2d(ub(:, i, row), seta, state % cfg % flux_scheme, gam, fp_k, fm_k)
            stencil(:, k) = fp_k
          end do
          call state % reconstruct(stencil, fpL)
          do k = 1, sw
            row = j - so - 1 - (k - 1)
            call fvs_vacuum_guard_2d(ub(:, i, row), gam, 'fdm-curv-fvs eta-sweep')
            call split_contravariant_2d(ub(:, i, row), seta, state % cfg % flux_scheme, gam, fp_k, fm_k)
            stencil(:, k) = fm_k
          end do
          call state % reconstruct(stencil, fmR)
          gy(:, j) = fpL + fmR
        end do
        do j = 1, ny
          state % resid(:, i, j) = state % resid(:, i, j) &
              & - (gy(:, j + 1) - gy(:, j)) / mesh % jac(i, j)
        end do
      end do
    end associate
  end subroutine compute_resid_fdm_curv_fvs_2d

  !> Allocate (or re-allocate on a dimension change) the FVS split-flux scratch
  !! fp_all/fm_all, halo-padded like ub (1-h:nx+h, 1-h:ny+h).  Allocated lazily
  !! from compute_resid_fdm_fvs_2d so non-FVS paths never pay the memory cost.
  subroutine ensure_fvs_scratch_2d(state, nx, ny, h)
    type(solver_state_2d_t), intent(inout) :: state
    integer, intent(in) :: nx, ny, h
    integer :: astat
    logical :: need

    need = .false.
    if (allocated(state % fp_all)) then
      if (size(state % fp_all, 2) /= nx + 2 * h .or. &
          size(state % fp_all, 3) /= ny + 2 * h) then
        deallocate (state % fp_all, state % fm_all, stat=astat)
        if (astat /= 0) error stop 'spatial_discretization_2d: fvs scratch deallocation failed'
        need = .true.
      end if
    else
      need = .true.
    end if
    if (need) then
      allocate (state % fp_all(neq2d, 1 - h:nx + h, 1 - h:ny + h), stat=astat)
      if (astat /= 0) error stop 'spatial_discretization_2d: fp_all allocation failed'
      allocate (state % fm_all(neq2d, 1 - h:nx + h, 1 - h:ny + h), stat=astat)
      if (astat /= 0) error stop 'spatial_discretization_2d: fm_all allocation failed'
    end if
  end subroutine ensure_fvs_scratch_2d

  !> Vacuum/positivity guard for the FVS split-flux precompute.  The FVS split
  !! routines evaluate sqrt(gam*p/rho); a non-positive density or pressure would
  !! produce NaN that then propagates silently through the reconstructed face
  !! fluxes and the residual.  Detect it here and fail clearly instead.
  !! Mirrors the 1D FVS guard in spatial_discretization.f90 (block-guarded
  !! check + parallel_fatal).  Operands are guarded with nested ifs because
  !! Fortran .and. does not short-circuit (rho is a divisor in p).
  subroutine fvs_vacuum_guard_2d(q, gam, ctx)
    real(wp), intent(in) :: q(neq2d), gam
    character(len=*), intent(in) :: ctx
    real(wp) :: rho_i, p_i
    rho_i = q(1)
    if (rho_i <= 0.0_wp) then
      call parallel_fatal('compute_resid_fdm_fvs_2d: non-positive density in FVS precompute ('//ctx//')')
    else
      p_i = (q(4) - 0.5_wp * (q(2)**2 + q(3)**2) / rho_i) * (gam - 1.0_wp)
      if (p_i <= 0.0_wp) &
        call parallel_fatal('compute_resid_fdm_fvs_2d: non-positive pressure in FVS precompute ('//ctx//')')
    end if
  end subroutine fvs_vacuum_guard_2d

  !> FDM + flux-vector-splitting residual (Shu-style conservative finite
  !! difference).  Per direction: split the physical flux at every node
  !! (interior + halo) into F^+/F^-, reconstruct F^+ with the LEFT-biased
  !! stencil and F^- with the RIGHT-biased (mirror) stencil, sum the two
  !! reconstructions at each face, then take the conservative difference.
  !!
  !! The stencil offsets are IDENTICAL to the uniform FDS sweeps in
  !! compute_resid_2d -- only the operand differs: FDS reconstructs the
  !! conserved state (left -> qL, right -> qR); FVS reconstructs the split
  !! fluxes (left -> F^+, right -> F^-).  Matching offsets is what preserves
  !! the high-order accuracy (Task 6.3 convergence gate).
  !!
  !! No positivity limiter is applied: FVS reconstructs fluxes, not states, so
  !! the face-state Zhang-Shu guard used by the FDS path does not apply here
  !! (mirrors the 1D FVS path, which also leaves split fluxes unlimited).
  subroutine compute_resid_fdm_fvs_2d(state)
    type(solver_state_2d_t), intent(inout) :: state
    integer :: i, j, k, nx, ny, sw, so, h, col, row
    real(wp) :: fpL(neq2d), fmR(neq2d), gam

    nx = state % nx_local; ny = state % ny_local
    sw = state % stencil_width; so = state % stencil_start_offset
    h = state % halo_width
    gam = state % cfg % gam
    call ensure_resid_scratch_2d(state, nx, ny, sw)
    call ensure_fvs_scratch_2d(state, nx, ny, h)

    associate (fx => state % rs_fx, gy => state % rs_gy, stencil => state % rs_stencil, &
               fp => state % fp_all, fm => state % fm_all, ub => state % ub)
      call apply_halos_2d(state)
      state % resid(:, 1:nx, 1:ny) = 0.0_wp

      ! ---- x sweep: split the x-flux at every node read by the x stencils ----
      ! Precompute over i = 1-h..nx+h (all columns the face stencils reach) and
      ! j = 1..ny (only interior rows are referenced in x).  Corner halos are
      ! never touched, so non-physical corner ghosts cannot poison split_flux_2d.
      do j = 1, ny
        do i = 1 - h, nx + h
          call fvs_vacuum_guard_2d(ub(:, i, j), gam, 'x-sweep')
          call split_flux_2d(ub(:, i, j), 1, state % cfg % flux_scheme, gam, &
              & fp(:, i, j), fm(:, i, j))
        end do
      end do
      do j = 1, ny
        do i = 1, nx + 1
          ! Left-biased stencil on F^+ (same offsets as the FDS qL stencil).
          do k = 1, sw
            col = i + so + (k - 1)
            stencil(:, k) = fp(:, col, j)
          end do
          call state % reconstruct(stencil, fpL)

          ! Right-biased (mirror) stencil on F^- (same offsets as FDS qR).
          do k = 1, sw
            col = i - so - 1 - (k - 1)
            stencil(:, k) = fm(:, col, j)
          end do
          call state % reconstruct(stencil, fmR)

          fx(:, i) = fpL + fmR
        end do
        do i = 1, nx
          state % resid(:, i, j) = state % resid(:, i, j) &
              & - (fx(:, i + 1) - fx(:, i)) / state % dx
        end do
      end do

      ! ---- y sweep: split the y-flux at every node read by the y stencils ----
      ! Precompute over i = 1..nx (interior columns) and j = 1-h..ny+h.
      do j = 1 - h, ny + h
        do i = 1, nx
          call fvs_vacuum_guard_2d(ub(:, i, j), gam, 'y-sweep')
          call split_flux_2d(ub(:, i, j), 2, state % cfg % flux_scheme, gam, &
              & fp(:, i, j), fm(:, i, j))
        end do
      end do
      do i = 1, nx
        do j = 1, ny + 1
          ! Left-biased stencil on F^+ along y.
          do k = 1, sw
            row = j + so + (k - 1)
            stencil(:, k) = fp(:, i, row)
          end do
          call state % reconstruct(stencil, fpL)

          ! Right-biased (mirror) stencil on F^- along y.
          do k = 1, sw
            row = j - so - 1 - (k - 1)
            stencil(:, k) = fm(:, i, row)
          end do
          call state % reconstruct(stencil, fmR)

          gy(:, j) = fpL + fmR
        end do
        do j = 1, ny
          state % resid(:, i, j) = state % resid(:, i, j) &
              & - (gy(:, j + 1) - gy(:, j)) / state % dy
        end do
      end do
    end associate
  end subroutine compute_resid_fdm_fvs_2d

end module spatial_discretization_2d
