!> @file euler_riemann_2d.f90
!> @brief 4-component 2D Euler face fluxes (Rusanov, HLLC) with face-normal rotation.
!!
!! Q = (rho, rho*u, rho*v, E). For a face with normal in direction dir (1=x, 2=y),
!! normal-momentum index mn and tangential mt:
!!   dir=1 -> mn=2 (rho*u), mt=3 (rho*v); dir=2 -> mn=3 (rho*v), mt=2 (rho*u).
!! Pressure uses BOTH velocity components (full KE) -- this is why the 3-component
!! 1D solvers cannot be reused directly.
module euler_riemann_2d
  use precision, only: wp
  use solver_state_2d, only: neq2d
  use option_registry, only: flux_rusanov
  implicit none
  private
  public :: face_flux_2d, face_flux_2d_normal, wall_flux_2d_normal

contains

  !> Inviscid solid-wall (no-penetration) face flux through unit normal nrm.
  !! At a solid wall v.n = 0, so the Euler flux collapses to the pressure term:
  !!   F = [ 0, p*nx, p*ny, 0 ],
  !! with p the interior-side pressure and (nx,ny) the boundary-face normal.
  !! Mass and energy flux are EXACTLY zero by construction (a "weak"/pressure
  !! wall), so no spurious mass crosses the wall and no spurious entropy is
  !! generated by an acoustic over-pressure -- unlike a ghost-reflection wall,
  !! whose Riemann star pressure carries an O(rho*c*|v_n|) error wherever the
  !! reconstructed interior normal velocity is non-zero (always, on a curved
  !! wall). `q` is the interior-side (reconstructed) conserved state; `gam` the
  !! ratio of specific heats. nrm need not be normalised in sign: the residual's
  !! +xi/+eta divergence orientation makes the resulting wall pressure force act
  !! into the domain at both min and max edges.
  pure subroutine wall_flux_2d_normal(q, nrm, gam, flux)
    real(wp), intent(in) :: q(neq2d), nrm(2), gam
    real(wp), intent(out) :: flux(neq2d)
    real(wp) :: p
    p = (gam - 1.0_wp) * (q(4) - 0.5_wp * (q(2)**2 + q(3)**2) / q(1))
    p = max(p, 0.0_wp)   ! defense-in-depth: a non-physical face must not pull inward
    flux(1) = 0.0_wp
    flux(2) = p * nrm(1)
    flux(3) = p * nrm(2)
    flux(4) = 0.0_wp
  end subroutine wall_flux_2d_normal

  pure subroutine prim_split(q, dir, gam, rho, un, ut, p, c, mn, mt)
    real(wp), intent(in) :: q(neq2d), gam
    integer, intent(in) :: dir
    real(wp), intent(out) :: rho, un, ut, p, c
    integer, intent(out) :: mn, mt
    real(wp) :: u, v
    rho = q(1)
    u = q(2) / rho
    v = q(3) / rho
    p = (gam - 1.0_wp) * (q(4) - 0.5_wp * rho * (u * u + v * v))
    ! Clamp the sound-speed radicand: a non-physical (negative-pressure)
    ! reconstructed face would otherwise give c = NaN that poisons the whole
    ! HLLC flux (sL/sR/star region). Mirrors the c_roe clamp in hllc_flux_2d;
    ! production paths already abort upstream on non-positive p (FVS/FDS guards),
    ! so this is defense-in-depth that keeps the flux finite for any input.
    c = sqrt(max(gam * p / rho, 0.0_wp))
    if (dir == 1) then
      un = u; ut = v; mn = 2; mt = 3
    else
      un = v; ut = u; mn = 3; mt = 2
    end if
  end subroutine prim_split

  pure subroutine normal_phys_flux(q, dir, gam, f)
    real(wp), intent(in) :: q(neq2d), gam
    integer, intent(in) :: dir
    real(wp), intent(out) :: f(neq2d)
    real(wp) :: rho, un, ut, p, c
    integer :: mn, mt
    call prim_split(q, dir, gam, rho, un, ut, p, c, mn, mt)
    f(1) = rho * un
    f(mn) = rho * un * un + p
    f(mt) = rho * un * ut
    f(4) = un * (q(4) + p)
  end subroutine normal_phys_flux

  subroutine face_flux_2d(qL, qR, dir, scheme, gam, flux)
    real(wp), intent(in) :: qL(neq2d), qR(neq2d), gam
    integer, intent(in) :: dir
    character(len=*), intent(in) :: scheme
    real(wp), intent(out) :: flux(neq2d)
    if (dir /= 1 .and. dir /= 2) error stop 'euler_riemann_2d: dir must be 1 or 2'
    select case (trim(scheme))
    case (flux_rusanov)
      call rusanov_flux_2d(qL, qR, dir, gam, flux)
    case default
      call hllc_flux_2d(qL, qR, dir, gam, flux)
    end select
  end subroutine face_flux_2d

  !> Face flux through an arbitrary UNIT normal nrm=(nx,ny). Rotates the
  !! conserved momentum into the face frame, reuses the axis solver (dir=1),
  !! then rotates the momentum flux back to (x,y). nrm must be a unit vector.
  subroutine face_flux_2d_normal(qL, qR, nrm, scheme, gam, flux)
    real(wp), intent(in) :: qL(neq2d), qR(neq2d), nrm(2), gam
    character(len=*), intent(in) :: scheme
    real(wp), intent(out) :: flux(neq2d)
    real(wp) :: qLr(neq2d), qRr(neq2d), fr(neq2d), nx, ny
    nx = nrm(1); ny = nrm(2)
    ! Rotate conserved momentum into (normal, tangential).
    qLr(1) = qL(1)
    qLr(2) = qL(2) * nx + qL(3) * ny
    qLr(3) = -qL(2) * ny + qL(3) * nx
    qLr(4) = qL(4)
    qRr(1) = qR(1)
    qRr(2) = qR(2) * nx + qR(3) * ny
    qRr(3) = -qR(2) * ny + qR(3) * nx
    qRr(4) = qR(4)
    ! Proven axis solver treats index 2 as normal-momentum, 3 as tangential.
    call face_flux_2d(qLr, qRr, 1, scheme, gam, fr)
    ! Rotate the momentum flux back to (x, y).
    flux(1) = fr(1)
    flux(2) = fr(2) * nx - fr(3) * ny
    flux(3) = fr(2) * ny + fr(3) * nx
    flux(4) = fr(4)
  end subroutine face_flux_2d_normal

  pure subroutine rusanov_flux_2d(qL, qR, dir, gam, flux)
    real(wp), intent(in) :: qL(neq2d), qR(neq2d), gam
    integer, intent(in) :: dir
    real(wp), intent(out) :: flux(neq2d)
    real(wp) :: fL(neq2d), fR(neq2d)
    real(wp) :: rhoL, unL, utL, pL, cL, rhoR, unR, utR, pR, cR, smax
    integer :: mn, mt
    call prim_split(qL, dir, gam, rhoL, unL, utL, pL, cL, mn, mt)
    call prim_split(qR, dir, gam, rhoR, unR, utR, pR, cR, mn, mt)
    call normal_phys_flux(qL, dir, gam, fL)
    call normal_phys_flux(qR, dir, gam, fR)
    smax = max(abs(unL) + cL, abs(unR) + cR)
    flux = 0.5_wp * (fL + fR) - 0.5_wp * smax * (qR - qL)
  end subroutine rusanov_flux_2d

  pure subroutine hllc_flux_2d(qL, qR, dir, gam, flux)
    real(wp), intent(in) :: qL(neq2d), qR(neq2d), gam
    integer, intent(in) :: dir
    real(wp), intent(out) :: flux(neq2d)
    real(wp) :: rhoL, unL, utL, pL, cL, rhoR, unR, utR, pR, cR
    real(wp) :: fL(neq2d), fR(neq2d)
    real(wp) :: sL, sR, sStar, factor, uStar(neq2d)
    real(wp) :: rt, un_roe, ut_roe, c_roe, hL, hR, h_roe
    integer :: mn, mt
    call prim_split(qL, dir, gam, rhoL, unL, utL, pL, cL, mn, mt)
    call prim_split(qR, dir, gam, rhoR, unR, utR, pR, cR, mn, mt)
    call normal_phys_flux(qL, dir, gam, fL)
    call normal_phys_flux(qR, dir, gam, fR)

    hL = (qL(4) + pL) / rhoL
    hR = (qR(4) + pR) / rhoR
    rt = sqrt(rhoR / rhoL)
    un_roe = (unL + rt * unR) / (1.0_wp + rt)
    ut_roe = (utL + rt * utR) / (1.0_wp + rt)
    h_roe = (hL + rt * hR) / (1.0_wp + rt)
    ! Clamp the Roe radicand at 0: for two physical (p>0) states Roe's identity
    ! keeps it >= 0, but a non-physical reconstructed face (negative pressure)
    ! can drive it negative, giving c_roe = NaN that poisons sL/sR and the whole
    ! flux. max(...,0) is a no-op for valid states and lets a bad face yield a
    ! finite flux (the upstream FDS positivity guard then aborts cleanly).
    c_roe = sqrt(max((gam - 1.0_wp) * (h_roe - 0.5_wp * (un_roe * un_roe + ut_roe * ut_roe)), 0.0_wp))
    sL = min(unL - cL, un_roe - c_roe)
    sR = max(unR + cR, un_roe + c_roe)

    if (sL >= 0.0_wp) then
      flux = fL; return
    end if
    if (sR <= 0.0_wp) then
      flux = fR; return
    end if
    sStar = (pR - pL + rhoL * unL * (sL - unL) - rhoR * unR * (sR - unR)) / &
            (rhoL * (sL - unL) - rhoR * (sR - unR))
    if (sStar >= 0.0_wp) then
      factor = rhoL * (sL - unL) / (sL - sStar)
      uStar(1) = factor
      uStar(mn) = factor * sStar
      uStar(mt) = factor * utL
      uStar(4) = factor * (qL(4) / rhoL + (sStar - unL) * (sStar + pL / (rhoL * (sL - unL))))
      flux = fL + sL * (uStar - qL)
    else
      factor = rhoR * (sR - unR) / (sR - sStar)
      uStar(1) = factor
      uStar(mn) = factor * sStar
      uStar(mt) = factor * utR
      uStar(4) = factor * (qR(4) / rhoR + (sStar - unR) * (sStar + pR / (rhoR * (sR - unR))))
      flux = fR + sR * (uStar - qR)
    end if
  end subroutine hllc_flux_2d

end module euler_riemann_2d
