diff --git a/docs/documentation/equations.md b/docs/documentation/equations.md index 0c120da36f..4db6899a40 100644 --- a/docs/documentation/equations.md +++ b/docs/documentation/equations.md @@ -359,6 +359,26 @@ Temperature is obtained from the internal energy by Newton iteration: \f[\frac{\hat{h}_m}{R_u\,T} = \frac{C_0}{T} + \sum_{r=1}^{5} \frac{C_r\,T^{r-1}}{r}\f] +### 3.3 Thermodynamics Interface + +Every solver path obtains thermodynamic quantities through one small set of operations rather than through EOS arithmetic repeated at each call site. The operations, with their realization in `src/common/m_variables_conversion.fpp` unless noted: + +- **Pressure from total energy** \f$p(\rho E, \ldots)\f$: `s_compute_pressure()`. +- **Total energy from pressure** \f$\rho E(p, \ldots)\f$: reconstructed in `s_convert_primitive_to_conservative_variables()`, and in the Riemann solvers for the wave speeds. +- **Speed of sound** \f$c\f$: `s_compute_speed_of_sound()`. +- **Temperature** \f$T\f$: `get_temperature()`, use-associated from the Pyrometheus-generated `m_thermochem.f90` (chemistry only; the stiffened-gas closure is temperature-free). +- **Mixture thermodynamic parameters** \f$(\Gamma, \Pi_\infty, q_v)\f$: `s_convert_to_mixture_variables()`. +- **State admissibility**: `f_validate_state()`, the predicate that gates the sound-speed floor and is shared by the simulation and post-process paths. + +Two adapters supply these operations, selected at compile time: + +- **Stiffened gas** (Section 3.1) is the closure for `chemistry = .false.`. +- **Ideal-gas mixture** (Section 3.2), generated by Pyrometheus, is selected automatically when `chemistry = .true.`. + +Selection is a compile-time branch on the `chemistry` parameter (`#:if not chemistry` versus `#:else`), not a runtime dispatch, so the device routines carry no indirection and the CPU-only and GPU builds resolve the same code. Adding an EOS family is a matter of adding a branch to these operations, not editing every call site. + +Extending `f_validate_state()` follows a fixed policy: per-cell scalars enter as arguments at the rung that needs them, EOS parameters and configuration flow in through module state, and each call site keeps its own floor. + --- ## 4. Viscous Stress Tensor (`viscous = .true.`) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index a230416b5e..4792bb4d8d 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -35,6 +35,7 @@ module m_variables_conversion s_compute_species_fraction, & #ifndef MFC_PRE_PROCESS s_compute_speed_of_sound, & + f_validate_state, & s_compute_fast_magnetosonic_speed, & #endif s_finalize_variables_conversion_module @@ -1286,14 +1287,28 @@ contains c = (H - 5.e-1*vel_sum - qv/rho)/gamma end if - if (mixture_err .and. c < 0._wp) then - c = 100._wp*sgm_eps - else + if (f_validate_state(c)) then c = sqrt(c) + else + c = 100._wp*sgm_eps end if end if end subroutine s_compute_speed_of_sound + + !> Admissibility predicate of the thermodynamics interface: true when the squared sound speed yields a real value. With + !! mixture_err off it always passes, so a negative argument reaches sqrt and returns NaN exactly as the pre-interface code did. + !! Shared by the simulation and post-process sound-speed paths, which each keep their own floor. + pure function f_validate_state(c_squared) result(admissible) + + $:GPU_ROUTINE(function_name='f_validate_state', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: c_squared + logical :: admissible + + admissible = .not. (mixture_err .and. c_squared < 0._wp) + + end function f_validate_state #endif #ifndef MFC_PRE_PROCESS diff --git a/src/post_process/m_derived_variables.fpp b/src/post_process/m_derived_variables.fpp index 05f1fc1805..521c701979 100644 --- a/src/post_process/m_derived_variables.fpp +++ b/src/post_process/m_derived_variables.fpp @@ -111,10 +111,10 @@ contains & k)/blkmod1 + (1._wp - q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k))/blkmod2))) end if - if (mixture_err .and. q_sf(i, j, k) < 0._wp) then - q_sf(i, j, k) = 1.e-16_wp - else + if (f_validate_state(q_sf(i, j, k))) then q_sf(i, j, k) = sqrt(q_sf(i, j, k)) + else + q_sf(i, j, k) = 1.e-16_wp end if end do end do