!> @file euler_physics_2d.f90
!> @brief 2D Euler physics: flux vectors F (x) and G (y), eigensystem, wave speed.
!!
!! Conserved state Q = (rho, rho*u, rho*v, E).
!!   p = (gam-1)(E - 0.5*rho*(u^2 + v^2))
!!   F = (rho u, rho u^2 + p, rho u v,   u(E+p))
!!   G = (rho v, rho u v,     rho v^2+p, v(E+p))
module euler_physics_2d
  use precision, only: wp
  use solver_state_2d, only: neq2d, solver_state_2d_t
  use parallel_reductions, only: par_max_real
  use mpi_runtime, only: parallel_fatal
  implicit none
  private
  public :: primitives_2d, compute_euler_flux_2d, compute_eigensystem_2d, compute_max_wave_speed_2d, &
            split_flux_2d, split_contravariant_2d

contains

  !> Extract primitives (rho, u, v, p) from a 2D conserved state.
  pure subroutine primitives_2d(q, gam, rho, u, v, p)
    real(wp), intent(in) :: q(:), gam
    real(wp), intent(out) :: rho, u, v, p
    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))
  end subroutine primitives_2d

  !> Compute both directional flux vectors F (x) and G (y) for state Q.
  pure subroutine compute_euler_flux_2d(q, fx, gy, gam)
    real(wp), intent(in) :: q(:)
    real(wp), intent(out) :: fx(:), gy(:)
    real(wp), intent(in) :: gam
    real(wp) :: rho, u, v, p

    call primitives_2d(q, gam, rho, u, v, p)

    fx(1) = q(2)
    fx(2) = rho * u * u + p
    fx(3) = rho * u * v
    fx(4) = u * (q(4) + p)

    gy(1) = q(3)
    gy(2) = rho * u * v
    gy(3) = rho * v * v + p
    gy(4) = v * (q(4) + p)
  end subroutine compute_euler_flux_2d

  !> Right-eigenvector matrix K and its inverse for the flux Jacobian in
  !! direction `dir` (1 = x: eigenvalues u-c,u,u,u+c; 2 = y: v-c,v,v,v+c).
  !! Columns of K, conserved variables Q=(rho,rho u,rho v,E):
  !!   x:  [1, u-c, v, H-u c] [1, u, v, q2] [0, 0, 1, v] [1, u+c, v, H+u c]
  !!   y:  [1, u, v-c, H-v c] [1, u, v, q2] [0, 1, 0, u] [1, u, v+c, H+v c]
  !! where H=(E+p)/rho, q2=0.5*(u^2+v^2). K^{-1} via LAPACK getrf/getri.
  subroutine compute_eigensystem_2d(q, dir, r_mat, r_inv, gam)
    real(wp), intent(in) :: q(neq2d)
    integer, intent(in) :: dir
    real(wp), intent(out) :: r_mat(neq2d, neq2d), r_inv(neq2d, neq2d)
    real(wp), intent(in) :: gam
    real(wp) :: rho, u, v, p, c, h_tot, q2

    if (dir /= 1 .and. dir /= 2) error stop 'euler_physics_2d: dir must be 1 (x) or 2 (y)'

    call primitives_2d(q, gam, rho, u, v, p)
    c = sqrt(gam * p / rho)
    h_tot = (q(4) + p) / rho
    q2 = 0.5_wp * (u * u + v * v)

    if (dir == 1) then
      ! eigenvalue order: u-c, u (entropy), u (shear), u+c
      r_mat(:, 1) = [1.0_wp, u - c, v, h_tot - u * c]
      r_mat(:, 2) = [1.0_wp, u, v, q2]
      r_mat(:, 3) = [0.0_wp, 0.0_wp, 1.0_wp, v]
      r_mat(:, 4) = [1.0_wp, u + c, v, h_tot + u * c]
    else
      ! eigenvalue order: v-c, v (entropy), v (shear), v+c
      r_mat(:, 1) = [1.0_wp, u, v - c, h_tot - v * c]
      r_mat(:, 2) = [1.0_wp, u, v, q2]
      r_mat(:, 3) = [0.0_wp, 1.0_wp, 0.0_wp, u]
      r_mat(:, 4) = [1.0_wp, u, v + c, h_tot + v * c]
    end if

    call invert_4x4(r_mat, r_inv)
  end subroutine compute_eigensystem_2d

  !> Inverse of a 4x4 matrix via LAPACK LU (dgetrf) + inversion (dgetri).
  subroutine invert_4x4(a, ainv)
    real(wp), intent(in) :: a(4, 4)
    real(wp), intent(out) :: ainv(4, 4)
    real(wp) :: work(16)
    integer :: ipiv(4), info
    interface
      subroutine dgetrf(m, n, a_in, lda, ipiv_in, info_out)
        integer, intent(in) :: m, n, lda
        double precision, intent(inout) :: a_in(lda, *)
        integer, intent(out) :: ipiv_in(*)
        integer, intent(out) :: info_out
      end subroutine dgetrf
      subroutine dgetri(n, a_in, lda, ipiv_in, work_in, lwork, info_out)
        integer, intent(in) :: n, lda, lwork
        double precision, intent(inout) :: a_in(lda, *)
        integer, intent(in) :: ipiv_in(*)
        double precision, intent(out) :: work_in(*)
        integer, intent(out) :: info_out
      end subroutine dgetri
    end interface

    ainv = a
    call dgetrf(4, 4, ainv, 4, ipiv, info)
    if (info /= 0) error stop 'euler_physics_2d: dgetrf failed (singular eigenvector matrix)'
    call dgetri(4, ainv, 4, ipiv, work, 16, info)
    if (info /= 0) error stop 'euler_physics_2d: dgetri failed'
  end subroutine invert_4x4

  !> Global maximum signal speed max_{i,j}( max(|u|,|v|) + c ) over interior cells.
  !!
  !! Loops the per-rank LOCAL interior (1:nx_local, 1:ny_local) — `ub` is
  !! allocated with local extent — then reduces with par_max_real so every rank
  !! returns the same global maximum. Mirrors compute_dt_2d. (At np=1
  !! nx_local==nx, so the serial answer is unchanged.)
  function compute_max_wave_speed_2d(state) result(s_max)
    type(solver_state_2d_t), intent(in) :: state
    real(wp) :: s_max
    integer :: i, j
    real(wp) :: rho, u, v, p, c

    s_max = 0.0_wp
    do j = 1, state % ny_local
      do i = 1, state % nx_local
        call primitives_2d(state % ub(:, i, j), state % cfg % gam, rho, u, v, p)
        c = sqrt(state % cfg % gam * p / rho)
        s_max = max(s_max, max(abs(u), abs(v)) + c)
      end do
    end do
    s_max = par_max_real(s_max)
  end function compute_max_wave_speed_2d

  !> Return the directional Euler flux: fx for dir==1, gy for dir==2.
  pure subroutine directional_flux(q, dir, gam, f)
    real(wp), intent(in) :: q(neq2d), gam
    integer, intent(in) :: dir
    real(wp), intent(out) :: f(neq2d)
    real(wp) :: fx(neq2d), gy(neq2d)
    call compute_euler_flux_2d(q, fx, gy, gam)
    if (dir == 1) then
      f = fx
    else
      f = gy
    end if
  end subroutine directional_flux

  !> Lax-Friedrichs directional splitting: F± = ½(F ± α Q), α = |vel_dir| + c.
  !! Guarantees fp + fm == F exactly (same α applied to both halves).
  pure subroutine lax_friedrichs_split_2d(q, dir, gam, fp, fm)
    real(wp), intent(in) :: q(neq2d), gam
    integer, intent(in) :: dir
    real(wp), intent(out) :: fp(neq2d), fm(neq2d)
    real(wp) :: rho, u, v, p, c, alpha, fx(neq2d), gy(neq2d), f(neq2d)
    call primitives_2d(q, gam, rho, u, v, p)
    c = sqrt(gam * p / rho)
    call compute_euler_flux_2d(q, fx, gy, gam)
    if (dir == 1) then
      alpha = abs(u) + c; f = fx
    else
      alpha = abs(v) + c; f = gy
    end if
    fp = 0.5_wp * (f + alpha * q)
    fm = 0.5_wp * (f - alpha * q)
  end subroutine lax_friedrichs_split_2d

  !> Steger-Warming splitting: F± = R·diag(λ±)·R⁻¹·Q.
  !! Exact for the homogeneous Euler flux (F = A·Q); fp+fm == F to machine precision.
  !! Eigenvalue order: x → (u-c, u, u, u+c); y → (v-c, v, v, v+c).
  !! NOT pure: calls compute_eigensystem_2d which uses LAPACK.
  subroutine steger_warming_split_2d(q, dir, gam, fp, fm)
    real(wp), intent(in) :: q(neq2d), gam
    integer, intent(in) :: dir
    real(wp), intent(out) :: fp(neq2d), fm(neq2d)
    real(wp) :: rmat(neq2d, neq2d), rinv(neq2d, neq2d)
    real(wp) :: rho, u, v, p, c, lam(neq2d), lp(neq2d), lm(neq2d), tmp(neq2d)
    call primitives_2d(q, gam, rho, u, v, p)
    c = sqrt(gam * p / rho)
    call compute_eigensystem_2d(q, dir, rmat, rinv, gam)
    if (dir == 1) then
      lam = [u - c, u, u, u + c]
    else
      lam = [v - c, v, v, v + c]
    end if
    lp = 0.5_wp * (lam + abs(lam))
    lm = 0.5_wp * (lam - abs(lam))
    tmp = matmul(rinv, q)          ! characteristic amplitudes: R^{-1} Q
    fp = matmul(rmat, lp * tmp)    ! R diag(λ+) R^{-1} Q
    fm = matmul(rmat, lm * tmp)    ! R diag(λ-) R^{-1} Q
  end subroutine steger_warming_split_2d

  !> van Leer splitting: Mach-split mass flux carrying transverse momentum.
  !! Pressure split: p± = p(M±1)²(2∓M)/4 (subsonic); supersonic branches are
  !! fully upwind (fp==F, fm==0 for M≥1; fp==0, fm==F for M≤-1).
  !! Uses `in` (normal) and `it` (transverse) momentum indices to handle both
  !! directions uniformly.
  pure subroutine van_leer_split_2d(q, dir, gam, fp, fm)
    real(wp), intent(in) :: q(neq2d), gam
    integer, intent(in) :: dir
    real(wp), intent(out) :: fp(neq2d), fm(neq2d)
    real(wp) :: rho, u, v, p, c, h_tot, vn, vt, mach, mdp, mdm, psp, psm
    integer :: in, it   ! normal / transverse momentum component indices
    call primitives_2d(q, gam, rho, u, v, p)
    c = sqrt(gam * p / rho)
    h_tot = (q(4) + p) / rho
    if (dir == 1) then
      vn = u; vt = v; in = 2; it = 3
    else
      vn = v; vt = u; in = 3; it = 2
    end if
    mach = vn / c
    if (mach >= 1.0_wp) then
      call directional_flux(q, dir, gam, fp); fm = 0.0_wp; return
    else if (mach <= -1.0_wp) then
      call directional_flux(q, dir, gam, fm); fp = 0.0_wp; return
    end if
    mdp = rho * c * (mach + 1.0_wp)**2 / 4.0_wp
    psp = p * (mach + 1.0_wp)**2 * (2.0_wp - mach) / 4.0_wp
    fp(1) = mdp; fp(in) = mdp * vn + psp; fp(it) = mdp * vt; fp(4) = mdp * h_tot
    mdm = -rho * c * (mach - 1.0_wp)**2 / 4.0_wp
    psm = p * (mach - 1.0_wp)**2 * (2.0_wp + mach) / 4.0_wp
    fm(1) = mdm; fm(in) = mdm * vn + psm; fm(it) = mdm * vt; fm(4) = mdm * h_tot
  end subroutine van_leer_split_2d

  !> Contravariant flux-vector splitting along a general metric/area vector
  !! `smetric` = S_xi or S_eta (size 2).  Splits the CONTRAVARIANT flux
  !!   Fhat = smetric(1)*F + smetric(2)*G   (F,G physical Euler fluxes)
  !! into fp/fm with eigenvalues U, U, U +- c|S| where the contravariant velocity
  !! is U = smetric(1)*u + smetric(2)*v and |S| = norm2(smetric).
  !!
  !! Realized by ROTATING the state into the metric-normal frame n = S/|S|,
  !! applying the proven axis split_flux_2d (dir=1, eigenvalues u_n -+ c), rotating
  !! the split fluxes back, and scaling by |S|.  Because the rotation is
  !! orthogonal and |S| > 0, this is algebraically IDENTICAL to the metric-normal
  !! Steger-Warming/van-Leer/Lax-Friedrichs split (eigenvalues |S|*(u_n -+ c) =
  !! U -+ c|S|), so it inherits:
  !!   * consistency  fp + fm == Fhat  (rotational invariance of the Euler flux);
  !!   * axis reduction: with smetric=(s,0), s>0, n=(1,0) and the rotation is the
  !!     identity, so the result is exactly s * split_flux_2d(...,dir=1,...).
  !! NOT pure: may reach steger_warming_split_2d (LAPACK) via split_flux_2d.
  subroutine split_contravariant_2d(q, smetric, scheme, gam, fp, fm)
    real(wp), intent(in) :: q(neq2d), smetric(2), gam
    character(len=*), intent(in) :: scheme
    real(wp), intent(out) :: fp(neq2d), fm(neq2d)
    real(wp) :: smag, nx, ny, qrot(neq2d), fprot(neq2d), fmrot(neq2d)

    smag = sqrt(smetric(1) * smetric(1) + smetric(2) * smetric(2))
    ! Rank-local data inside the per-substage residual: tear down
    ! collectively rather than stranding peers in the next halo exchange /
    ! reduction (audit 2026-07-06 N4).
    if (smag <= 0.0_wp) call parallel_fatal('euler_physics_2d: split_contravariant_2d degenerate metric vector')
    nx = smetric(1) / smag
    ny = smetric(2) / smag

    ! Rotate the momentum into the metric-normal frame (rho, E invariant).
    qrot(1) = q(1)
    qrot(2) = nx * q(2) + ny * q(3)
    qrot(3) = -ny * q(2) + nx * q(3)
    qrot(4) = q(4)

    call split_flux_2d(qrot, 1, scheme, gam, fprot, fmrot)

    ! Rotate the split fluxes back and scale by |S| (eigenvalues -> U -+ c|S|).
    fp(1) = smag * fprot(1)
    fp(2) = smag * (nx * fprot(2) - ny * fprot(3))
    fp(3) = smag * (ny * fprot(2) + nx * fprot(3))
    fp(4) = smag * fprot(4)
    fm(1) = smag * fmrot(1)
    fm(2) = smag * (nx * fmrot(2) - ny * fmrot(3))
    fm(3) = smag * (ny * fmrot(2) + nx * fmrot(3))
    fm(4) = smag * fmrot(4)
  end subroutine split_contravariant_2d

  !> Dispatch to one of the three FVS routines by scheme name.
  !! dir=1 → x-splitting (uses F, eigenvalues u±c); dir=2 → y-splitting (uses G, v±c).
  !! scheme: 'lax_friedrichs' | 'steger_warming' | 'van_leer'.
  !! NOT pure: may call steger_warming_split_2d (which calls LAPACK).
  subroutine split_flux_2d(q, dir, scheme, gam, fp, fm)
    real(wp), intent(in) :: q(neq2d), gam
    integer, intent(in) :: dir
    character(len=*), intent(in) :: scheme
    real(wp), intent(out) :: fp(neq2d), fm(neq2d)
    select case (trim(scheme))
    case ('lax_friedrichs'); call lax_friedrichs_split_2d(q, dir, gam, fp, fm)
    case ('steger_warming'); call steger_warming_split_2d(q, dir, gam, fp, fm)
    case ('van_leer'); call van_leer_split_2d(q, dir, gam, fp, fm)
    case default; error stop 'split_flux_2d: unknown FVS scheme "'//trim(scheme)//'"'
    end select
  end subroutine split_flux_2d

end module euler_physics_2d