!> @file euler_physics.f90 !> @brief Euler equation physics: flux function, flux splitting, and eigensystem. module euler_physics use precision, only: wp use option_registry, only: flux_lax_friedrichs, flux_steger_warming, flux_van_leer, & flux_ausm_plus, flux_hll, flux_hllc, flux_roe, & flux_scheme_names, join_token_list implicit none private public :: compute_euler_flux, compute_eigensystem, compute_max_wave_speed public :: init_flux_scheme contains ! --------------------------------------------------------------------------- !> Bind the flux dispatch pointers in @p state to the requested scheme. !! !! Valid scheme names: !! 'lax_friedrichs' — Local Lax-Friedrichs splitting [default] !! 'steger_warming' — Steger-Warming eigenvalue splitting !! 'van_leer' — van Leer smooth polynomial splitting !! 'ausm_plus' — AUSM+ FDS scheme, AUSM family (Liou 1996) !! 'hll' — HLL approximate Riemann solver (Harten, Lax & van Leer 1983) !! 'hllc' — HLLC solver with contact restoration (Toro et al. 1994) !! 'roe' — Roe approximate Riemann solver (Roe 1981) !! !! @param[inout] state Solver instance; use_fds, flux_split and !! fds_solver are written. !! @param[in] scheme Name of the flux scheme. ! --------------------------------------------------------------------------- subroutine init_flux_scheme(state, scheme) use solver_state, only: solver_state_t type(solver_state_t), intent(inout) :: state character(len=*), intent(in) :: scheme select case (trim(scheme)) case (flux_lax_friedrichs) state % use_fds = .false. state % flux_split => lax_friedrichs_split state % split_both => lax_friedrichs_split_both case (flux_steger_warming) state % use_fds = .false. state % flux_split => steger_warming_split state % split_both => steger_warming_split_both case (flux_van_leer) state % use_fds = .false. state % flux_split => van_leer_split state % split_both => van_leer_split_both case (flux_ausm_plus) state % use_fds = .true. state % fds_solver => ausm_plus case (flux_hll) state % use_fds = .true. state % fds_solver => hll_flux case (flux_hllc) state % use_fds = .true. state % fds_solver => hllc_flux case (flux_roe) state % use_fds = .true. state % fds_solver => roe_flux case default error stop 'euler_physics: unknown flux scheme "'//trim(scheme)// & '"; valid: '//trim(join_token_list(flux_scheme_names)) end select end subroutine init_flux_scheme ! --------------------------------------------------------------------------- !> Extract primitive variables (rho, u, p) from a conserved state vector. !! !! Inverse of the conserved-to-primitive map for the 1-D Euler equations: !! rho = q(1), u = q(2)/rho, p = (q(3) - ½ρu²)(γ-1). !! !! @param[in] q Conserved state (ρ, ρu, E) !! @param[in] gam Ratio of specific heats γ !! @param[out] rho Density !! @param[out] u Velocity !! @param[out] p Pressure ! --------------------------------------------------------------------------- pure subroutine extract_primitives(q, gam, rho, u, p) real(wp), intent(in) :: q(:) real(wp), intent(in) :: gam real(wp), intent(out) :: rho, u, p rho = q(1) u = q(2) / rho p = (q(3) - 0.5_wp * rho * u**2) * (gam - 1.0_wp) end subroutine extract_primitives ! --------------------------------------------------------------------------- !> Compute Roe-averaged velocity, enthalpy, and sound speed. !! !! Uses √ρ-weighted averages (Roe 1981): !! ũ = (√ρ_L u_L + √ρ_R u_R)/(√ρ_L + √ρ_R), H̃ similarly. !! c̃ = √((γ-1)(H̃ - ũ²/2)) !! !! @param[in] rho_L, u_L, h_L Left-state density, velocity, total specific enthalpy !! @param[in] rho_R, u_R, h_R Right-state density, velocity, total specific enthalpy !! @param[in] gam Ratio of specific heats γ !! @param[out] u_roe Roe-averaged velocity !! @param[out] h_roe Roe-averaged specific enthalpy !! @param[out] c_roe Roe-averaged sound speed ! --------------------------------------------------------------------------- pure subroutine compute_roe_averages(rho_L, u_L, h_L, rho_R, u_R, h_R, gam, & u_roe, h_roe, c_roe) real(wp), intent(in) :: rho_L, u_L, h_L, rho_R, u_R, h_R, gam real(wp), intent(out) :: u_roe, h_roe, c_roe real(wp) :: sqrt_rho_L, sqrt_rho_R, denom sqrt_rho_L = sqrt(rho_L) sqrt_rho_R = sqrt(rho_R) denom = sqrt_rho_L + sqrt_rho_R u_roe = (sqrt_rho_L * u_L + sqrt_rho_R * u_R) / denom h_roe = (sqrt_rho_L * h_L + sqrt_rho_R * h_R) / denom c_roe = sqrt((gam - 1.0_wp) * (h_roe - 0.5_wp * u_roe**2)) end subroutine compute_roe_averages ! --------------------------------------------------------------------------- !> Compute the Euler flux vector F(Q) for a given conserved state. !! !! F = (rho*u, rho*u^2 + p, u*(E + p)) !! See course notes, Eq. (2). !! !! @param[in] q Conserved state vector (rho, rho*u, E) !! @param[out] flux Euler flux vector !! @param[in] gam Ratio of specific heats ! --------------------------------------------------------------------------- pure subroutine compute_euler_flux(q, flux, gam) real(wp), intent(in) :: q(:) real(wp), intent(out) :: flux(:) real(wp), intent(in) :: gam real(wp) :: rho, u, p rho = q(1) u = q(2) / rho ! velocity p = (q(3) - 0.5_wp * rho * u**2) * (gam - 1.0_wp) ! pressure, Eq. (3)-(4) flux(1) = q(2) ! rho * u flux(2) = rho * u**2 + p ! rho * u^2 + p flux(3) = u * (q(3) + p) ! u * (E + p) end subroutine compute_euler_flux ! --------------------------------------------------------------------------- !> Lax-Friedrichs flux splitting. !! !! F^{+} = 0.5 * (F + alpha * Q), F^{-} = 0.5 * (F - alpha * Q) !! where alpha = max|lambda| = |u| + c (local wave speed estimate). !! See course notes, Sec. "Flux Splitting". !! !! @param[in] q Conserved state vector !! @param[out] flux Split flux (positive or negative part) !! @param[in] sign +1 for F^+, -1 for F^- !! @param[in] gam Ratio of specific heats ! --------------------------------------------------------------------------- pure subroutine lax_friedrichs_split(q, flux, sign, gam) real(wp), intent(in) :: q(:) real(wp), intent(out) :: flux(:) integer, intent(in) :: sign real(wp), intent(in) :: gam real(wp) :: rho, u, p, a, alpha call extract_primitives(q, gam, rho, u, p) a = sqrt(gam * p / rho) ! speed of sound alpha = abs(u) + a ! max eigenvalue magnitude ! Full Euler flux flux(1) = q(2) flux(2) = rho * u**2 + p flux(3) = u * (q(3) + p) ! Lax-Friedrichs splitting: F^{+-} = 0.5 * (F +- alpha * Q) flux = 0.5_wp * (flux + real(sign, wp) * alpha * q) end subroutine lax_friedrichs_split ! --------------------------------------------------------------------------- !> Fused Lax-Friedrichs flux splitting: returns F^+ and F^- in one call. !! !! Equivalent to calling lax_friedrichs_split twice (sign=+1 and sign=-1), !! but extracts primitives and computes the Euler flux only once per cell. !! !! @param[in] q Conserved state vector. !! @param[out] fp Positive split flux F^+. !! @param[out] fm Negative split flux F^-. !! @param[in] gam Ratio of specific heats. ! --------------------------------------------------------------------------- pure subroutine lax_friedrichs_split_both(q, fp, fm, gam) real(wp), intent(in) :: q(:), gam real(wp), intent(out) :: fp(:), fm(:) real(wp) :: rho, u, p, a, alpha, f_euler(3) call extract_primitives(q, gam, rho, u, p) a = sqrt(gam * p / rho) alpha = abs(u) + a f_euler(1) = q(2) f_euler(2) = rho * u**2 + p f_euler(3) = u * (q(3) + p) fp = 0.5_wp * (f_euler + alpha * q) fm = 0.5_wp * (f_euler - alpha * q) end subroutine lax_friedrichs_split_both ! --------------------------------------------------------------------------- !> Fused Steger-Warming flux splitting: returns F^+ and F^- in one call. !! !! Extracts primitives and eigenvalues once; computes both positive and !! negative split eigenvalues λᵢ± in a single pass. !! !! @param[in] q Conserved state vector. !! @param[out] fp Positive split flux F^+. !! @param[out] fm Negative split flux F^-. !! @param[in] gam Ratio of specific heats. ! --------------------------------------------------------------------------- pure subroutine steger_warming_split_both(q, fp, fm, gam) real(wp), intent(in) :: q(:), gam real(wp), intent(out) :: fp(:), fm(:) real(wp) :: rho, u, p, c, h_total real(wp) :: lam1, lam2, lam3 real(wp) :: lp1, lp2, lp3 ! λᵢ+ = (λᵢ + |λᵢ|) / 2 real(wp) :: lm1, lm2, lm3 ! λᵢ- = (λᵢ - |λᵢ|) / 2 real(wp) :: coeff ! rho / (2 * gam) call extract_primitives(q, gam, rho, u, p) c = sqrt(gam * p / rho) h_total = (q(3) + p) / rho lam1 = u - c; lam2 = u; lam3 = u + c lp1 = 0.5_wp * (lam1 + abs(lam1)) lp2 = 0.5_wp * (lam2 + abs(lam2)) lp3 = 0.5_wp * (lam3 + abs(lam3)) lm1 = lam1 - lp1 ! = (lam1 - |lam1|) / 2 lm2 = lam2 - lp2 lm3 = lam3 - lp3 coeff = rho / (2.0_wp * gam) fp(1) = coeff * (2.0_wp * (gam - 1.0_wp) * lp2 + lp1 + lp3) fp(2) = coeff * (2.0_wp * (gam - 1.0_wp) * u * lp2 + (u - c) * lp1 + (u + c) * lp3) fp(3) = coeff * ((gam - 1.0_wp) * u**2 * lp2 + (h_total - u * c) * lp1 & + (h_total + u * c) * lp3) fm(1) = coeff * (2.0_wp * (gam - 1.0_wp) * lm2 + lm1 + lm3) fm(2) = coeff * (2.0_wp * (gam - 1.0_wp) * u * lm2 + (u - c) * lm1 + (u + c) * lm3) fm(3) = coeff * ((gam - 1.0_wp) * u**2 * lm2 + (h_total - u * c) * lm1 & + (h_total + u * c) * lm3) end subroutine steger_warming_split_both ! --------------------------------------------------------------------------- !> Fused van Leer flux splitting: returns F^+ and F^- in one call. !! !! Extracts primitives once and evaluates both split fluxes, handling all !! Mach-number regimes (supersonic positive, subsonic, supersonic negative). !! !! @param[in] q Conserved state vector. !! @param[out] fp Positive split flux F^+. !! @param[out] fm Negative split flux F^-. !! @param[in] gam Ratio of specific heats. ! --------------------------------------------------------------------------- pure subroutine van_leer_split_both(q, fp, fm, gam) real(wp), intent(in) :: q(:), gam real(wp), intent(out) :: fp(:), fm(:) real(wp) :: rho, u, p, c, mach, h_total, mdot_p, mdot_m, p_sp_p, p_sp_m real(wp) :: f_euler(3) call extract_primitives(q, gam, rho, u, p) c = sqrt(gam * p / rho) mach = u / c h_total = (q(3) + p) / rho ! F^+ if (mach >= 1.0_wp) then ! fully supersonic rightward: F+ = F, F- = 0 f_euler(1) = q(2); f_euler(2) = rho * u**2 + p; f_euler(3) = u * (q(3) + p) fp = f_euler fm = 0.0_wp return else if (mach <= -1.0_wp) then ! fully supersonic leftward: F+ = 0, F- = F f_euler(1) = q(2); f_euler(2) = rho * u**2 + p; f_euler(3) = u * (q(3) + p) fp = 0.0_wp fm = f_euler return end if ! Subsonic: compute both splits without branching on sign mdot_p = rho * c * (mach + 1.0_wp)**2 / 4.0_wp p_sp_p = p * (mach + 1.0_wp)**2 * (2.0_wp - mach) / 4.0_wp fp(1) = mdot_p fp(2) = mdot_p * u + p_sp_p fp(3) = mdot_p * h_total mdot_m = -rho * c * (mach - 1.0_wp)**2 / 4.0_wp p_sp_m = p * (mach - 1.0_wp)**2 * (2.0_wp + mach) / 4.0_wp fm(1) = mdot_m fm(2) = mdot_m * u + p_sp_m fm(3) = mdot_m * h_total end subroutine van_leer_split_both ! --------------------------------------------------------------------------- !> Compute the right eigenvector matrix K and its inverse K^{-1} !! of the Euler Jacobian at a given state. !! !! The eigenvectors are used for the characteristic projection in !! the WENO reconstruction of systems. !! See course notes, Eq. (10) for eigenvalues, Eq. (11) for K, !! and Eq. (14) for K^{-1}. !! !! @param[in] q_avg Averaged conserved state (e.g. arithmetic average) !! @param[out] r_mat Right eigenvector matrix K (3x3) !! @param[out] r_inv Left eigenvector matrix K^{-1} (3x3) !! @param[in] gam Ratio of specific heats ! --------------------------------------------------------------------------- pure subroutine compute_eigensystem(q_avg, r_mat, r_inv, gam) use solver_state, only: neq real(wp), intent(in) :: q_avg(neq) real(wp), intent(out) :: r_mat(neq, neq) real(wp), intent(out) :: r_inv(neq, neq) real(wp), intent(in) :: gam real(wp) :: rho, u, p, c, h_total real(wp) :: inv_2c2 ! (γ-1)/(2c²): factor appearing in every row of K^{-1} call extract_primitives(q_avg, gam, rho, u, p) c = sqrt(gam * p / rho) ! sound speed, a h_total = (q_avg(3) + p) / rho ! total specific enthalpy, H = (E+p)/rho ! -- Right eigenvector matrix K, Eq. (11) -- ! ! K^(1) K^(2) K^(3) ! [ 1 1 1 ] ! K = [ u-c u u+c ] ! [ H-u*c 0.5*u^2 H+u*c ] r_mat(1, 1) = 1.0_wp r_mat(2, 1) = u - c r_mat(3, 1) = h_total - u * c r_mat(1, 2) = 1.0_wp r_mat(2, 2) = u r_mat(3, 2) = 0.5_wp * u**2 r_mat(1, 3) = 1.0_wp r_mat(2, 3) = u + c r_mat(3, 3) = h_total + u * c ! -- Left eigenvector matrix K^{-1}, Eq. (14) -- ! Closed-form inverse avoids calling a general LU solver. inv_2c2 = (gam - 1.0_wp) / (2.0_wp * c**2) r_inv(1, 1) = 0.5_wp * (inv_2c2 * u**2 + u / c) r_inv(1, 2) = -0.5_wp * (2.0_wp * inv_2c2 * u + 1.0_wp / c) r_inv(1, 3) = inv_2c2 r_inv(2, 1) = 1.0_wp - inv_2c2 * u**2 r_inv(2, 2) = 2.0_wp * inv_2c2 * u r_inv(2, 3) = -2.0_wp * inv_2c2 r_inv(3, 1) = 0.5_wp * (inv_2c2 * u**2 - u / c) r_inv(3, 2) = -0.5_wp * (2.0_wp * inv_2c2 * u - 1.0_wp / c) r_inv(3, 3) = inv_2c2 end subroutine compute_eigensystem ! --------------------------------------------------------------------------- !> Steger-Warming flux vector splitting. !! !! F^{±} = rho/(2*gamma) * A^{±}, where !! A^{±}₁ = 2(γ-1) λ₂± + λ₁± + λ₃± !! A^{±}₂ = 2(γ-1) u λ₂± + (u-c) λ₁± + (u+c) λ₃± !! A^{±}₃ = (γ-1) u² λ₂± + (H-uc) λ₁± + (H+uc) λ₃± !! and λᵢ± = (λᵢ ± |λᵢ|) / 2 are the split eigenvalues. !! See Steger & Warming (1981), J. Comput. Phys. 40, 263–293. !! !! @param[in] q Conserved state vector !! @param[out] flux Split flux (positive or negative part) !! @param[in] sign +1 for F^+, -1 for F^- !! @param[in] gam Ratio of specific heats ! --------------------------------------------------------------------------- pure subroutine steger_warming_split(q, flux, sign, gam) real(wp), intent(in) :: q(:) real(wp), intent(out) :: flux(:) integer, intent(in) :: sign real(wp), intent(in) :: gam real(wp) :: rho, u, p, c, h_total real(wp) :: lam1, lam2, lam3 ! eigenvalues: u-c, u, u+c real(wp) :: lp1, lp2, lp3 ! split eigenvalues λᵢ± call extract_primitives(q, gam, rho, u, p) c = sqrt(gam * p / rho) ! speed of sound h_total = (q(3) + p) / rho ! total specific enthalpy H lam1 = u - c lam2 = u lam3 = u + c ! λᵢ± = (λᵢ ± |λᵢ|) / 2 lp1 = 0.5_wp * (lam1 + real(sign, wp) * abs(lam1)) lp2 = 0.5_wp * (lam2 + real(sign, wp) * abs(lam2)) lp3 = 0.5_wp * (lam3 + real(sign, wp) * abs(lam3)) flux(1) = (rho / (2.0_wp * gam)) * & (2.0_wp * (gam - 1.0_wp) * lp2 + lp1 + lp3) flux(2) = (rho / (2.0_wp * gam)) * & (2.0_wp * (gam - 1.0_wp) * u * lp2 + (u - c) * lp1 + (u + c) * lp3) flux(3) = (rho / (2.0_wp * gam)) * & ((gam - 1.0_wp) * u**2 * lp2 + (h_total - u * c) * lp1 + (h_total + u * c) * lp3) end subroutine steger_warming_split ! --------------------------------------------------------------------------- !> van Leer flux vector splitting. !! !! For subsonic flow |M| < 1 (M = u/c): !! ṁ± = ±ρc (M±1)² / 4 !! p± = p (1±M)² (2∓M) / 4 !! F± = ( ṁ±, ṁ± u + p±, ṁ± H ) !! For |M| ≥ 1 the supersonic limit applies: F+ = F (M>1), F- = 0 (M>1). !! See van Leer (1982), ICASE Report 82-30. !! !! @param[in] q Conserved state vector !! @param[out] flux Split flux (positive or negative part) !! @param[in] sign +1 for F^+, -1 for F^- !! @param[in] gam Ratio of specific heats ! --------------------------------------------------------------------------- pure subroutine van_leer_split(q, flux, sign, gam) real(wp), intent(in) :: q(:) real(wp), intent(out) :: flux(:) integer, intent(in) :: sign real(wp), intent(in) :: gam real(wp) :: rho, u, p, c, mach, h_total, mdot, p_sp call extract_primitives(q, gam, rho, u, p) c = sqrt(gam * p / rho) mach = u / c h_total = (q(3) + p) / rho ! = c²/(γ-1) + u²/2 if (sign == 1) then ! --- F^+ --- if (mach >= 1.0_wp) then ! fully supersonic rightward: F+ = F flux(1) = q(2) flux(2) = rho * u**2 + p flux(3) = u * (q(3) + p) return else if (mach <= -1.0_wp) then ! fully supersonic leftward: F+ = 0 flux = 0.0_wp return end if mdot = rho * c * (mach + 1.0_wp)**2 / 4.0_wp p_sp = p * (mach + 1.0_wp)**2 * (2.0_wp - mach) / 4.0_wp else ! --- F^- --- if (mach <= -1.0_wp) then ! fully supersonic leftward: F- = F flux(1) = q(2) flux(2) = rho * u**2 + p flux(3) = u * (q(3) + p) return else if (mach >= 1.0_wp) then ! fully supersonic rightward: F- = 0 flux = 0.0_wp return end if mdot = -rho * c * (mach - 1.0_wp)**2 / 4.0_wp p_sp = p * (mach - 1.0_wp)**2 * (2.0_wp + mach) / 4.0_wp end if flux(1) = mdot flux(2) = mdot * u + p_sp flux(3) = mdot * h_total end subroutine van_leer_split ! --------------------------------------------------------------------------- !> AUSM+ FDS scheme (AUSM family). !! !! Computes the numerical flux at a cell interface given the reconstructed !! left and right conserved states Q_L and Q_R. !! !! Interface sound speed (Liou 1996): !! c* = sqrt(2(γ-1)/(γ+1) · H), c_½ = min(c*²_L/max(c*_L,|u_L|), !! c*²_R/max(c*_R,|u_R|)) !! !! Split Mach polynomials at Mᵢ = uᵢ/c_½ (i = L, R): !! M± = ±(M±1)²/4 for |M|≤1, (M±|M|)/2 otherwise !! !! Face flux: F_½ = ṁ φ_{upwind} + [0, p_½, 0]ᵀ !! ṁ = c_½ (M⁺(M_L) ρ_L + M⁻(M_R) ρ_R) !! p_½ = P⁺(M_L) p_L + P⁻(M_R) p_R !! φ = [1, u, H]ᵀ, selected from L or R by sign(ṁ) !! See Liou (1996), J. Comput. Phys. 129, 364–382. !! !! @param[in] q_L Left reconstructed conserved state !! @param[in] q_R Right reconstructed conserved state !! @param[out] flux Numerical face flux !! @param[in] gam Ratio of specific heats ! --------------------------------------------------------------------------- pure subroutine ausm_plus(q_L, q_R, flux, gam) real(wp), intent(in) :: q_L(:), q_R(:) real(wp), intent(out) :: flux(:) real(wp), intent(in) :: gam real(wp) :: rho_L, u_L, p_L, h_L real(wp) :: rho_R, u_R, p_R, h_R real(wp) :: c_s_L, c_s_R, c_half real(wp) :: mach_L, mach_R real(wp) :: mplus, mminus, mdot, p_face real(wp) :: pplus, pminus ! Decode left and right states call extract_primitives(q_L, gam, rho_L, u_L, p_L) h_L = (q_L(3) + p_L) / rho_L ! total specific enthalpy call extract_primitives(q_R, gam, rho_R, u_R, p_R) h_R = (q_R(3) + p_R) / rho_R ! Critical sound speed: c*² = 2(γ-1)/(γ+1) · H ! Interface sound speed: c_½ = min(c*²_L/max(c*_L,|u_L|), c*²_R/max(c*_R,|u_R|)) c_s_L = sqrt(2.0_wp * (gam - 1.0_wp) / (gam + 1.0_wp) * h_L) c_s_R = sqrt(2.0_wp * (gam - 1.0_wp) / (gam + 1.0_wp) * h_R) c_half = min(c_s_L**2 / max(c_s_L, abs(u_L)), & c_s_R**2 / max(c_s_R, abs(u_R))) mach_L = u_L / c_half mach_R = u_R / c_half ! M+ from left state if (abs(mach_L) <= 1.0_wp) then mplus = 0.25_wp * (mach_L + 1.0_wp)**2 else mplus = 0.5_wp * (mach_L + abs(mach_L)) end if ! M- from right state if (abs(mach_R) <= 1.0_wp) then mminus = -0.25_wp * (mach_R - 1.0_wp)**2 else mminus = 0.5_wp * (mach_R - abs(mach_R)) end if ! Interface mass flux: ṁ = c_½ (M⁺ ρ_L + M⁻ ρ_R) mdot = c_half * (mplus * rho_L + mminus * rho_R) ! P+ from left state if (abs(mach_L) <= 1.0_wp) then pplus = 0.25_wp * p_L * (mach_L + 1.0_wp)**2 * (2.0_wp - mach_L) else if (mach_L > 1.0_wp) then pplus = p_L else pplus = 0.0_wp end if ! P- from right state if (abs(mach_R) <= 1.0_wp) then pminus = 0.25_wp * p_R * (mach_R - 1.0_wp)**2 * (2.0_wp + mach_R) else if (mach_R < -1.0_wp) then pminus = p_R else pminus = 0.0_wp end if p_face = pplus + pminus ! Convective flux: ṁ φ upwinded by sign of ṁ if (mdot >= 0.0_wp) then flux(1) = mdot flux(2) = mdot * u_L + p_face flux(3) = mdot * h_L else flux(1) = mdot flux(2) = mdot * u_R + p_face flux(3) = mdot * h_R end if end subroutine ausm_plus ! --------------------------------------------------------------------------- !> HLL approximate Riemann solver. !! !! Uses Einfeldt (HLLE) wave-speed estimates based on the minimum/maximum of !! Davis and Roe-averaged wave speeds, providing robust positivity preservation. !! !! Wave speeds: !! S_L = min(u_L - c_L, ũ - c̃), S_R = max(u_R + c_R, ũ + c̃) !! !! Numerical flux: !! F_HLL = F_L if S_L ≥ 0 !! F_HLL = F_R if S_R ≤ 0 !! F_HLL = (S_R F_L - S_L F_R + S_L S_R (Q_R - Q_L)) / (S_R - S_L) otherwise !! !! See Harten, Lax & van Leer (1983), SIAM Rev. 25, 35–61. !! See Einfeldt (1988), SIAM J. Numer. Anal. 25, 294–318. !! !! @param[in] q_L Left reconstructed conserved state !! @param[in] q_R Right reconstructed conserved state !! @param[out] flux Numerical face flux !! @param[in] gam Ratio of specific heats ! --------------------------------------------------------------------------- pure subroutine hll_flux(q_L, q_R, flux, gam) real(wp), intent(in) :: q_L(:), q_R(:) real(wp), intent(out) :: flux(:) real(wp), intent(in) :: gam real(wp) :: rho_L, u_L, p_L, c_L real(wp) :: rho_R, u_R, p_R, c_R real(wp) :: f_L(3), f_R(3) real(wp) :: u_roe, h_L, h_R, h_roe, c_roe real(wp) :: s_L, s_R ! Decode left and right states call extract_primitives(q_L, gam, rho_L, u_L, p_L) c_L = sqrt(gam * p_L / rho_L) h_L = (q_L(3) + p_L) / rho_L call extract_primitives(q_R, gam, rho_R, u_R, p_R) c_R = sqrt(gam * p_R / rho_R) h_R = (q_R(3) + p_R) / rho_R ! Roe-averaged velocity, enthalpy, and sound speed (√ρ weighting) call compute_roe_averages(rho_L, u_L, h_L, rho_R, u_R, h_R, gam, u_roe, h_roe, c_roe) ! Einfeldt wave-speed estimates s_L = min(u_L - c_L, u_roe - c_roe) s_R = max(u_R + c_R, u_roe + c_roe) ! Euler fluxes at left and right states call compute_euler_flux(q_L, f_L, gam) call compute_euler_flux(q_R, f_R, gam) ! HLL flux selection if (s_L >= 0.0_wp) then flux = f_L else if (s_R <= 0.0_wp) then flux = f_R else flux = (s_R * f_L - s_L * f_R + s_L * s_R * (q_R - q_L)) / (s_R - s_L) end if end subroutine hll_flux ! --------------------------------------------------------------------------- !> HLLC approximate Riemann solver with contact-wave restoration. !! !! Extends HLL by including the contact/shear wave at speed S*, yielding !! sharper resolution of contact discontinuities and shear layers. !! !! Uses the same Einfeldt wave-speed estimates as HLL: !! S_L = min(u_L - c_L, ũ - c̃), S_R = max(u_R + c_R, ũ + c̃) !! !! Contact speed: !! S* = (p_R - p_L + ρ_L u_L (S_L - u_L) - ρ_R u_R (S_R - u_R)) !! / (ρ_L (S_L - u_L) - ρ_R (S_R - u_R)) !! !! Intermediate state (K = L or R): !! U*_K = ρ_K (S_K - u_K)/(S_K - S*) · [1, S*, E_K/ρ_K + (S*-u_K)(S* + p_K/(ρ_K(S_K-u_K)))] !! !! Numerical flux: !! F_L if S_L ≥ 0 !! F_L + S_L (U*_L - Q_L) if S_L ≤ 0 ≤ S* !! F_R + S_R (U*_R - Q_R) if S* ≤ 0 ≤ S_R !! F_R if S_R ≤ 0 !! !! See Toro, Spruce & Speares (1994), Shock Waves 4, 25–34. !! !! @param[in] q_L Left reconstructed conserved state !! @param[in] q_R Right reconstructed conserved state !! @param[out] flux Numerical face flux !! @param[in] gam Ratio of specific heats ! --------------------------------------------------------------------------- pure subroutine hllc_flux(q_L, q_R, flux, gam) real(wp), intent(in) :: q_L(:), q_R(:) real(wp), intent(out) :: flux(:) real(wp), intent(in) :: gam real(wp) :: rho_L, u_L, p_L, c_L, e_L real(wp) :: rho_R, u_R, p_R, c_R, e_R real(wp) :: f_L(3), f_R(3) real(wp) :: u_roe, h_L, h_R, h_roe, c_roe real(wp) :: s_L, s_R, s_star real(wp) :: factor, u_star(3) ! Decode left and right states call extract_primitives(q_L, gam, rho_L, u_L, p_L) c_L = sqrt(gam * p_L / rho_L) e_L = q_L(3) / rho_L ! specific total energy E/ρ h_L = e_L + p_L / rho_L call extract_primitives(q_R, gam, rho_R, u_R, p_R) c_R = sqrt(gam * p_R / rho_R) e_R = q_R(3) / rho_R h_R = e_R + p_R / rho_R ! Roe-averaged velocity, enthalpy, and sound speed (√ρ weighting) call compute_roe_averages(rho_L, u_L, h_L, rho_R, u_R, h_R, gam, u_roe, h_roe, c_roe) ! Einfeldt wave-speed estimates s_L = min(u_L - c_L, u_roe - c_roe) s_R = max(u_R + c_R, u_roe + c_roe) ! Euler fluxes at left and right states call compute_euler_flux(q_L, f_L, gam) call compute_euler_flux(q_R, f_R, gam) if (s_L >= 0.0_wp) then flux = f_L return end if if (s_R <= 0.0_wp) then flux = f_R return end if ! Contact wave speed S* (Toro 1994, Eq. 10) s_star = (p_R - p_L & + rho_L * u_L * (s_L - u_L) & - rho_R * u_R * (s_R - u_R)) & / (rho_L * (s_L - u_L) - rho_R * (s_R - u_R)) if (s_star >= 0.0_wp) then ! Left star state: F*_L = F_L + S_L (U*_L - Q_L) factor = rho_L * (s_L - u_L) / (s_L - s_star) u_star(1) = factor u_star(2) = factor * s_star u_star(3) = factor * (e_L + (s_star - u_L) * (s_star + p_L / (rho_L * (s_L - u_L)))) flux = f_L + s_L * (u_star - q_L) else ! Right star state: F*_R = F_R + S_R (U*_R - Q_R) factor = rho_R * (s_R - u_R) / (s_R - s_star) u_star(1) = factor u_star(2) = factor * s_star u_star(3) = factor * (e_R + (s_star - u_R) * (s_star + p_R / (rho_R * (s_R - u_R)))) flux = f_R + s_R * (u_star - q_R) end if end subroutine hllc_flux ! --------------------------------------------------------------------------- !> Roe approximate Riemann solver with Harten entropy fix. !! !! Computes the face flux using Roe-averaged states and characteristic !! wave decomposition. An entropy fix prevents non-physical expansion shocks !! at sonic points by smoothing acoustic eigenvalues that change sign. !! !! Roe averages (√ρ-weighted): !! ũ = (√ρ_L u_L + √ρ_R u_R) / (√ρ_L + √ρ_R) !! H̃ = (√ρ_L H_L + √ρ_R H_R) / (√ρ_L + √ρ_R) !! c̃ = √((γ-1)(H̃ - ũ²/2)), ρ̃ = √(ρ_L ρ_R) !! !! Harten entropy fix for acoustic eigenvalues λ₁ = ũ-c̃, λ₃ = ũ+c̃ !! (Harten & Hyman 1983, J. Comput. Phys. 50): !! δ = 2 ε c̃ (ε = 0.1) !! |λ|_fix = (λ² + δ²)/(2δ) if |λ| < δ, else |λ| !! !! Roe flux: !! F_Roe = ½(F_L + F_R) - ½ Σₖ |λₖ|_fix Δwₖ rₖ !! !! See Roe (1981), J. Comput. Phys. 43, 357–372. !! See Harten & Hyman (1983), J. Comput. Phys. 50, 235–269. !! !! @param[in] q_L Left reconstructed conserved state !! @param[in] q_R Right reconstructed conserved state !! @param[out] flux Numerical face flux !! @param[in] gam Ratio of specific heats ! --------------------------------------------------------------------------- pure subroutine roe_flux(q_L, q_R, flux, gam) real(wp), intent(in) :: q_L(:), q_R(:) real(wp), intent(out) :: flux(:) real(wp), intent(in) :: gam real(wp) :: rho_L, u_L, p_L, h_L real(wp) :: rho_R, u_R, p_R, h_R real(wp) :: f_L(3), f_R(3) real(wp) :: rho_roe, u_roe, h_roe, c_roe real(wp) :: lam1, lam2, lam3 ! eigenvalues: ũ-c̃, ũ, ũ+c̃ real(wp) :: abs_lam1, abs_lam3 ! entropy-fixed |λ| real(wp) :: delta ! entropy-fix threshold real(wp) :: dp, du, drho real(wp) :: dw1, dw2, dw3 ! wave strengths ! Entropy-fix width δ = 2ε·c̃; ε=0.1 is O(Δx) and tuned empirically. ! See Harten & Hyman (1983), J. Comput. Phys. 50, Eq. (2.11). real(wp), parameter :: eps_entropy = 0.1_wp ! Decode left and right states call extract_primitives(q_L, gam, rho_L, u_L, p_L) h_L = (q_L(3) + p_L) / rho_L call extract_primitives(q_R, gam, rho_R, u_R, p_R) h_R = (q_R(3) + p_R) / rho_R ! Roe averages (√ρ weighting); ρ̃ = √(ρ_L ρ_R) call compute_roe_averages(rho_L, u_L, h_L, rho_R, u_R, h_R, gam, u_roe, h_roe, c_roe) rho_roe = sqrt(rho_L * rho_R) ! Eigenvalues lam1 = u_roe - c_roe lam2 = u_roe lam3 = u_roe + c_roe ! Harten entropy fix for acoustic waves (λ₁, λ₃) ! See Harten & Hyman (1983), Eq. (2.11) delta = 2.0_wp * eps_entropy * c_roe if (abs(lam1) < delta) then abs_lam1 = (lam1**2 + delta**2) / (2.0_wp * delta) else abs_lam1 = abs(lam1) end if if (abs(lam3) < delta) then abs_lam3 = (lam3**2 + delta**2) / (2.0_wp * delta) else abs_lam3 = abs(lam3) end if ! Primitive differences dp = p_R - p_L du = u_R - u_L drho = rho_R - rho_L ! Wave strengths Δw (differences in characteristic variables) ! Derived from K⁻¹ ΔQ where K is the Roe-average right-eigenvector matrix dw1 = (dp - rho_roe * c_roe * du) / (2.0_wp * c_roe**2) dw2 = drho - dp / c_roe**2 dw3 = (dp + rho_roe * c_roe * du) / (2.0_wp * c_roe**2) ! Euler fluxes at left and right states call compute_euler_flux(q_L, f_L, gam) call compute_euler_flux(q_R, f_R, gam) ! Roe flux: F = ½(F_L + F_R) - ½ Σ |λₖ| Δwₖ rₖ ! Right eigenvectors evaluated at Roe-averaged state flux(1) = 0.5_wp * (f_L(1) + f_R(1)) & - 0.5_wp * (abs_lam1 * dw1 * 1.0_wp & + abs(lam2) * dw2 * 1.0_wp & + abs_lam3 * dw3 * 1.0_wp) flux(2) = 0.5_wp * (f_L(2) + f_R(2)) & - 0.5_wp * (abs_lam1 * dw1 * (u_roe - c_roe) & + abs(lam2) * dw2 * u_roe & + abs_lam3 * dw3 * (u_roe + c_roe)) flux(3) = 0.5_wp * (f_L(3) + f_R(3)) & - 0.5_wp * (abs_lam1 * dw1 * (h_roe - u_roe * c_roe) & + abs(lam2) * dw2 * (0.5_wp * u_roe**2) & + abs_lam3 * dw3 * (h_roe + u_roe * c_roe)) end subroutine roe_flux ! --------------------------------------------------------------------------- !> Return the global maximum wave speed max_i(|u_i| + c_i). !! !! Scans all grid points and returns the largest signal speed, which is !! used by the driver to compute a CFL-stable time step: !! dt = cfl * dx / compute_max_wave_speed(state) !! !! @param[in] state Solver state (reads state%ub, state%n_pt, state%cfg%gam) !! @return s_max Maximum wave speed over all grid points [m/s] ! --------------------------------------------------------------------------- function compute_max_wave_speed(state) result(s_max) use solver_state, only: solver_state_t type(solver_state_t), intent(in) :: state real(wp) :: s_max integer :: i real(wp) :: rho, u, p, c s_max = 0.0_wp do i = 1, state % n_pt call extract_primitives(state % ub(:, i), state % cfg % gam, rho, u, p) c = sqrt(state % cfg % gam * p / rho) s_max = max(s_max, abs(u) + c) end do end function compute_max_wave_speed end module euler_physics