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