Refactor: establish the thermodynamics interface

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.`)

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

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