diff --git a/cmake/ParamsCodegen.cmake b/cmake/ParamsCodegen.cmake index 98d9362d4..a27b08188 100644 --- a/cmake/ParamsCodegen.cmake +++ b/cmake/ParamsCodegen.cmake @@ -9,7 +9,7 @@ file(GLOB_RECURSE _mfc_gen_inputs "${CMAKE_CURRENT_SOURCE_DIR}/toolchain/mfc/params/*.py" ) -# Enumerate the 15 generated .fpp files explicitly so ninja can track them as +# Enumerate the 18 generated .fpp files explicitly so ninja can track them as # build-time outputs and so HANDLE_SOURCES does not need a configure-time GLOB # of ${CMAKE_BINARY_DIR}/include// (which fails when the dir is empty). set(_mfc_gen_inc "${CMAKE_BINARY_DIR}/include") @@ -17,6 +17,7 @@ set(_mfc_gen_files_pre_process "${_mfc_gen_inc}/pre_process/generated_namelist.fpp" "${_mfc_gen_inc}/pre_process/generated_decls.fpp" "${_mfc_gen_inc}/pre_process/generated_constants.fpp" + "${_mfc_gen_inc}/pre_process/generated_eos.fpp" "${_mfc_gen_inc}/pre_process/generated_bcast.fpp" "${_mfc_gen_inc}/pre_process/generated_case_opt_decls.fpp" ) @@ -24,6 +25,7 @@ set(_mfc_gen_files_simulation "${_mfc_gen_inc}/simulation/generated_namelist.fpp" "${_mfc_gen_inc}/simulation/generated_decls.fpp" "${_mfc_gen_inc}/simulation/generated_constants.fpp" + "${_mfc_gen_inc}/simulation/generated_eos.fpp" "${_mfc_gen_inc}/simulation/generated_bcast.fpp" "${_mfc_gen_inc}/simulation/generated_case_opt_decls.fpp" ) @@ -31,6 +33,7 @@ set(_mfc_gen_files_post_process "${_mfc_gen_inc}/post_process/generated_namelist.fpp" "${_mfc_gen_inc}/post_process/generated_decls.fpp" "${_mfc_gen_inc}/post_process/generated_constants.fpp" + "${_mfc_gen_inc}/post_process/generated_eos.fpp" "${_mfc_gen_inc}/post_process/generated_bcast.fpp" "${_mfc_gen_inc}/post_process/generated_case_opt_decls.fpp" ) diff --git a/docs/documentation/contributing.md b/docs/documentation/contributing.md index dd3c0f33c..1db7d18b6 100644 --- a/docs/documentation/contributing.md +++ b/docs/documentation/contributing.md @@ -87,7 +87,7 @@ between them is narrow. mfc.sh (env bootstrap, venv, module loading, lock) └─ toolchain/mfc/build.py (config slugs, cmake invocation) └─ CMakeLists.txt + cmake/{GPU,Fypp,ParamsCodegen,MFCTargets}.cmake - └─ toolchain/mfc/params/generators/cmake_gen.py (writes 15 generated .fpp includes) + └─ toolchain/mfc/params/generators/cmake_gen.py (writes 18 generated .fpp includes) ``` **`mfc.sh` → `build.py`.** `mfc.sh` is a thin shell wrapper that activates the Python @@ -100,8 +100,8 @@ mode, debug, chemistry, MPI. Staging and install trees are namespaced by slug u **CMake layer.** `cmake/Fypp.cmake` defines `HANDLE_SOURCES`, which sets up one `add_custom_command` per `.fpp` file to run Fypp at build time. `cmake/ParamsCodegen.cmake` registers a single ninja-tracked `add_custom_command` (DEPENDS all `params/*.py`) that -invokes `cmake_gen.py` and writes the 15 generated includes under -`build/include//`. There is no configure-time generation: all 15 files are build +invokes `cmake_gen.py` and writes the 18 generated includes under +`build/include//`. There is no configure-time generation: all 18 files are build outputs, so changing any `params/*.py` triggers only a targeted rebuild, not a full reconfigure. @@ -466,11 +466,14 @@ If an array is allocated inside an `if` block, its deallocation must follow the ### How to Add an Equation of State -Every stiffened-gas expression lives in `src/common/m_variables_conversion.fpp`. Adding a second EOS -means supplying these, not grepping for `gammas`: +The equation-of-state operators live in `src/common/m_eos.fpp`; the mixture closure rules that +combine them (`s_compute_mixture_coefficients`, `s_compute_speed_of_sound` and their variants) +stay in `src/common/m_variables_conversion.fpp`. Adding a second EOS means supplying these, not +grepping for `gammas`: | Operator | Gives | |---|---| +| `s_reference_curve` | the reference curve \f$p_{ref}, e_{ref}\f$ and \f$\Gamma_G\f$ of a state-dependent family - one `case` per family, and nothing else | | `s_compute_mixture_coefficients` / `_dt` | mixture \f$\Gamma, \Pi_\infty, q_v\f$ from the phase fractions, and their time derivative | | `f_pressure` / `s_compute_energy` | \f$p(e)\f$ and \f$E(p)\f$ | | `f_bulk_modulus` | \f$K(p)\f$ - every sound speed in MFC is \f$K/\rho\f$, differing only in how phases are mixed | @@ -489,6 +492,35 @@ forms the user supplies (see @ref sec-stored-forms), and `isentrope_n`/`isentrop as \f$p + B = \textrm{const}\,\rho^n\f$, derived once at start-up. Convert with the `f_isentrope_*` operators rather than open-coding either relation. +Both are resolved once in `s_initialize_eos_module`. A state-dependent family skips this and computes +its coefficients per cell from `s_reference_curve` instead. + +**Adding a family** - a new state-dependent parameter set alongside Mie-Gruneisen, JWL and Vinet - +starts at `toolchain/mfc/params/eos_families.py`: one `EosFamily` entry. + +From that entry, these are generated and need no hand edit: the Fortran `eos_*` constants, the +`physical_parameters` parameter registration in `definitions.py`, the validator's per-family +required/optional parameter sets and its initial-state coefficient call, both the +`f_is_state_dependent` and `f_has_isentropic_reference` family tests, both layers of +`s_initialize_eos_module`, and the `any_state_dependent_eos` case-optimization flag. + +Three things stay hand-written alongside the entry, each checked against the registry by a test +that fails on disagreement: the `case` body in `s_reference_curve` (the per-family mathematics), +the mirror function in `toolchain/mfc/eos.py`, and the family's fields on `physical_parameters` in +`src/common/m_derived_types.fpp`. Those four edits are the whole job in the normal case. + +Two further edits are possible and are guarded by the compiler rather than by a test - they fail +loudly at build with a clear error, not silently at runtime: + +- `m_eos.fpp`'s `use m_constants, only: ...` list must name the new `eos_` constant. The + generated family predicates reference it, so omitting it fails to compile. +- the `eos_coefficients` type in `src/common/m_derived_types.fpp` needs a new field only if the + family stores a coefficient the existing ones do not provide. The generated assignment in + `s_initialize_eos_module` fails to compile if the field is missing. + +One caveat: `f_has_isentropic_reference`'s `gruneisen_a == 0._wp` conjunct is a runtime test, not a +family property, and stays hand-written in `m_eos.fpp`. + ### How to Add a Test Case **Step 1: Create a case file** diff --git a/docs/documentation/equations.md b/docs/documentation/equations.md index aaba7934a..21edbb533 100644 --- a/docs/documentation/equations.md +++ b/docs/documentation/equations.md @@ -25,7 +25,7 @@ where: The parameter `model_eqns` (1, 2, or 3) selects the governing equation set. -**Key source files:** `src/simulation/m_rhs.fpp` (RHS evaluation), `src/common/m_variables_conversion.fpp` (EOS and variable conversion). +**Key source files:** `src/simulation/m_rhs.fpp` (RHS evaluation), `src/common/m_eos.fpp` (equations of state), `src/common/m_variables_conversion.fpp` (variable conversion and mixture rules). --- diff --git a/docs/module_categories.json b/docs/module_categories.json index 23a703fa2..f3b0857f2 100644 --- a/docs/module_categories.json +++ b/docs/module_categories.json @@ -20,6 +20,7 @@ { "category": "Physics Models", "modules": [ + "m_eos", "m_viscous", "m_hb_function", "m_surface_tension", diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 3e5eb2fad..d2f84239c 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -113,16 +113,12 @@ module m_constants integer, parameter :: BC_NO_SLIP_WALL = -16 integer, parameter :: BC_DIRICHLET = -17 - ! Synthetic turbulence array size limits - !> Equation of state per fluid. Values must match _EOS_NAMES in toolchain/mfc/params/definitions.py; compound registry keys - !! cannot be auto-generated, so these are hand-written. - integer, parameter :: eos_stiffened_gas = 1 - integer, parameter :: eos_ideal_gas = 2 - integer, parameter :: eos_mie_gruneisen = 3 - integer, parameter :: eos_jwl = 4 - integer, parameter :: eos_vinet = 5 + !> Equation of state per fluid (eos_stiffened_gas, eos_ideal_gas, eos_mie_gruneisen, eos_jwl, eos_vinet) is generated from + !! EOS_FAMILIES in toolchain/mfc/params/eos_families.py; see generated_constants.fpp below. integer, parameter :: eos_rk4_steps = 8 !< fixed-step RK4 along a phasic isentrope or a reference temperature integer, parameter :: ode_isentrope = 1, ode_reference_temperature = 2 !< the two ODEs s_rk4 integrates + + ! Synthetic turbulence array size limits integer, parameter :: num_synth_shells_max = 50 !< Max energy shells for synthetic turbulence integer, parameter :: num_turb_sources_max = 10 !< Max Gaussian forcing zones for synthetic turbulence diff --git a/src/common/m_eos.fpp b/src/common/m_eos.fpp new file mode 100644 index 000000000..8f5a99c28 --- /dev/null +++ b/src/common/m_eos.fpp @@ -0,0 +1,527 @@ +!> +!! @file +!! @brief Contains module m_eos + +#:include 'macros.fpp' +#:include 'case.fpp' +#:include 'generated_eos.fpp' + +!> @brief Equations of state in Gamma/Pi form, rho e = Gamma(rho) p + Pi(rho). +!! +!! Stiffened and ideal gas keep constant coefficients, resolved once at start-up. The +!! state-dependent families (Mie-Gruneisen, JWL, Vinet) evaluate theirs per cell from a +!! reference curve. Mixture closure rules -- Wood's law, the six-equation mean, the bubbly +!! branch -- are not equations of state and live in m_variables_conversion. +!! +!! This module is a leaf: it directly uses only m_derived_types, m_constants, and +!! m_global_parameters_common. Adding an EOS family means one case in s_reference_curve. +module m_eos + + use m_derived_types + use m_global_parameters_common + use m_mpi_common, only: s_prohibit_abort + use m_constants, only: eos_ideal_gas, eos_mie_gruneisen, eos_jwl, eos_vinet, eos_rk4_steps, ode_isentrope, & + & ode_reference_temperature, sgm_eps, verysmall, dflt_real + + implicit none + + private + + public :: s_initialize_eos_module, s_finalize_eos_module, f_pressure, f_bulk_modulus, f_relativistic_enthalpy, & + & f_isentrope_exponent, f_isentrope_pressure, f_sg_thermal, f_is_state_dependent, s_phase_coefficients, & + & s_phase_pressure_on_isentrope, s_phase_temperature, s_phase_density_on_isentrope, s_phase_internal_energy, & + & s_phase_bulk_modulus + +contains + + !> Resolve every fluid's EOS coefficients once, before any conversion runs. + impure subroutine s_initialize_eos_module() + + integer :: i + logical :: state_dependent !< Whether this case's fluids need a density-dependent EOS + + @:ALLOCATE(gammas (1:num_fluids)) + @:ALLOCATE(eoss (1:num_fluids)) + @:ALLOCATE(isentrope_n (1:num_fluids)) + @:ALLOCATE(pi_infs(1:num_fluids)) + @:ALLOCATE(isentrope_B(1:num_fluids)) + @:ALLOCATE(cvs (1:num_fluids)) + @:ALLOCATE(qvs (1:num_fluids)) + @:ALLOCATE(qvps (1:num_fluids)) + + state_dependent = .false. + do i = 1, num_fluids + gammas(i) = fluid_pp(i)%gamma + isentrope_n(i) = f_isentrope_exponent(gammas(i)) + + ! Each EOS supplies its own coefficients. Resolved once here, not per cell: a branch in the mixture loop costs + ! registers in the Riemann kernels. An EOS whose coefficients depend on state must move to per-cell evaluation. + select case (fluid_pp(i)%eos) + case (eos_ideal_gas) + pi_infs(i) = 0._wp + case default + pi_infs(i) = fluid_pp(i)%pi_inf + end select + isentrope_B(i) = f_isentrope_pressure(pi_infs(i), gammas(i)) + cvs(i) = fluid_pp(i)%cv + qvs(i) = fluid_pp(i)%qv + qvps(i) = fluid_pp(i)%qvp + eoss(i) = fluid_pp(i)%eos + ! Every fluid's single-source coefficients, mu_max among them: where a cubic Hugoniot fit turns over. + ! mu(u_p) peaks where c0 = s2 u_p^2 + 2 s3 u_p^3, and past it no shock state exists, so the Newton below + ! would wander. Solved once here, on the host. + @:EOS_INIT_COEFFS(i) + ! One reference state and Gruneisen closure for every family; the user-facing names keep their prefix. + @:EOS_INIT_REFERENCE_STATE(i) + if (f_is_state_dependent(i)) state_dependent = .true. + end do + #:if MFC_CASE_OPTIMIZATION + ! Baked in at build time, so a case that changed its EOS family since the build would silently + ! run the wrong branch. The namelist still carries fluid_pp%eos, so check the two agree. + @:PROHIBIT(state_dependent .neqv. any_state_dependent_eos, & + & "This case's equations of state do not match the ones this case-optimized binary was built for. Rebuild.") + #:else + any_state_dependent_eos = state_dependent + #:endif + $:GPU_UPDATE(device='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, eoss, eos_coeffs]') + #:if not MFC_CASE_OPTIMIZATION + $:GPU_UPDATE(device='[any_state_dependent_eos]') + #:endif + + end subroutine s_initialize_eos_module + + !> Deallocate the fluid property arrays allocated in s_initialize_eos_module. + impure subroutine s_finalize_eos_module() + + @:DEALLOCATE(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, eoss) + + end subroutine s_finalize_eos_module + + !> The reference curve of a state-dependent EOS at rho: p_ref, e_ref, their d/drho, and Gamma_G with its d/drho. A new family + !! adds one case here and nothing else. + subroutine s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, G0, dG0) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: rho + integer, intent(in) :: i + real(wp), intent(out) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0 + real(wp) :: mu, d, V, ea, eb, up, us, dus, dup_dmu, x, ex, dp_dmu, de_dmu + integer :: iter + + mu = rho/eos_coeffs(i)%rho0 - 1._wp + ! Past the fit's turnover there is no shock state to find; clamp rather than let the Newton below wander + ! off and return a silently wrong pressure. mu_max is huge for the linear fit, so this is a no-op there. + ! Bounded here so the step finishes and the host-side check in s_write_run_time_information can report it; + ! past the turnover there is no shock state and the Newton below would wander. + if (eoss(i) == eos_mie_gruneisen .and. mu > eos_coeffs(i)%mu_max) mu = eos_coeffs(i)%mu_max + select case (eoss(i)) + case (eos_mie_gruneisen) + ! Hugoniot reference u_s = c0 + s u_p + s2 u_p^2 + s3 u_p^3, with p_H = rho0 u_s u_p and the Hugoniot + ! energy e_H = p_H mu/(2 rho0 (1 + mu)); linear on release. Pole at mu = 1/(s - 1) for the linear fit; + ! the validator refuses initial states outside the EOS. + if (mu < 0._wp) then + p_ref = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2*mu + dp_dmu = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2 + else if (eos_coeffs(i)%s2 == 0._wp .and. eos_coeffs(i)%s3 == 0._wp) then + d = 1._wp - (eos_coeffs(i)%s - 1._wp)*mu + p_ref = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2*mu*(1._wp + mu)/(d*d) + dp_dmu = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2*((1._wp + 2._wp*mu)*d + 2._wp*(eos_coeffs(i)%s - 1._wp)*mu*(1._wp & + & + mu))/(d*d*d) + else + ! u_p solves u_s(u_p) mu = u_p (1 + mu): Newton from the linear fit, then implicit differentiation + up = eos_coeffs(i)%c0*mu/(1._wp - (eos_coeffs(i)%s - 1._wp)*mu) + $:GPU_LOOP(parallelism='[seq]') + do iter = 1, 8 + us = eos_coeffs(i)%c0 + up*(eos_coeffs(i)%s + up*(eos_coeffs(i)%s2 + up*eos_coeffs(i)%s3)) + dus = eos_coeffs(i)%s + up*(2._wp*eos_coeffs(i)%s2 + 3._wp*eos_coeffs(i)%s3*up) + up = up - (us*mu - up*(1._wp + mu))/(dus*mu - (1._wp + mu)) + end do + us = eos_coeffs(i)%c0 + up*(eos_coeffs(i)%s + up*(eos_coeffs(i)%s2 + up*eos_coeffs(i)%s3)) + dus = eos_coeffs(i)%s + up*(2._wp*eos_coeffs(i)%s2 + 3._wp*eos_coeffs(i)%s3*up) + dup_dmu = (up - us)/(dus*mu - (1._wp + mu)) + p_ref = eos_coeffs(i)%rho0*us*up + dp_dmu = eos_coeffs(i)%rho0*(dus*up + us)*dup_dmu + end if + e_ref = p_ref*mu/(2._wp*eos_coeffs(i)%rho0*(1._wp + mu)) + de_dmu = (dp_dmu*mu*(1._wp + mu) + p_ref)/(2._wp*eos_coeffs(i)%rho0*(1._wp + mu)**2) + dp_drho = dp_dmu/eos_coeffs(i)%rho0 + de_drho = de_dmu/eos_coeffs(i)%rho0 + case (eos_jwl) + ! JWL: p_ref = A exp(-R1 V) + B exp(-R2 V), V = rho0/rho. The curve is itself an isentrope, so de_ref = -p_ref d(1/rho). + V = eos_coeffs(i)%rho0/rho + ea = eos_coeffs(i)%a*exp(-eos_coeffs(i)%r1*V) + eb = eos_coeffs(i)%b*exp(-eos_coeffs(i)%r2*V) + p_ref = ea + eb + e_ref = (ea/eos_coeffs(i)%r1 + eb/eos_coeffs(i)%r2)/eos_coeffs(i)%rho0 + dp_drho = (eos_coeffs(i)%rho0/rho**2)*(eos_coeffs(i)%r1*ea + eos_coeffs(i)%r2*eb) + de_drho = p_ref/rho**2 + case (eos_vinet) + ! Vinet cold curve: p_c = 3 K0 (1 - x)/x^2 exp(eta (1 - x)), x = (rho0/rho)^(1/3), eta = 3 (K0' - 1)/2, + ! an isentrope like JWL (its energy integrates in closed form). + d = 1.5_wp*(eos_coeffs(i)%k0p - 1._wp) + x = (eos_coeffs(i)%rho0/rho)**(1._wp/3._wp) + ex = exp(d*(1._wp - x)) + p_ref = 3._wp*eos_coeffs(i)%k0*(1._wp - x)/x**2*ex + e_ref = 9._wp*eos_coeffs(i)%k0/(eos_coeffs(i)%rho0*d**2)*(1._wp - (1._wp - d*(1._wp - x))*ex) + dp_drho = 3._wp*eos_coeffs(i)%k0*ex*(-1._wp/x**2 - 2._wp*(1._wp - x)/x**3 - d*(1._wp - x)/x**2)*(-x/(3._wp*rho)) + de_drho = p_ref/rho**2 + end select + G0 = eos_coeffs(i)%gruneisen0 + eos_coeffs(i)%gruneisen_a*mu + dG0 = eos_coeffs(i)%gruneisen_a/eos_coeffs(i)%rho0 + + end subroutine s_reference_curve + + !> Whether the EOS of fluid i is a family whose coefficients vary with density. + function f_is_state_dependent(i) result(yes) + + $:GPU_ROUTINE(function_name='f_is_state_dependent', parallelism='[seq]', cray_inline=True) + + integer, intent(in) :: i + logical :: yes + + yes = ${EOS_IS_STATE_DEPENDENT('i')}$ + + end function f_is_state_dependent + + !> True when fluid i's reference curve is itself an isentrope (de_ref = -p_ref d(1/rho), which holds for JWL and Vinet but not + !! for the Mie-Gruneisen Hugoniot) and its Gruneisen coefficient is constant. Those two together make the isentrope through any + !! state closed-form, so it never has to be integrated. + function f_has_isentropic_reference(i) result(yes) + + $:GPU_ROUTINE(function_name='f_has_isentropic_reference', parallelism='[seq]', cray_inline=True) + + integer, intent(in) :: i + logical :: yes + + yes = (${EOS_HAS_ISENTROPIC_REFERENCE('i')}$) .and. eos_coeffs(i)%gruneisen_a == 0._wp + + end function f_has_isentropic_reference + + !> The largest compression a cubic Hugoniot fit can represent. mu(u_p) = u_p/(u_s - u_p) rises, peaks where c0 = s2 u_p^2 + 2 s3 + !! u_p^3, and falls after; only the rising branch is a physical shock. Returns a huge value for the linear fit, which never + !! turns over. Host-side: called once per fluid at initialization. + impure function f_hugoniot_compression_limit(c0, s, s2, s3) result(mu_max) + + real(wp), intent(in) :: c0, s, s2, s3 + real(wp) :: mu_max, up, f, df, us + integer :: iter + + if (s2 == 0._wp .and. s3 == 0._wp) then + mu_max = huge(1._wp) + return + end if + + ! Newton on c0 - s2 u^2 - 2 s3 u^3 = 0, from a guess that brackets the physical range + up = c0 + do iter = 1, 100 + f = c0 - s2*up**2 - 2._wp*s3*up**3 + df = -2._wp*s2*up - 6._wp*s3*up**2 + if (abs(df) < verysmall) exit + up = max(up - f/df, verysmall) + end do + us = c0 + up*(s + up*(s2 + up*s3)) + mu_max = up/max(us - up, verysmall) + + end function f_hugoniot_compression_limit + + !> Gamma, Pi, dPi/drho and dGamma/drho of fluid i at density rho, the coefficients of rho e = Gamma p + Pi(rho). Stiffened and + !! ideal gas keep the constants resolved at init, bit for bit. + subroutine s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: rho + integer, intent(in) :: i + real(wp), intent(out) :: gamma, pi_inf, dpi, dgamma + real(wp) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0 + + if (.not. f_is_state_dependent(i)) then + gamma = gammas(i) + pi_inf = pi_infs(i) + dpi = 0._wp + dgamma = 0._wp + return + end if + call s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, G0, dG0) + gamma = 1._wp/G0 + pi_inf = rho*e_ref - p_ref/G0 + dpi = e_ref + rho*de_drho - dp_drho/G0 + p_ref*dG0/G0**2 + dgamma = -dG0/G0**2 + + end subroutine s_eos_coefficients + + !> Exponent of the stiffened-gas isentrope p + B = const rho**n. Precomputed per fluid as isentrope_n. + function f_isentrope_exponent(gamma) result(n) + + $:GPU_ROUTINE(function_name='f_isentrope_exponent', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: gamma + real(wp) :: n + + n = 1._wp/gamma + 1._wp + + end function f_isentrope_exponent + + !> Reference pressure of that isentrope. Precomputed per fluid as isentrope_B. + function f_isentrope_pressure(pi_inf, gamma) result(B) + + $:GPU_ROUTINE(function_name='f_isentrope_pressure', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: pi_inf, gamma + real(wp) :: B + + B = pi_inf/(1._wp + gamma) + + end function f_isentrope_pressure + + !> Stiffened-gas thermal law p + B = (n - 1)*cv*rho*T. Pass rho to get T, or T to get rho. + function f_sg_thermal(pres, rho_or_T, n, B, cv) result(T_or_rho) + + $:GPU_ROUTINE(function_name='f_sg_thermal', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: pres, rho_or_T, n, B, cv + real(wp) :: T_or_rho + + T_or_rho = (pres + B)/((n - 1._wp)*cv*rho_or_T) + + end function f_sg_thermal + + !> Coefficients of phase i at its own density alpha_rho/alpha: the per-cell dispatch when some fluid's EOS is state dependent, + !! the constants resolved at init otherwise (bit for bit). + subroutine s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma) + + $:GPU_ROUTINE(function_name='s_phase_coefficients', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: alpha_rho, alpha + integer, intent(in) :: i + real(wp), intent(out) :: rho, gamma, pi_inf, dpi, dgamma + + rho = max(alpha_rho, sgm_eps)/max(alpha, sgm_eps) + if (any_state_dependent_eos) then + call s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma) + else + gamma = gammas(i) + pi_inf = pi_infs(i) + dpi = 0._wp + dgamma = 0._wp + end if + + end subroutine s_phase_coefficients + + !> c^2 = [((Gamma + 1) p + Pi)/rho - dPi/drho - p dGamma/drho]/Gamma, the frozen speed of one phase. + function f_c2_from_coefficients(rho, pres, gamma, pi_inf, dpi, dgamma) result(c2) + + $:GPU_ROUTINE(function_name='f_c2_from_coefficients', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: rho, pres, gamma, pi_inf, dpi, dgamma + real(wp) :: c2 + + c2 = (((gamma + 1._wp)*pres + pi_inf)/rho - dpi - pres*dgamma)/gamma + + end function f_c2_from_coefficients + + !> Frozen sound speed squared of one phase at (rho, p) from its own coefficients. These helpers are subroutines, not functions: + !! a device function that calls a device subroutine is a pattern no other backend-tested code in MFC uses. + subroutine s_phase_c2(rho, pres, i, c2) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: rho, pres + integer, intent(in) :: i + real(wp), intent(out) :: c2 + real(wp) :: gamma, pi_inf, dpi, dgamma + + call s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma) + c2 = f_c2_from_coefficients(rho, pres, gamma, pi_inf, dpi, dgamma) + + end subroutine s_phase_c2 + + !> Slope of the ODE `kind` for fluid i: dp/drho = c^2 along an isentrope (x = rho, y = p), or the reference temperature dT/dV = + !! (de_ref/dV + p_ref)/c_v - Gamma_G T/V (x = V, y = T), the Maxwell relation applied to e = e_ref + c_v (T - T_ref). + subroutine s_ode_slope(kind, i, x, y, dydx) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: kind, i + real(wp), intent(in) :: x, y + real(wp), intent(out) :: dydx + real(wp) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0 + + if (kind == ode_isentrope) then + call s_phase_c2(x, y, i, dydx) + else + call s_reference_curve(1._wp/x, i, p_ref, e_ref, dp_drho, de_drho, G0, dG0) + dydx = (p_ref - de_drho/x**2)/cvs(i) - G0*y/x + end if + + end subroutine s_ode_slope + + !> Fixed-step classical RK4 for the ODE `kind` from (x0, y0) to x1. + subroutine s_rk4(kind, i, x0, y0, x1, y) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: kind, i + real(wp), intent(in) :: x0, y0, x1 + real(wp), intent(out) :: y + real(wp) :: x, h, k1, k2, k3, k4 + integer :: step + + x = x0 + y = y0 + h = (x1 - x0)/eos_rk4_steps + $:GPU_LOOP(parallelism='[seq]') + do step = 1, eos_rk4_steps + call s_ode_slope(kind, i, x, y, k1) + call s_ode_slope(kind, i, x + 0.5_wp*h, y + 0.5_wp*h*k1, k2) + call s_ode_slope(kind, i, x + 0.5_wp*h, y + 0.5_wp*h*k2, k3) + call s_ode_slope(kind, i, x + h, y + h*k3, k4) + y = y + h*(k1 + 2._wp*(k2 + k3) + k4)/6._wp + x = x + h + end do + + end subroutine s_rk4 + + !> Pressure of phase i after the isentropic density change rho -> xi rho: closed form for the constant-coefficient families, + !! integrated for a state-dependent EOS (the star states it serves are close to rho). + subroutine s_phase_pressure_on_isentrope(pres, rho, xi, i, p_isen) + + $:GPU_ROUTINE(function_name='s_phase_pressure_on_isentrope', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: pres, rho, xi + integer, intent(in) :: i + real(wp), intent(out) :: p_isen + real(wp) :: p_ref_from, p_ref_to, e_ref, dp_drho, de_drho, G0, dG0 + + if (.not. f_is_state_dependent(i)) then + p_isen = (pres + isentrope_B(i))*xi**isentrope_n(i) - isentrope_B(i) + else if (f_has_isentropic_reference(i)) then + ! Exact: the offset from an isentropic reference obeys dDelta/Delta = Gamma drho/rho, so + ! p - p_ref scales as (rho'/rho)**(1 + Gamma). Integrating it instead costs a decimal per + ! doubling of the expansion and turns the pressure negative past roughly twentyfold. + call s_reference_curve(rho, i, p_ref_from, e_ref, dp_drho, de_drho, G0, dG0) + call s_reference_curve(xi*rho, i, p_ref_to, e_ref, dp_drho, de_drho, G0, dG0) + p_isen = p_ref_to + (pres - p_ref_from)*xi**(1._wp + eos_coeffs(i)%gruneisen0) + else + call s_rk4(ode_isentrope, i, rho, pres, xi*rho, p_isen) + end if + + end subroutine s_phase_pressure_on_isentrope + + !> Temperature of phase i at (rho, p): the stiffened-gas relation, or T_ref(rho) + (e - e_ref)/c_v. + subroutine s_phase_temperature(rho, pres, i, T) + + $:GPU_ROUTINE(function_name='s_phase_temperature', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: rho, pres + integer, intent(in) :: i + real(wp), intent(out) :: T + real(wp) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0, T0, T_ref + + if (f_is_state_dependent(i)) then + call s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, G0, dG0) + T0 = eos_coeffs(i)%t0 + call s_rk4(ode_reference_temperature, i, 1._wp/eos_coeffs(i)%rho0, T0, 1._wp/rho, T_ref) + T = T_ref + (pres - p_ref)/(rho*G0*cvs(i)) + else + T = (pres + isentrope_B(i))/((isentrope_n(i) - 1._wp)*cvs(i)*rho) + end if + + end subroutine s_phase_temperature + + !> Density of phase i on the isentrope through (rho_from, p_from) at p_to, and c^2 there: Newton on the pressure integrator, + !! whose slope is c^2. The relaxation's own Newton wraps this, so a few steps suffice. + subroutine s_phase_density_on_isentrope(i, rho_from, p_from, p_to, rho_to, c2_to) + + $:GPU_ROUTINE(function_name='s_phase_density_on_isentrope', parallelism='[seq]') + + integer, intent(in) :: i + real(wp), intent(in) :: rho_from, p_from, p_to + real(wp), intent(out) :: rho_to, c2_to + real(wp) :: p_at, c2_at + integer :: iter + + rho_to = rho_from + $:GPU_LOOP(parallelism='[seq]') + do iter = 1, 4 + call s_phase_pressure_on_isentrope(p_from, rho_from, rho_to/rho_from, i, p_at) + call s_phase_c2(rho_to, p_at, i, c2_at) + rho_to = rho_to - (p_at - p_to)/c2_at + end do + call s_phase_c2(rho_to, p_to, i, c2_to) + + end subroutine s_phase_density_on_isentrope + + !> Internal energy per unit volume of phase i at pressure pres: alpha (Gamma p + Pi) + alpha_rho qv, with the coefficients at + !! the phase's own density. + subroutine s_phase_internal_energy(pres, alpha, alpha_rho, i, e_phase) + + $:GPU_ROUTINE(function_name='s_phase_internal_energy', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: pres, alpha, alpha_rho + integer, intent(in) :: i + real(wp), intent(out) :: e_phase + real(wp) :: rho, gamma, pi_inf, dpi, dgamma + + call s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma) + e_phase = alpha*(gamma*pres + pi_inf) + alpha_rho*qvs(i) + + end subroutine s_phase_internal_energy + + !> Bulk modulus rho c^2 of phase i at pressure pres: f_bulk_modulus for a constant-coefficient fluid, bit for bit, minus the + !! reference-curve terms rho (dPi/drho + p dGamma/drho)/Gamma otherwise. + subroutine s_phase_bulk_modulus(pres, alpha, alpha_rho, i, blkmod) + + $:GPU_ROUTINE(function_name='s_phase_bulk_modulus', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: alpha_rho, alpha, pres + integer, intent(in) :: i + real(wp), intent(out) :: blkmod + real(wp) :: rho, gamma, pi_inf, dpi, dgamma + + call s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma) + blkmod = f_bulk_modulus(pres, gamma, pi_inf) - rho*(dpi + pres*dgamma)/gamma + + end subroutine s_phase_bulk_modulus + + !> Pressure of a stiffened gas from its internal energy density - the inverse of s_compute_energy. Callers subtract the kinetic, + !! magnetic and elastic energy first; none of those are equation-of-state terms. + function f_pressure(e_int, gamma, pi_inf, qv) result(pres) + + $:GPU_ROUTINE(function_name='f_pressure', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: e_int, gamma, pi_inf, qv + real(wp) :: pres + + pres = (e_int - pi_inf - qv)/gamma + + end function f_pressure + + !> Isentropic bulk modulus. Takes coefficients rather than a fluid index, so a mixture - whose effective gamma and pi_inf come + !! from s_compute_mixture_coefficients - is the same call as a single fluid. Elastic callers add their own shear term. + function f_bulk_modulus(pres, gamma, pi_inf) result(blkmod) + + $:GPU_ROUTINE(function_name='f_bulk_modulus', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: pres, gamma, pi_inf + real(wp) :: blkmod + + blkmod = ((gamma + 1._wp)*pres + pi_inf)/gamma + + end function f_bulk_modulus + + !> Relativistic specific enthalpy, h = 1 + (Gamma + 1)p/rho. Ideal gas only: the stiffness does not appear, so a fluid with a + !! nonzero pi_inf is not represented here (the validator refuses that combination). + function f_relativistic_enthalpy(pres, rho, gamma) result(H) + + $:GPU_ROUTINE(function_name='f_relativistic_enthalpy', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: pres, rho, gamma + real(wp) :: H + + H = 1._wp + (gamma + 1._wp)*pres/rho + + end function f_relativistic_enthalpy + +end module m_eos diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index a76a35c90..d8a083a75 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/src/common/m_global_parameters_common.fpp @@ -48,8 +48,7 @@ module m_global_parameters_common !> @} !> @name Material properties derived from fluid_pp - !> @{ One declaration is shared by all executables and initialized by m_variables_conversion after the case parameters have been - !! read. + !> @{ One declaration is shared by all executables and initialized by m_eos after the case parameters have been read. !> gammas is the stored form 1/(gamma - 1), not the ratio of specific heats; isentrope_n and isentrope_B are the same EOS !! written as p + B = const*rho**n. real(wp), allocatable, dimension(:) :: gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps diff --git a/src/common/m_phase_change.fpp b/src/common/m_phase_change.fpp index 0da551cb2..b9bd88aa7 100644 --- a/src/common/m_phase_change.fpp +++ b/src/common/m_phase_change.fpp @@ -12,6 +12,7 @@ module m_phase_change use m_global_parameters use m_mpi_proxy use m_variables_conversion + use m_eos use ieee_arithmetic use m_helper_basic use m_constants, only: model_eqns_6eq diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index a787a9727..80b19ed70 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -9,6 +9,7 @@ module m_variables_conversion use m_derived_types + use m_eos use m_global_parameters use m_mpi_proxy use m_helper_basic @@ -25,12 +26,9 @@ module m_variables_conversion & s_convert_mixture_to_mixture_variables, s_convert_species_to_mixture_variables, & & s_convert_species_to_mixture_variables_kernel, s_convert_conservative_to_primitive_variables, & & s_convert_primitive_to_conservative_variables, s_convert_primitive_to_flux_variables, s_compute_pressure, & - & s_compute_species_fraction, s_compute_mixture_coefficients, s_compute_energy, s_compute_speed_of_sound, f_bulk_modulus, & - & f_pressure, s_phase_internal_energy, f_isentrope_exponent, f_isentrope_pressure, f_sg_thermal, & + & s_compute_species_fraction, s_compute_mixture_coefficients, s_compute_energy, s_compute_speed_of_sound, & & s_compute_mixture_coefficients_dt, s_compute_speed_of_sound_avg, s_compute_fast_magnetosonic_speed, f_elastic_energy, & - & f_hypoelastic_energy, f_relativistic_enthalpy, s_eos_coefficients, s_phase_coefficients, s_phase_pressure_on_isentrope, & - & s_phase_temperature, f_is_state_dependent, s_phase_bulk_modulus, s_phase_density_on_isentrope, & - & s_finalize_variables_conversion_module, gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps + & f_hypoelastic_energy, s_finalize_variables_conversion_module, gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps real(wp), allocatable, dimension(:) :: Gs_vc integer, allocatable, dimension(:) :: bubrs_vc @@ -254,7 +252,6 @@ contains logical, optional, intent(in) :: enforce_density_floor, preserve_qbmm_number integer, optional, intent(in) :: lagrange_beta_index logical :: allocate_mixture_fields - logical :: state_dependent !< Whether this case's fluids need a density-dependent EOS allocate_mixture_fields = .false. if (present(store_mixture_fields)) allocate_mixture_fields = store_mixture_fields @@ -268,86 +265,12 @@ contains $:GPU_ENTER_DATA(copyin='[is1b, is1e, is2b, is2e, is3b, is3e]') $:GPU_UPDATE(device='[enforce_density_floor_vc, preserve_qbmm_number_vc, lagrange_beta_index_vc]') - @:ALLOCATE(gammas (1:num_fluids)) - @:ALLOCATE(eoss (1:num_fluids)) - @:ALLOCATE(isentrope_n (1:num_fluids)) - @:ALLOCATE(pi_infs(1:num_fluids)) - @:ALLOCATE(isentrope_B(1:num_fluids)) - @:ALLOCATE(cvs (1:num_fluids)) - @:ALLOCATE(qvs (1:num_fluids)) - @:ALLOCATE(qvps (1:num_fluids)) @:ALLOCATE(Gs_vc (1:num_fluids)) - state_dependent = .false. do i = 1, num_fluids - gammas(i) = fluid_pp(i)%gamma - isentrope_n(i) = f_isentrope_exponent(gammas(i)) - - ! Each EOS supplies its own coefficients. Resolved once here, not per cell: a branch in the mixture loop costs - ! registers in the Riemann kernels. An EOS whose coefficients depend on state must move to per-cell evaluation. - select case (fluid_pp(i)%eos) - case (eos_ideal_gas) - pi_infs(i) = 0._wp - case default - pi_infs(i) = fluid_pp(i)%pi_inf - end select Gs_vc(i) = fluid_pp(i)%G - isentrope_B(i) = f_isentrope_pressure(pi_infs(i), gammas(i)) - cvs(i) = fluid_pp(i)%cv - qvs(i) = fluid_pp(i)%qv - qvps(i) = fluid_pp(i)%qvp - eoss(i) = fluid_pp(i)%eos - eos_coeffs(i)%c0 = fluid_pp(i)%mg_c0 - eos_coeffs(i)%s = fluid_pp(i)%mg_s - eos_coeffs(i)%s2 = fluid_pp(i)%mg_s2 - eos_coeffs(i)%s3 = fluid_pp(i)%mg_s3 - ! Where a cubic Hugoniot fit turns over: mu(u_p) peaks where c0 = s2 u_p^2 + 2 s3 u_p^3, and past it - ! no shock state exists, so the Newton below would wander. Solved once here, on the host. - eos_coeffs(i)%mu_max = f_hugoniot_compression_limit(fluid_pp(i)%mg_c0, fluid_pp(i)%mg_s, fluid_pp(i)%mg_s2, & - & fluid_pp(i)%mg_s3) - eos_coeffs(i)%a = fluid_pp(i)%jwl_a - eos_coeffs(i)%b = fluid_pp(i)%jwl_b - eos_coeffs(i)%r1 = fluid_pp(i)%jwl_r1 - eos_coeffs(i)%r2 = fluid_pp(i)%jwl_r2 - eos_coeffs(i)%k0 = fluid_pp(i)%vinet_k0 - eos_coeffs(i)%k0p = fluid_pp(i)%vinet_k0p - ! One reference state and Gruneisen closure for every family; the user-facing names keep their prefix. - select case (fluid_pp(i)%eos) - case (eos_mie_gruneisen) - eos_coeffs(i)%rho0 = fluid_pp(i)%mg_rho0 - eos_coeffs(i)%t0 = fluid_pp(i)%mg_t0 - eos_coeffs(i)%gruneisen0 = fluid_pp(i)%mg_gruneisen - eos_coeffs(i)%gruneisen_a = fluid_pp(i)%mg_gruneisen_a - case (eos_jwl) - eos_coeffs(i)%rho0 = fluid_pp(i)%jwl_rho0 - eos_coeffs(i)%t0 = fluid_pp(i)%jwl_t0 - eos_coeffs(i)%gruneisen0 = fluid_pp(i)%jwl_omega - eos_coeffs(i)%gruneisen_a = 0._wp - case default - eos_coeffs(i)%rho0 = dflt_real - eos_coeffs(i)%t0 = dflt_real - eos_coeffs(i)%gruneisen0 = dflt_real - eos_coeffs(i)%gruneisen_a = 0._wp - case (eos_vinet) - eos_coeffs(i)%rho0 = fluid_pp(i)%vinet_rho0 - eos_coeffs(i)%t0 = fluid_pp(i)%vinet_t0 - eos_coeffs(i)%gruneisen0 = fluid_pp(i)%vinet_gruneisen - eos_coeffs(i)%gruneisen_a = fluid_pp(i)%vinet_gruneisen_a - end select - if (f_is_state_dependent(i)) state_dependent = .true. end do - #:if MFC_CASE_OPTIMIZATION - ! Baked in at build time, so a case that changed its EOS family since the build would silently - ! run the wrong branch. The namelist still carries fluid_pp%eos, so check the two agree. - @:PROHIBIT(state_dependent .neqv. any_state_dependent_eos, & - & "This case's equations of state do not match the ones this case-optimized binary was built for. Rebuild.") - #:else - any_state_dependent_eos = state_dependent - #:endif - $:GPU_UPDATE(device='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, eos_coeffs]') - #:if not MFC_CASE_OPTIMIZATION - $:GPU_UPDATE(device='[any_state_dependent_eos]') - #:endif + $:GPU_UPDATE(device='[Gs_vc]') @:ALLOCATE(Res_vc(1:2, 1:max(1, Re_size_max))) Res_vc = dflt_real @@ -1246,7 +1169,7 @@ contains if (allocated(rho_sf)) deallocate (rho_sf, gamma_sf, pi_inf_sf) - @:DEALLOCATE(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss) + @:DEALLOCATE(Gs_vc) if (allocated(bubrs_vc)) then @:DEALLOCATE(bubrs_vc) end if @@ -1367,394 +1290,6 @@ contains end subroutine s_compute_energy - !> The reference curve of a state-dependent EOS at rho: p_ref, e_ref, their d/drho, and Gamma_G with its d/drho. A new family - !! adds one case here and nothing else. - subroutine s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, G0, dG0) - - $:GPU_ROUTINE(parallelism='[seq]') - - real(wp), intent(in) :: rho - integer, intent(in) :: i - real(wp), intent(out) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0 - real(wp) :: mu, d, V, ea, eb, up, us, dus, dup_dmu, x, ex, dp_dmu, de_dmu - integer :: iter - - mu = rho/eos_coeffs(i)%rho0 - 1._wp - ! Past the fit's turnover there is no shock state to find; clamp rather than let the Newton below wander - ! off and return a silently wrong pressure. mu_max is huge for the linear fit, so this is a no-op there. - ! Bounded here so the step finishes and the host-side check in s_write_run_time_information can report it; - ! past the turnover there is no shock state and the Newton below would wander. - if (eoss(i) == eos_mie_gruneisen .and. mu > eos_coeffs(i)%mu_max) mu = eos_coeffs(i)%mu_max - select case (eoss(i)) - case (eos_mie_gruneisen) - ! Hugoniot reference u_s = c0 + s u_p + s2 u_p^2 + s3 u_p^3, with p_H = rho0 u_s u_p and the Hugoniot - ! energy e_H = p_H mu/(2 rho0 (1 + mu)); linear on release. Pole at mu = 1/(s - 1) for the linear fit; - ! the validator refuses initial states outside the EOS. - if (mu < 0._wp) then - p_ref = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2*mu - dp_dmu = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2 - else if (eos_coeffs(i)%s2 == 0._wp .and. eos_coeffs(i)%s3 == 0._wp) then - d = 1._wp - (eos_coeffs(i)%s - 1._wp)*mu - p_ref = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2*mu*(1._wp + mu)/(d*d) - dp_dmu = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2*((1._wp + 2._wp*mu)*d + 2._wp*(eos_coeffs(i)%s - 1._wp)*mu*(1._wp & - & + mu))/(d*d*d) - else - ! u_p solves u_s(u_p) mu = u_p (1 + mu): Newton from the linear fit, then implicit differentiation - up = eos_coeffs(i)%c0*mu/(1._wp - (eos_coeffs(i)%s - 1._wp)*mu) - $:GPU_LOOP(parallelism='[seq]') - do iter = 1, 8 - us = eos_coeffs(i)%c0 + up*(eos_coeffs(i)%s + up*(eos_coeffs(i)%s2 + up*eos_coeffs(i)%s3)) - dus = eos_coeffs(i)%s + up*(2._wp*eos_coeffs(i)%s2 + 3._wp*eos_coeffs(i)%s3*up) - up = up - (us*mu - up*(1._wp + mu))/(dus*mu - (1._wp + mu)) - end do - us = eos_coeffs(i)%c0 + up*(eos_coeffs(i)%s + up*(eos_coeffs(i)%s2 + up*eos_coeffs(i)%s3)) - dus = eos_coeffs(i)%s + up*(2._wp*eos_coeffs(i)%s2 + 3._wp*eos_coeffs(i)%s3*up) - dup_dmu = (up - us)/(dus*mu - (1._wp + mu)) - p_ref = eos_coeffs(i)%rho0*us*up - dp_dmu = eos_coeffs(i)%rho0*(dus*up + us)*dup_dmu - end if - e_ref = p_ref*mu/(2._wp*eos_coeffs(i)%rho0*(1._wp + mu)) - de_dmu = (dp_dmu*mu*(1._wp + mu) + p_ref)/(2._wp*eos_coeffs(i)%rho0*(1._wp + mu)**2) - dp_drho = dp_dmu/eos_coeffs(i)%rho0 - de_drho = de_dmu/eos_coeffs(i)%rho0 - case (eos_jwl) - ! JWL: p_ref = A exp(-R1 V) + B exp(-R2 V), V = rho0/rho. The curve is itself an isentrope, so de_ref = -p_ref d(1/rho). - V = eos_coeffs(i)%rho0/rho - ea = eos_coeffs(i)%a*exp(-eos_coeffs(i)%r1*V) - eb = eos_coeffs(i)%b*exp(-eos_coeffs(i)%r2*V) - p_ref = ea + eb - e_ref = (ea/eos_coeffs(i)%r1 + eb/eos_coeffs(i)%r2)/eos_coeffs(i)%rho0 - dp_drho = (eos_coeffs(i)%rho0/rho**2)*(eos_coeffs(i)%r1*ea + eos_coeffs(i)%r2*eb) - de_drho = p_ref/rho**2 - case (eos_vinet) - ! Vinet cold curve: p_c = 3 K0 (1 - x)/x^2 exp(eta (1 - x)), x = (rho0/rho)^(1/3), eta = 3 (K0' - 1)/2, - ! an isentrope like JWL (its energy integrates in closed form). - d = 1.5_wp*(eos_coeffs(i)%k0p - 1._wp) - x = (eos_coeffs(i)%rho0/rho)**(1._wp/3._wp) - ex = exp(d*(1._wp - x)) - p_ref = 3._wp*eos_coeffs(i)%k0*(1._wp - x)/x**2*ex - e_ref = 9._wp*eos_coeffs(i)%k0/(eos_coeffs(i)%rho0*d**2)*(1._wp - (1._wp - d*(1._wp - x))*ex) - dp_drho = 3._wp*eos_coeffs(i)%k0*ex*(-1._wp/x**2 - 2._wp*(1._wp - x)/x**3 - d*(1._wp - x)/x**2)*(-x/(3._wp*rho)) - de_drho = p_ref/rho**2 - end select - G0 = eos_coeffs(i)%gruneisen0 + eos_coeffs(i)%gruneisen_a*mu - dG0 = eos_coeffs(i)%gruneisen_a/eos_coeffs(i)%rho0 - - end subroutine s_reference_curve - - !> Whether the EOS of fluid i is a family whose coefficients vary with density. - function f_is_state_dependent(i) result(yes) - - $:GPU_ROUTINE(function_name='f_is_state_dependent', parallelism='[seq]', cray_inline=True) - - integer, intent(in) :: i - logical :: yes - - yes = eoss(i) == eos_mie_gruneisen .or. eoss(i) == eos_jwl .or. eoss(i) == eos_vinet - - end function f_is_state_dependent - - !> True when fluid i's reference curve is itself an isentrope (de_ref = -p_ref d(1/rho), which holds for JWL and Vinet but not - !! for the Mie-Gruneisen Hugoniot) and its Gruneisen coefficient is constant. Those two together make the isentrope through any - !! state closed-form, so it never has to be integrated. - function f_has_isentropic_reference(i) result(yes) - - $:GPU_ROUTINE(function_name='f_has_isentropic_reference', parallelism='[seq]', cray_inline=True) - - integer, intent(in) :: i - logical :: yes - - yes = (eoss(i) == eos_jwl .or. eoss(i) == eos_vinet) .and. eos_coeffs(i)%gruneisen_a == 0._wp - - end function f_has_isentropic_reference - - !> The largest compression a cubic Hugoniot fit can represent. mu(u_p) = u_p/(u_s - u_p) rises, peaks where c0 = s2 u_p^2 + 2 s3 - !! u_p^3, and falls after; only the rising branch is a physical shock. Returns a huge value for the linear fit, which never - !! turns over. Host-side: called once per fluid at initialization. - impure function f_hugoniot_compression_limit(c0, s, s2, s3) result(mu_max) - - real(wp), intent(in) :: c0, s, s2, s3 - real(wp) :: mu_max, up, f, df, us - integer :: iter - - if (s2 == 0._wp .and. s3 == 0._wp) then - mu_max = huge(1._wp) - return - end if - - ! Newton on c0 - s2 u^2 - 2 s3 u^3 = 0, from a guess that brackets the physical range - up = c0 - do iter = 1, 100 - f = c0 - s2*up**2 - 2._wp*s3*up**3 - df = -2._wp*s2*up - 6._wp*s3*up**2 - if (abs(df) < verysmall) exit - up = max(up - f/df, verysmall) - end do - us = c0 + up*(s + up*(s2 + up*s3)) - mu_max = up/max(us - up, verysmall) - - end function f_hugoniot_compression_limit - - !> Gamma, Pi, dPi/drho and dGamma/drho of fluid i at density rho, the coefficients of rho e = Gamma p + Pi(rho). Stiffened and - !! ideal gas keep the constants resolved at init, bit for bit. - subroutine s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma) - - $:GPU_ROUTINE(parallelism='[seq]') - - real(wp), intent(in) :: rho - integer, intent(in) :: i - real(wp), intent(out) :: gamma, pi_inf, dpi, dgamma - real(wp) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0 - - if (.not. f_is_state_dependent(i)) then - gamma = gammas(i) - pi_inf = pi_infs(i) - dpi = 0._wp - dgamma = 0._wp - return - end if - call s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, G0, dG0) - gamma = 1._wp/G0 - pi_inf = rho*e_ref - p_ref/G0 - dpi = e_ref + rho*de_drho - dp_drho/G0 + p_ref*dG0/G0**2 - dgamma = -dG0/G0**2 - - end subroutine s_eos_coefficients - - !> Exponent of the stiffened-gas isentrope p + B = const rho**n. Precomputed per fluid as isentrope_n. - function f_isentrope_exponent(gamma) result(n) - - $:GPU_ROUTINE(function_name='f_isentrope_exponent', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: gamma - real(wp) :: n - - n = 1._wp/gamma + 1._wp - - end function f_isentrope_exponent - - !> Reference pressure of that isentrope. Precomputed per fluid as isentrope_B. - function f_isentrope_pressure(pi_inf, gamma) result(B) - - $:GPU_ROUTINE(function_name='f_isentrope_pressure', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: pi_inf, gamma - real(wp) :: B - - B = pi_inf/(1._wp + gamma) - - end function f_isentrope_pressure - - !> Stiffened-gas thermal law p + B = (n - 1)*cv*rho*T. Pass rho to get T, or T to get rho. - function f_sg_thermal(pres, rho_or_T, n, B, cv) result(T_or_rho) - - $:GPU_ROUTINE(function_name='f_sg_thermal', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: pres, rho_or_T, n, B, cv - real(wp) :: T_or_rho - - T_or_rho = (pres + B)/((n - 1._wp)*cv*rho_or_T) - - end function f_sg_thermal - - !> Coefficients of phase i at its own density alpha_rho/alpha: the per-cell dispatch when some fluid's EOS is state dependent, - !! the constants resolved at init otherwise (bit for bit). - subroutine s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma) - - $:GPU_ROUTINE(function_name='s_phase_coefficients', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: alpha_rho, alpha - integer, intent(in) :: i - real(wp), intent(out) :: rho, gamma, pi_inf, dpi, dgamma - - rho = max(alpha_rho, sgm_eps)/max(alpha, sgm_eps) - if (any_state_dependent_eos) then - call s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma) - else - gamma = gammas(i) - pi_inf = pi_infs(i) - dpi = 0._wp - dgamma = 0._wp - end if - - end subroutine s_phase_coefficients - - !> c^2 = [((Gamma + 1) p + Pi)/rho - dPi/drho - p dGamma/drho]/Gamma, the frozen speed of one phase. - function f_c2_from_coefficients(rho, pres, gamma, pi_inf, dpi, dgamma) result(c2) - - $:GPU_ROUTINE(function_name='f_c2_from_coefficients', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: rho, pres, gamma, pi_inf, dpi, dgamma - real(wp) :: c2 - - c2 = (((gamma + 1._wp)*pres + pi_inf)/rho - dpi - pres*dgamma)/gamma - - end function f_c2_from_coefficients - - !> Frozen sound speed squared of one phase at (rho, p) from its own coefficients. These helpers are subroutines, not functions: - !! a device function that calls a device subroutine is a pattern no other backend-tested code in MFC uses. - subroutine s_phase_c2(rho, pres, i, c2) - - $:GPU_ROUTINE(parallelism='[seq]') - - real(wp), intent(in) :: rho, pres - integer, intent(in) :: i - real(wp), intent(out) :: c2 - real(wp) :: gamma, pi_inf, dpi, dgamma - - call s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma) - c2 = f_c2_from_coefficients(rho, pres, gamma, pi_inf, dpi, dgamma) - - end subroutine s_phase_c2 - - !> Slope of the ODE `kind` for fluid i: dp/drho = c^2 along an isentrope (x = rho, y = p), or the reference temperature dT/dV = - !! (de_ref/dV + p_ref)/c_v - Gamma_G T/V (x = V, y = T), the Maxwell relation applied to e = e_ref + c_v (T - T_ref). - subroutine s_ode_slope(kind, i, x, y, dydx) - - $:GPU_ROUTINE(parallelism='[seq]') - - integer, intent(in) :: kind, i - real(wp), intent(in) :: x, y - real(wp), intent(out) :: dydx - real(wp) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0 - - if (kind == ode_isentrope) then - call s_phase_c2(x, y, i, dydx) - else - call s_reference_curve(1._wp/x, i, p_ref, e_ref, dp_drho, de_drho, G0, dG0) - dydx = (p_ref - de_drho/x**2)/cvs(i) - G0*y/x - end if - - end subroutine s_ode_slope - - !> Fixed-step classical RK4 for the ODE `kind` from (x0, y0) to x1. - subroutine s_rk4(kind, i, x0, y0, x1, y) - - $:GPU_ROUTINE(parallelism='[seq]') - - integer, intent(in) :: kind, i - real(wp), intent(in) :: x0, y0, x1 - real(wp), intent(out) :: y - real(wp) :: x, h, k1, k2, k3, k4 - integer :: step - - x = x0 - y = y0 - h = (x1 - x0)/eos_rk4_steps - $:GPU_LOOP(parallelism='[seq]') - do step = 1, eos_rk4_steps - call s_ode_slope(kind, i, x, y, k1) - call s_ode_slope(kind, i, x + 0.5_wp*h, y + 0.5_wp*h*k1, k2) - call s_ode_slope(kind, i, x + 0.5_wp*h, y + 0.5_wp*h*k2, k3) - call s_ode_slope(kind, i, x + h, y + h*k3, k4) - y = y + h*(k1 + 2._wp*(k2 + k3) + k4)/6._wp - x = x + h - end do - - end subroutine s_rk4 - - !> Pressure of phase i after the isentropic density change rho -> xi rho: closed form for the constant-coefficient families, - !! integrated for a state-dependent EOS (the star states it serves are close to rho). - subroutine s_phase_pressure_on_isentrope(pres, rho, xi, i, p_isen) - - $:GPU_ROUTINE(parallelism='[seq]') - - real(wp), intent(in) :: pres, rho, xi - integer, intent(in) :: i - real(wp), intent(out) :: p_isen - real(wp) :: p_ref_from, p_ref_to, e_ref, dp_drho, de_drho, G0, dG0 - - if (.not. f_is_state_dependent(i)) then - p_isen = (pres + isentrope_B(i))*xi**isentrope_n(i) - isentrope_B(i) - else if (f_has_isentropic_reference(i)) then - ! Exact: the offset from an isentropic reference obeys dDelta/Delta = Gamma drho/rho, so - ! p - p_ref scales as (rho'/rho)**(1 + Gamma). Integrating it instead costs a decimal per - ! doubling of the expansion and turns the pressure negative past roughly twentyfold. - call s_reference_curve(rho, i, p_ref_from, e_ref, dp_drho, de_drho, G0, dG0) - call s_reference_curve(xi*rho, i, p_ref_to, e_ref, dp_drho, de_drho, G0, dG0) - p_isen = p_ref_to + (pres - p_ref_from)*xi**(1._wp + eos_coeffs(i)%gruneisen0) - else - call s_rk4(ode_isentrope, i, rho, pres, xi*rho, p_isen) - end if - - end subroutine s_phase_pressure_on_isentrope - - !> Temperature of phase i at (rho, p): the stiffened-gas relation, or T_ref(rho) + (e - e_ref)/c_v. - subroutine s_phase_temperature(rho, pres, i, T) - - $:GPU_ROUTINE(parallelism='[seq]') - - real(wp), intent(in) :: rho, pres - integer, intent(in) :: i - real(wp), intent(out) :: T - real(wp) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0, T0, T_ref - - if (f_is_state_dependent(i)) then - call s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, G0, dG0) - T0 = eos_coeffs(i)%t0 - call s_rk4(ode_reference_temperature, i, 1._wp/eos_coeffs(i)%rho0, T0, 1._wp/rho, T_ref) - T = T_ref + (pres - p_ref)/(rho*G0*cvs(i)) - else - T = (pres + isentrope_B(i))/((isentrope_n(i) - 1._wp)*cvs(i)*rho) - end if - - end subroutine s_phase_temperature - - !> Density of phase i on the isentrope through (rho_from, p_from) at p_to, and c^2 there: Newton on the pressure integrator, - !! whose slope is c^2. The relaxation's own Newton wraps this, so a few steps suffice. - subroutine s_phase_density_on_isentrope(i, rho_from, p_from, p_to, rho_to, c2_to) - - $:GPU_ROUTINE(function_name='s_phase_density_on_isentrope', parallelism='[seq]') - - integer, intent(in) :: i - real(wp), intent(in) :: rho_from, p_from, p_to - real(wp), intent(out) :: rho_to, c2_to - real(wp) :: p_at, c2_at - integer :: iter - - rho_to = rho_from - $:GPU_LOOP(parallelism='[seq]') - do iter = 1, 4 - call s_phase_pressure_on_isentrope(p_from, rho_from, rho_to/rho_from, i, p_at) - call s_phase_c2(rho_to, p_at, i, c2_at) - rho_to = rho_to - (p_at - p_to)/c2_at - end do - call s_phase_c2(rho_to, p_to, i, c2_to) - - end subroutine s_phase_density_on_isentrope - - !> Internal energy per unit volume of phase i at pressure pres: alpha (Gamma p + Pi) + alpha_rho qv, with the coefficients at - !! the phase's own density. - subroutine s_phase_internal_energy(pres, alpha, alpha_rho, i, e_phase) - - $:GPU_ROUTINE(function_name='s_phase_internal_energy', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: pres, alpha, alpha_rho - integer, intent(in) :: i - real(wp), intent(out) :: e_phase - real(wp) :: rho, gamma, pi_inf, dpi, dgamma - - call s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma) - e_phase = alpha*(gamma*pres + pi_inf) + alpha_rho*qvs(i) - - end subroutine s_phase_internal_energy - - !> Bulk modulus rho c^2 of phase i at pressure pres: f_bulk_modulus for a constant-coefficient fluid, bit for bit, minus the - !! reference-curve terms rho (dPi/drho + p dGamma/drho)/Gamma otherwise. - subroutine s_phase_bulk_modulus(pres, alpha, alpha_rho, i, blkmod) - - $:GPU_ROUTINE(function_name='s_phase_bulk_modulus', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: alpha_rho, alpha, pres - integer, intent(in) :: i - real(wp), intent(out) :: blkmod - real(wp) :: rho, gamma, pi_inf, dpi, dgamma - - call s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma) - blkmod = f_bulk_modulus(pres, gamma, pi_inf) - rho*(dpi + pres*dgamma)/gamma - - end subroutine s_phase_bulk_modulus - !> Elastic strain energy of one stress component, doubled for a shear component: the tensor stores it once, the energy counts !! both off-diagonal entries. Zero without a shear modulus. function f_elastic_energy(tau, G, is_shear) result(dE) @@ -1789,45 +1324,6 @@ contains end function f_hypoelastic_energy - !> Pressure of a stiffened gas from its internal energy density - the inverse of s_compute_energy. Callers subtract the kinetic, - !! magnetic and elastic energy first; none of those are equation-of-state terms. - function f_pressure(e_int, gamma, pi_inf, qv) result(pres) - - $:GPU_ROUTINE(function_name='f_pressure', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: e_int, gamma, pi_inf, qv - real(wp) :: pres - - pres = (e_int - pi_inf - qv)/gamma - - end function f_pressure - - !> Isentropic bulk modulus. Takes coefficients rather than a fluid index, so a mixture - whose effective gamma and pi_inf come - !! from s_compute_mixture_coefficients - is the same call as a single fluid. Elastic callers add their own shear term. - function f_bulk_modulus(pres, gamma, pi_inf) result(blkmod) - - $:GPU_ROUTINE(function_name='f_bulk_modulus', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: pres, gamma, pi_inf - real(wp) :: blkmod - - blkmod = ((gamma + 1._wp)*pres + pi_inf)/gamma - - end function f_bulk_modulus - - !> Relativistic specific enthalpy, h = 1 + (Gamma + 1)p/rho. Ideal gas only: the stiffness does not appear, so a fluid with a - !! nonzero pi_inf is not represented here (the validator refuses that combination). - function f_relativistic_enthalpy(pres, rho, gamma) result(H) - - $:GPU_ROUTINE(function_name='f_relativistic_enthalpy', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: pres, rho, gamma - real(wp) :: H - - H = 1._wp + (gamma + 1._wp)*pres/rho - - end function f_relativistic_enthalpy - !> Speed of sound of a thermodynamic state. Enthalpy is not an argument: for a real state H, |u|^2 and qv all cancel out of c^2 !! = ((Gamma + 1)p + Pi)/(Gamma rho). Averaged states, whose enthalpy is a free input, use the _avg variant. subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho) diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index b8606d4e5..795ae7e73 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -12,6 +12,7 @@ module m_data_output use m_compile_specific use m_helper use m_variables_conversion + use m_eos use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, format_silo, format_binary, precision_single implicit none diff --git a/src/post_process/m_derived_variables.fpp b/src/post_process/m_derived_variables.fpp index 49085c5b4..105512143 100644 --- a/src/post_process/m_derived_variables.fpp +++ b/src/post_process/m_derived_variables.fpp @@ -11,6 +11,7 @@ module m_derived_variables use m_mpi_proxy use m_helper_basic use m_variables_conversion + use m_eos use m_constants, only: model_eqns_gamma_law implicit none diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 86a00c001..8db04e1d7 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -17,6 +17,7 @@ module m_start_up use m_boundary_common use m_boundary_io use m_variables_conversion + use m_eos use m_data_input use m_data_output use m_derived_variables @@ -792,6 +793,7 @@ contains call s_initialize_mpi_common_module(exchange_all_chemistry_temperatures_in=.true., use_rdma_transport_in=.false.) end if call s_initialize_boundary_common_module() + call s_initialize_eos_module() call s_initialize_variables_conversion_module(store_mixture_fields=.true., lagrange_beta_index=beta_idx) call s_initialize_data_input_module() call s_initialize_derived_variables_module() @@ -996,6 +998,7 @@ contains call s_finalize_derived_variables_module() call s_finalize_data_input_module() call s_finalize_variables_conversion_module() + call s_finalize_eos_module() if (num_procs > 1) then call s_finalize_mpi_proxy_module() call s_finalize_mpi_common_module() diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 04e399737..476327b8d 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -11,6 +11,7 @@ module m_start_up use m_global_parameters use m_mpi_proxy use m_mpi_common + use m_eos, only: s_initialize_eos_module, s_finalize_eos_module use m_variables_conversion use m_grid use m_initial_condition @@ -478,6 +479,7 @@ contains end if call s_initialize_mpi_common_module(exchange_all_chemistry_temperatures_in=.false., use_rdma_transport_in=.false.) call s_initialize_data_output_module() + call s_initialize_eos_module() call s_initialize_variables_conversion_module() call s_initialize_grid_module() call s_initialize_initial_condition_module() @@ -652,6 +654,7 @@ contains call s_finalize_mpi_common_module() call s_finalize_grid_module() call s_finalize_variables_conversion_module() + call s_finalize_eos_module() call s_finalize_data_output_module() call s_finalize_global_parameters_module() call s_finalize_assign_variables_module() diff --git a/src/simulation/m_acoustic_src.fpp b/src/simulation/m_acoustic_src.fpp index 24f9ccec3..5ca556ede 100644 --- a/src/simulation/m_acoustic_src.fpp +++ b/src/simulation/m_acoustic_src.fpp @@ -11,6 +11,7 @@ module m_acoustic_src use m_global_parameters use m_bubbles use m_variables_conversion + use m_eos use m_helper_basic use m_constants diff --git a/src/simulation/m_bubbles_EE.fpp b/src/simulation/m_bubbles_EE.fpp index 258b94db8..c1f954742 100644 --- a/src/simulation/m_bubbles_EE.fpp +++ b/src/simulation/m_bubbles_EE.fpp @@ -11,6 +11,7 @@ module m_bubbles_EE use m_global_parameters use m_mpi_proxy use m_variables_conversion + use m_eos use m_bubbles implicit none diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index 4ff7eb531..15a709c17 100644 --- a/src/simulation/m_bubbles_EL.fpp +++ b/src/simulation/m_bubbles_EL.fpp @@ -12,6 +12,7 @@ module m_bubbles_EL use m_bubbles_EL_kernels use m_bubbles use m_variables_conversion + use m_eos use m_compile_specific use m_boundary_common use m_helper_basic diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp index f6d6f0080..187fa70ce 100644 --- a/src/simulation/m_hypoelastic.fpp +++ b/src/simulation/m_hypoelastic.fpp @@ -11,7 +11,7 @@ module m_hypoelastic use m_global_parameters use m_finite_differences use m_helper - use m_variables_conversion, only: f_bulk_modulus, s_phase_bulk_modulus + use m_eos, only: s_phase_bulk_modulus implicit none diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 99a212260..8f7d24b0e 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -12,6 +12,7 @@ module m_ibm use m_global_parameters use m_mpi_proxy use m_variables_conversion + use m_eos use m_helper use m_helper_basic use m_constants diff --git a/src/simulation/m_igr.fpp b/src/simulation/m_igr.fpp index e8d06502f..846a039a9 100644 --- a/src/simulation/m_igr.fpp +++ b/src/simulation/m_igr.fpp @@ -11,6 +11,7 @@ module m_igr use m_derived_types use m_global_parameters use m_variables_conversion + use m_eos use m_mpi_proxy use m_helper use m_boundary_common diff --git a/src/simulation/m_pressure_relaxation.fpp b/src/simulation/m_pressure_relaxation.fpp index 71fb0d079..6c6dfe6dd 100644 --- a/src/simulation/m_pressure_relaxation.fpp +++ b/src/simulation/m_pressure_relaxation.fpp @@ -12,8 +12,8 @@ module m_pressure_relaxation use m_derived_types use m_global_parameters use m_mpi_proxy, only: s_mpi_abort - use m_variables_conversion, only: s_convert_species_to_mixture_variables_kernel, f_pressure, s_phase_internal_energy, & - & s_phase_coefficients, s_phase_density_on_isentrope, f_is_state_dependent + use m_variables_conversion, only: s_convert_species_to_mixture_variables_kernel + use m_eos, only: f_pressure, s_phase_internal_energy, s_phase_coefficients, s_phase_density_on_isentrope, f_is_state_dependent implicit none diff --git a/src/simulation/m_qbmm.fpp b/src/simulation/m_qbmm.fpp index fbd50093d..cde6830a7 100644 --- a/src/simulation/m_qbmm.fpp +++ b/src/simulation/m_qbmm.fpp @@ -12,6 +12,7 @@ module m_qbmm use m_global_parameters use m_mpi_proxy use m_variables_conversion + use m_eos use m_helper_basic use m_helper use m_constants, only: bubble_model_keller_miksis, bubble_model_rayleigh_plesset diff --git a/src/simulation/m_reactive_burn.fpp b/src/simulation/m_reactive_burn.fpp index ee312c8e6..73da758cb 100644 --- a/src/simulation/m_reactive_burn.fpp +++ b/src/simulation/m_reactive_burn.fpp @@ -15,7 +15,8 @@ module m_reactive_burn use m_global_parameters - use m_variables_conversion, only: s_phase_temperature, s_compute_mixture_coefficients, f_pressure + use m_variables_conversion, only: s_compute_mixture_coefficients + use m_eos, only: s_phase_temperature, f_pressure implicit none diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index edd500ef1..efe9c40b9 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -13,6 +13,7 @@ module m_rhs use m_global_parameters use m_mpi_proxy use m_variables_conversion + use m_eos use m_weno use m_constants, only: riemann_solver_hll, riemann_solver_hlld, model_eqns_6eq, int_comp_mthinc, recon_type_weno, & & recon_type_muscl diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index a5ffec343..94084f645 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -11,6 +11,7 @@ module m_riemann_solver_hll use m_derived_types use m_global_parameters use m_variables_conversion + use m_eos use m_constants, only: riemann_solver_hll, riemann_solver_hllc, riemann_solver_lax_friedrichs, avg_state_roe, & & avg_state_arithmetic, wave_speeds_direct, wave_speeds_pressure use m_chemistry diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index 621f00fff..7006785cb 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -11,6 +11,7 @@ module m_riemann_solver_hllc use m_derived_types use m_global_parameters use m_variables_conversion + use m_eos use m_bubbles use m_constants, only: riemann_solver_hll, riemann_solver_hllc, riemann_solver_lax_friedrichs, model_eqns_5eq, & & model_eqns_6eq, avg_state_roe, avg_state_arithmetic, wave_speeds_direct, wave_speeds_pressure diff --git a/src/simulation/m_riemann_solver_hypo_hlld.fpp b/src/simulation/m_riemann_solver_hypo_hlld.fpp index b4d955b07..2fb798e6a 100644 --- a/src/simulation/m_riemann_solver_hypo_hlld.fpp +++ b/src/simulation/m_riemann_solver_hypo_hlld.fpp @@ -11,6 +11,7 @@ module m_riemann_solver_hypo_hlld use m_derived_types use m_global_parameters use m_variables_conversion + use m_eos use m_riemann_state implicit none diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 6a5434f5b..3a2a7e139 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -13,6 +13,7 @@ module m_start_up use m_mpi_proxy use m_mpi_common use m_variables_conversion + use m_eos use m_weno use m_muscl use m_thinc @@ -836,6 +837,7 @@ contains end if call s_initialize_mpi_common_module(exchange_all_chemistry_temperatures_in=.false., use_rdma_transport_in=rdma_mpi) call s_initialize_mpi_proxy_module() + call s_initialize_eos_module() call s_initialize_variables_conversion_module(enforce_density_floor=.true., preserve_qbmm_number=.true.) if (grid_geometry == 3) call s_initialize_fftw_module() @@ -1141,6 +1143,7 @@ contains end if if (int_comp > 0) call s_finalize_thinc_module() call s_finalize_variables_conversion_module() + call s_finalize_eos_module() if (grid_geometry == 3) call s_finalize_fftw_module call s_finalize_mpi_common_module() call s_finalize_global_parameters_module() diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py index 6dd02de5b..0359e4523 100644 --- a/toolchain/mfc/case.py +++ b/toolchain/mfc/case.py @@ -9,6 +9,7 @@ from . import common from .analytic_expr import AnalyticExprError, fortranize_expr +from .params.eos_families import EOS_FAMILIES from .printer import cons from .run import case_dicts from .state import ARG @@ -32,6 +33,10 @@ def _suggest_similar_params(unknown_key: str, valid_keys: list, n: int = 3) -> l "Bz": "eqn_idx%B%end", } +# fluid_pp(i)%eos values whose coefficients depend on the local state; derived so a new +# family never leaves any_state_dependent_eos stale (see EOS_FAMILIES). +EOS_STATE_DEPENDENT_VALUES = frozenset(f.value for f in EOS_FAMILIES if f.state_dependent) + MIBM_ANALYTIC_VARS = ["vel(1)", "vel(2)", "vel(3)", "angular_vel(1)", "angular_vel(2)", "angular_vel(3)"] # "eqn_idx%B%end - 1" not "eqn_idx%B%beg + 1" must be used because 1D does not have Bx @@ -373,7 +378,7 @@ def __get_sim_fpp(self, print: bool) -> str: # Baking this lets the compiler drop the state-dependent EOS chain entirely. Left in the # call graph it costs registers, and so occupancy, in every kernel that can reach it. - eos_state_dependent = {3, 4, 5} # Mie-Gruneisen, JWL, Vinet; see eos_* in m_constants.fpp + eos_state_dependent = EOS_STATE_DEPENDENT_VALUES num_fluids_case = int(self.params.get("num_fluids", 1)) any_state_dependent_eos = 1 if any(int(self.params.get(f"fluid_pp({f})%eos", 1)) in eos_state_dependent for f in range(1, num_fluids_case + 1)) else 0 diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 1da790a90..ea0f9c273 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -20,6 +20,7 @@ from . import eos from .common import MFCException from .params.definitions import CONSTRAINTS +from .params.eos_families import EOS_FAMILIES from .params.namelist_parser import get_fortran_constants from .state import CFG @@ -960,20 +961,27 @@ def check_ibm(self): elif model_id is not None and model_id > 0: self.prohibit(True, f"patch_icpp({i})%model_id is set but geometry ({geometry}) is not an STL model (21)") - def _check_initial_states_inside_eos(self, num_fluids, eos_names): + def _eos_coefficient_args(self, i, family): + """Fluid i's arguments to family.coefficients_fn, after rho, in registry order. An optional + parameter defaults to 0.0, matching the Fortran `case default`; a required one is passed + as-is so a missing value raises TypeError and the parameter rules report it.""" + optional = {suffix for suffix, _math in family.optional} + args = [] + for suffix in family.coefficients_args: + value = self.get(f"fluid_pp({i})%{family.prefix}_{suffix}") + args.append((value or 0.0) if suffix in optional else value) + return args + + def _check_initial_states_inside_eos(self, num_fluids): """Every patch must start each state-dependent fluid where rho e > 0 and c^2 > 0; the solver has no clamp.""" num_patches = self.get("num_patches", 0) or 0 + families = {f.value: f for f in EOS_FAMILIES if f.state_dependent} for i in range(1, num_fluids + 1): - eos_ = self.get(f"fluid_pp({i})%eos") - g = lambda k: self.get(f"fluid_pp({i})%{k}") # noqa: E731 - if eos_ == eos_names["mie_gruneisen"]: - coefficients = lambda r: eos.eos_coefficients(r, g("mg_rho0"), g("mg_c0"), g("mg_s"), g("mg_gruneisen"), g("mg_gruneisen_a") or 0.0, g("mg_s2") or 0.0, g("mg_s3") or 0.0) # noqa: E731 - elif eos_ == eos_names["jwl"]: - coefficients = lambda r: eos.jwl_coefficients(r, g("jwl_rho0"), g("jwl_a"), g("jwl_b"), g("jwl_r1"), g("jwl_r2"), g("jwl_omega")) # noqa: E731 - elif eos_ == eos_names.get("vinet"): - coefficients = lambda r: eos.coefficients_from_curve(r, eos.vinet_reference(r, g("vinet_rho0"), g("vinet_k0"), g("vinet_k0p")), g("vinet_gruneisen")) # noqa: E731 - else: + family = families.get(self.get(f"fluid_pp({i})%eos")) + if family is None: continue + fn = getattr(eos, family.coefficients_fn) + coefficients = lambda r, f=family, fn=fn: fn(r, *self._eos_coefficient_args(i, f)) # noqa: E731 for j in range(1, num_patches + 1): ar, a, p = (self.get(f"patch_icpp({j})%alpha_rho({i})"), self.get(f"patch_icpp({j})%alpha({i})"), self.get(f"patch_icpp({j})%pres")) if not all(isinstance(x, (int, float)) for x in (ar, a, p)) or a <= 0: @@ -999,13 +1007,10 @@ def check_eos_selector(self): return eos_names = CONSTRAINTS["fluid_pp(1)%eos"]["names"] eos_ideal_gas = eos_names["ideal_gas"] - # The state-dependent families: selector value -> (parameter prefix, its parameters) - families = { - eos_names["mie_gruneisen"]: ("mg", ("rho0", "c0", "s", "gruneisen")), - eos_names["jwl"]: ("jwl", ("a", "b", "r1", "r2", "omega", "rho0")), - eos_names["vinet"]: ("vinet", ("k0", "k0p", "rho0", "gruneisen")), - } - optional = {"mg": ("gruneisen_a", "t0", "s2", "s3"), "jwl": ("t0",), "vinet": ("gruneisen_a", "t0")} + # The state-dependent families: selector value -> (parameter prefix, its required parameters) + families = {f.value: (f.prefix, tuple(n for n, _ in f.required)) for f in EOS_FAMILIES if f.state_dependent} + optional = {f.prefix: tuple(n for n, _ in f.optional) for f in EOS_FAMILIES if f.state_dependent} + state_dependent_names = ", ".join(f.suffix for f in EOS_FAMILIES if f.state_dependent) bub_fac = 1 if self.get("bubbles_euler", "F") == "T" else 0 state_dependent = {} for i in range(1, num_fluids + 1 + bub_fac): @@ -1061,16 +1066,16 @@ def check_eos_selector(self): if self.get("T_wrt", "F") == "T" or (i == 1 and self._is_numeric(rta) and rta > 0): t0 = self.get(f"fluid_pp({i})%{prefix}_t0") self.prohibit(t0 is None or t0 <= 0, f"the temperature of fluid {i} needs fluid_pp({i})%{prefix}_t0 > 0") - self._check_initial_states_inside_eos(num_fluids, eos_names) + self._check_initial_states_inside_eos(num_fluids) # The per-phase evaluation is wired through the 5-equation paths only; every feature below still # reads the stiffened-gas coefficients directly. - self.prohibit(self.get("model_eqns") not in (2, 3), "a state-dependent eos (mie_gruneisen, jwl, vinet) requires model_eqns = 2 or 3") - self.prohibit(self.get("riemann_solver") not in (1, 2, 5), "a state-dependent eos (mie_gruneisen, jwl, vinet) requires riemann_solver = 1, 2 or 5") - self.prohibit(self.get("wave_speeds") == 2, "a state-dependent eos (mie_gruneisen, jwl, vinet) requires wave_speeds = 1 (the PVRS estimate is stiffened-gas only)") + self.prohibit(self.get("model_eqns") not in (2, 3), f"a state-dependent eos ({state_dependent_names}) requires model_eqns = 2 or 3") + self.prohibit(self.get("riemann_solver") not in (1, 2, 5), f"a state-dependent eos ({state_dependent_names}) requires riemann_solver = 1, 2 or 5") + self.prohibit(self.get("wave_speeds") == 2, f"a state-dependent eos ({state_dependent_names}) requires wave_speeds = 1 (the PVRS estimate is stiffened-gas only)") for j in range(1, (self.get("num_patches") or 0) + 1): self.prohibit(self.get(f"patch_icpp({j})%hcid") in (202, 203), f"patch_icpp({j})%hcid = 202/203 read fluid_pp(1)%gamma, which a state-dependent eos does not set") for flag in ("bubbles_euler", "bubbles_lagrange", "igr", "relativity", "mhd", "chemistry", "relax"): - self.prohibit(self.get(flag, "F") == "T", f"a state-dependent eos (mie_gruneisen, jwl, vinet) is not supported with {flag} = T") + self.prohibit(self.get(flag, "F") == "T", f"a state-dependent eos ({state_dependent_names}) is not supported with {flag} = T") def check_stiffened_eos(self): """Checks constraints on stiffened equation of state fluids parameters""" @@ -2046,8 +2051,8 @@ def check_reactive_burn(self): # differing only in qv; violating these silently corrupts the mass/energy balance. self.prohibit(self.get("num_fluids") != 2, "reactive_burn requires num_fluids = 2 (reactant then product) to be set") # A state-dependent family carries its own curve; the shared-EOS check is a stiffened-gas one. - names = CONSTRAINTS["fluid_pp(1)%eos"]["names"] - state_dependent = any(self.get(f"fluid_pp({k})%eos") in (names["mie_gruneisen"], names["jwl"], names["vinet"]) for k in (1, 2)) + state_dependent_values = {f.value for f in EOS_FAMILIES if f.state_dependent} + state_dependent = any(self.get(f"fluid_pp({k})%eos") in state_dependent_values for k in (1, 2)) for prop in () if state_dependent else ("gamma", "pi_inf"): v1 = self.get(f"fluid_pp(1)%{prop}") v2 = self.get(f"fluid_pp(2)%{prop}") diff --git a/toolchain/mfc/eos.py b/toolchain/mfc/eos.py index eeb8db982..e344d9ba4 100644 --- a/toolchain/mfc/eos.py +++ b/toolchain/mfc/eos.py @@ -128,3 +128,9 @@ def vinet_reference(rho, rho0, k0, k0p): dx_drho = -x / (3.0 * rho) dp_dx = 3.0 * k0 * ex * ((-1.0) / x**2 - 2.0 * (1.0 - x) / x**3 - eta * (1.0 - x) / x**2) return p, e, dp_dx * dx_drho, p / rho**2 + + +def vinet_coefficients(rho, rho0, k0, k0p, gruneisen, gruneisen_a=0.0): + """Vinet cold curve with Gamma_G = Gamma_0 + a mu: what s_eos_coefficients returns.""" + mu = rho / rho0 - 1.0 + return coefficients_from_curve(rho, vinet_reference(rho, rho0, k0, k0p), gruneisen + gruneisen_a * mu, gruneisen_a / rho0) diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index ae1ecc028..ecbd679ea 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -7,6 +7,7 @@ import re from typing import Any, Dict +from .eos_families import EOS_FAMILIES from .namelist_parser import get_fortran_constants from .registry import REGISTRY, IndexedFamily from .schema import ParamDef, ParamType @@ -173,7 +174,7 @@ def _lookup_hint(name): # Schema Validation for Constraints and Dependencies # Uses rapidfuzz for "did you mean?" suggestions when typos are detected -_VALID_CONSTRAINT_KEYS = {"choices", "min", "max", "value_labels", "names"} +_VALID_CONSTRAINT_KEYS = {"choices", "min", "max", "value_labels", "names", "fortran_prefix"} _VALID_DEPENDENCY_KEYS = {"when_true", "when_set", "when_value"} _VALID_CONDITION_KEYS = {"requires", "recommends", "requires_value"} @@ -216,6 +217,11 @@ def _validate_constraint(param_name: str, constraint: Dict[str, Any]) -> None: raise ValueError(f"names for '{param_name}' map two names to the same value") if "choices" in constraint and set(names.values()) != set(constraint["choices"]): raise ValueError(f"names for '{param_name}' must cover exactly its choices {constraint['choices']}") + if "fortran_prefix" in constraint: + if "names" not in constraint: + raise ValueError(f"Constraint 'fortran_prefix' for '{param_name}' requires 'names'") + if not isinstance(constraint["fortran_prefix"], str) or not re.match(r"^[a-z0-9][a-z0-9_]*$", constraint["fortran_prefix"]): + raise ValueError(f"Constraint 'fortran_prefix' for '{param_name}' must be a lowercase identifier") def _validate_dependency(param_name: str, dependency: Dict[str, Any]) -> None: @@ -910,9 +916,10 @@ def _load(): for mm in range(-ll, ll + 1): _r(f"{px}sph_har_coeff({ll},{mm})", REAL) - # Values must match the hand-written eos_* constants in src/common/m_constants.fpp. - _EOS_NAMES = {"stiffened_gas": 1, "ideal_gas": 2, "mie_gruneisen": 3, "jwl": 4, "vinet": 5} - _EOS_VALUE_LABELS = {1: "stiffened-gas", 2: "ideal-gas", 3: "Mie-Gruneisen", 4: "JWL", 5: "Vinet"} + # Derived from EOS_FAMILIES (eos_families.py), which must match the eos_* constants in src/common/m_constants.fpp. + _EOS_NAMES = {f.suffix: f.value for f in EOS_FAMILIES} + _EOS_VALUE_LABELS = {f.value: f.label for f in EOS_FAMILIES} + _EOS_CHOICES = sorted(f.value for f in EOS_FAMILIES) # fluid_pp (10 fluids) # Members present in physical_parameters: gamma, pi_inf, Re, cv, qv, qvp, G. @@ -920,40 +927,15 @@ def _load(): # by upstream #1085/#1093 — they must NOT be registered (namelist read would crash). for f in range(1, NF + 1): px = f"fluid_pp({f})%" - CONSTRAINTS[f"fluid_pp({f})%eos"] = {"choices": [1, 2, 3, 4, 5], "value_labels": _EOS_VALUE_LABELS, "names": _EOS_NAMES} + CONSTRAINTS[f"fluid_pp({f})%eos"] = {"choices": _EOS_CHOICES, "value_labels": _EOS_VALUE_LABELS, "names": _EOS_NAMES, "fortran_prefix": "eos"} for a, sym in [("gamma", r"\f$\gamma_k\f$"), ("pi_inf", r"\f$\pi_{\infty,k}\f$"), ("cv", r"\f$c_{v,k}\f$"), ("qv", r"\f$q_{v,k}\f$"), ("qvp", r"\f$q'_{v,k}\f$")]: _r(f"{px}{a}", REAL, math=sym) _r(f"{px}eos", INT, math=r"\f$\mathrm{EOS}_k\f$") - for a, sym in [ - ("mg_rho0", r"\f$\rho_{0,k}\f$"), - ("mg_c0", r"\f$c_{0,k}\f$"), - ("mg_s", r"\f$s_k\f$"), - ("mg_gruneisen", r"\f$\Gamma_{G,k}\f$"), - ("mg_gruneisen_a", r"\f$a_k\f$"), - ("mg_t0", r"\f$T_{0,k}\f$"), - ("mg_s2", r"\f$s_{2,k}\f$"), - ("mg_s3", r"\f$s_{3,k}\f$"), - ]: - _r(f"{px}{a}", REAL, math=sym) - for a, sym in [ - ("jwl_a", r"\f$A_k\f$"), - ("jwl_b", r"\f$B_k\f$"), - ("jwl_r1", r"\f$R_{1,k}\f$"), - ("jwl_r2", r"\f$R_{2,k}\f$"), - ("jwl_omega", r"\f$\omega_k\f$"), - ("jwl_rho0", r"\f$\rho_{0,k}\f$"), - ("jwl_t0", r"\f$T_{0,k}\f$"), - ]: - _r(f"{px}{a}", REAL, math=sym) - for a, sym in [ - ("vinet_k0", r"\f$K_{0,k}\f$"), - ("vinet_k0p", r"\f$K'_{0,k}\f$"), - ("vinet_rho0", r"\f$\rho_{0,k}\f$"), - ("vinet_gruneisen", r"\f$\Gamma_{G,k}\f$"), - ("vinet_gruneisen_a", r"\f$a_k\f$"), - ("vinet_t0", r"\f$T_{0,k}\f$"), - ]: - _r(f"{px}{a}", REAL, math=sym) + for fam in EOS_FAMILIES: + if fam.prefix is None: + continue + for suffix, sym in fam.required + fam.optional: + _r(f"{px}{fam.prefix}_{suffix}", REAL, math=sym) _r(f"{px}G", REAL, {"hypoelasticity"}, math=r"\f$G_k\f$") _r(f"{px}Re(1)", REAL, {"viscosity"}, math=r"\f$\mathrm{Re}_k\f$ (shear)") _r(f"{px}Re(2)", REAL, {"viscosity"}, math=r"\f$\mathrm{Re}_k\f$ (bulk)") diff --git a/toolchain/mfc/params/eos_families.py b/toolchain/mfc/params/eos_families.py new file mode 100644 index 000000000..0e402fa08 --- /dev/null +++ b/toolchain/mfc/params/eos_families.py @@ -0,0 +1,158 @@ +"""The equation-of-state family registry: one entry per family, everything mechanical derives from it. + +See docs/superpowers/specs/2026-09-12-eos-family-registry-design.md. `required`/`optional` are +ordered (parameter suffix, Doxygen math symbol) pairs; the full parameter name is +`fluid_pp(i)%_`. `eos_coeffs` maps a state-dependent family's `eos_coeffs` fields +(src/common/m_derived_types.fpp) to where their value comes from: a parameter (`Param`), a fixed +Fortran literal (`FortranLiteral`), or a value a Fortran helper computes (`Computed`) — not every field is +a copy of a parameter, so a flat name->name dict cannot express all three. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Union + + +@dataclass(frozen=True) +class Param: + """eos_coeffs field copied from fluid_pp(i)%_.""" + + suffix: str + + +@dataclass(frozen=True) +class FortranLiteral: + """eos_coeffs field set to a fixed Fortran literal; no parameter backs it (e.g. JWL's gruneisen_a).""" + + fortran: str + + +@dataclass(frozen=True) +class Computed: + """eos_coeffs field derived by a Fortran helper, not copied or literal (e.g. mu_max).""" + + fortran_call: str + + +EosCoeffSource = Union[Param, FortranLiteral, Computed] + + +@dataclass(frozen=True) +class EosFamily: + value: int # fluid_pp(i)%eos selector + suffix: str # Fortran constant is eos_ + label: str # _EOS_VALUE_LABELS display text + prefix: str | None # parameter prefix (e.g. "mg"); None for stiffened_gas/ideal_gas + required: tuple[tuple[str, str], ...] = () + optional: tuple[tuple[str, str], ...] = () + state_dependent: bool = False + isentropic_reference: bool = False # reference curve is itself an isentrope (JWL, Vinet) + coefficients_fn: str | None = None # name of the mirror function in toolchain/mfc/eos.py + # Parameter suffixes to pass to coefficients_fn after rho, in call order. A suffix listed in + # `optional` is passed as `value or 0.0`; one listed in `required` is passed as-is. + coefficients_args: tuple[str, ...] = () + eos_coeffs: dict[str, EosCoeffSource] = field(default_factory=dict) + + +EOS_FAMILIES: tuple[EosFamily, ...] = ( + EosFamily(value=1, suffix="stiffened_gas", label="stiffened-gas", prefix=None), + EosFamily(value=2, suffix="ideal_gas", label="ideal-gas", prefix=None), + EosFamily( + value=3, + suffix="mie_gruneisen", + label="Mie-Gruneisen", + prefix="mg", + required=( + ("rho0", r"\f$\rho_{0,k}\f$"), + ("c0", r"\f$c_{0,k}\f$"), + ("s", r"\f$s_k\f$"), + ("gruneisen", r"\f$\Gamma_{G,k}\f$"), + ), + optional=( + ("gruneisen_a", r"\f$a_k\f$"), + ("t0", r"\f$T_{0,k}\f$"), + ("s2", r"\f$s_{2,k}\f$"), + ("s3", r"\f$s_{3,k}\f$"), + ), + state_dependent=True, + coefficients_fn="eos_coefficients", + coefficients_args=("rho0", "c0", "s", "gruneisen", "gruneisen_a", "s2", "s3"), + eos_coeffs={ + "c0": Param("c0"), + "s": Param("s"), + "s2": Param("s2"), + "s3": Param("s3"), + "rho0": Param("rho0"), + "t0": Param("t0"), + "gruneisen0": Param("gruneisen"), + "gruneisen_a": Param("gruneisen_a"), + "mu_max": Computed("f_hugoniot_compression_limit(mg_c0, mg_s, mg_s2, mg_s3)"), + }, + ), + EosFamily( + value=4, + suffix="jwl", + label="JWL", + prefix="jwl", + required=( + ("a", r"\f$A_k\f$"), + ("b", r"\f$B_k\f$"), + ("r1", r"\f$R_{1,k}\f$"), + ("r2", r"\f$R_{2,k}\f$"), + ("omega", r"\f$\omega_k\f$"), + ("rho0", r"\f$\rho_{0,k}\f$"), + ), + optional=(("t0", r"\f$T_{0,k}\f$"),), + state_dependent=True, + isentropic_reference=True, + coefficients_fn="jwl_coefficients", + coefficients_args=("rho0", "a", "b", "r1", "r2", "omega"), + eos_coeffs={ + "a": Param("a"), + "b": Param("b"), + "r1": Param("r1"), + "r2": Param("r2"), + "rho0": Param("rho0"), + "t0": Param("t0"), + "gruneisen0": Param("omega"), + "gruneisen_a": FortranLiteral("0._wp"), # JWL has no gruneisen_a parameter + }, + ), + EosFamily( + value=5, + suffix="vinet", + label="Vinet", + prefix="vinet", + required=( + ("k0", r"\f$K_{0,k}\f$"), + ("k0p", r"\f$K'_{0,k}\f$"), + ("rho0", r"\f$\rho_{0,k}\f$"), + ("gruneisen", r"\f$\Gamma_{G,k}\f$"), + ), + optional=( + ("gruneisen_a", r"\f$a_k\f$"), + ("t0", r"\f$T_{0,k}\f$"), + ), + state_dependent=True, + isentropic_reference=True, + coefficients_fn="vinet_coefficients", + coefficients_args=("rho0", "k0", "k0p", "gruneisen", "gruneisen_a"), + eos_coeffs={ + "k0": Param("k0"), + "k0p": Param("k0p"), + "rho0": Param("rho0"), + "t0": Param("t0"), + "gruneisen0": Param("gruneisen"), + "gruneisen_a": Param("gruneisen_a"), + }, + ), +) + + +# The `case default` arm of s_initialize_eos_module: what a fluid whose family is not +# state-dependent gets for each family-dispatched field. Not a family property, so it cannot +# live on an EosFamily. dflt_real marks "unused"; gruneisen_a must be a live zero because +# s_reference_curve reads it unconditionally. +EOS_COEFF_DEFAULT = "dflt_real" +EOS_COEFF_DEFAULTS: dict[str, str] = {"gruneisen_a": "0._wp"} diff --git a/toolchain/mfc/params/generators/fortran_gen.py b/toolchain/mfc/params/generators/fortran_gen.py index 26fcf1510..605c97a9f 100644 --- a/toolchain/mfc/params/generators/fortran_gen.py +++ b/toolchain/mfc/params/generators/fortran_gen.py @@ -1,5 +1,6 @@ """Fortran parameter code generator — namelist and scalar decl fragments per target.""" +import re from pathlib import Path from typing import List, Tuple @@ -259,19 +260,176 @@ def generate_constants_fpp() -> str: from ..definitions import CONSTRAINTS lines = [_HEADER.rstrip()] + emitted_prefixes = set() # e.g. fluid_pp(1..10)%eos all share fortran_prefix "eos"; emit once for param in sorted(CONSTRAINTS): - names = CONSTRAINTS[param].get("names") + constraint = CONSTRAINTS[param] + names = constraint.get("names") if not names: continue - # Compound keys (fluid_pp(1)%eos) do not form valid Fortran identifiers, so their - # constants are hand-written in m_constants.fpp and guarded by test_eos_selector.py. + # Compound keys (fluid_pp(1)%eos) do not form valid Fortran identifiers on their own; + # a "fortran_prefix" supplies the standalone name (eos_) to emit instead. + prefix = constraint.get("fortran_prefix") if "%" in param or "(" in param: + if prefix is None: + continue + else: + prefix = param + if prefix in emitted_prefixes: continue + emitted_prefixes.add(prefix) for name, value in sorted(names.items(), key=lambda kv: kv[1]): - lines.append(f"integer, parameter :: {param}_{name} = {value}") + lines.append(f"integer, parameter :: {prefix}_{name} = {value}") + return "\n".join(lines) + "\n" + + +_EOS_INDENT = " " # macro-body indent, matching src/common/include/macros.fpp +_EOS_CONT = " & " # continuation prefix for a wrapped assignment +_EOS_WRAP = 120 # wrap a rendered assignment past this; the hard Fortran free-form limit is 132 +_EOS_CALL_RE = re.compile(r"^(\w+)\((.*)\)$") + + +def _eos_case_fields() -> set: + """eos_coeffs fields that s_initialize_eos_module must dispatch on the family. + + A field more than one family writes needs a `case` arm to pick the right source; a field a + single family writes has one source and can be assigned for every fluid, which is what the + hand-written init did (mg_c0 and friends are copied regardless of family). + """ + from ..eos_families import EOS_FAMILIES + + writers: dict = {} + for family in EOS_FAMILIES: + for field_name in family.eos_coeffs: + writers[field_name] = writers.get(field_name, 0) + 1 + return {name for name, count in writers.items() if count > 1} + + +def _eos_computed_call(family, call: str, idx: str) -> str: + """Resolve a Computed field's call, qualifying each argument as fluid_pp(i)%.""" + match = _EOS_CALL_RE.match(call.strip()) + if match is None: + raise ValueError(f"EOS family {family.suffix}: Computed fortran_call {call!r} is not a plain () call, so its arguments cannot be qualified.") + known = {f"{family.prefix}_{suffix}" for suffix, _ in family.required + family.optional} + args = [] + for raw in match.group(2).split(","): + arg = raw.strip() + if arg not in known: + raise ValueError(f"EOS family {family.suffix}: Computed fortran_call argument {arg!r} is not one of the family's parameters ({', '.join(sorted(known))}).") + args.append(f"fluid_pp({idx})%{arg}") + return f"{match.group(1)}({', '.join(args)})" + + +def _eos_rhs(family, source, idx: str) -> str: + """The Fortran right-hand side an eos_coeffs source resolves to.""" + from ..eos_families import Computed, FortranLiteral, Param + + if isinstance(source, Param): + return f"fluid_pp({idx})%{family.prefix}_{source.suffix}" + if isinstance(source, FortranLiteral): + return source.fortran + if isinstance(source, Computed): + return _eos_computed_call(family, source.fortran_call, idx) + raise ValueError(f"EOS family {family.suffix}: unknown eos_coeffs source {source!r}.") + + +def _eos_assign(indent: str, field_name: str, rhs: str, idx: str) -> List[str]: + """`eos_coeffs(i)% = `, split at argument boundaries if it renders too long.""" + lhs = f"{indent}eos_coeffs({idx})%{field_name} = " + rendered = len(lhs.replace(f"({idx})", "(i)")) + len(rhs.replace(f"({idx})", "(i)")) + if rendered <= _EOS_WRAP: + return [lhs + rhs] + parts = rhs.split(", ") + lines, current = [], lhs + parts[0] + for part in parts[1:]: + candidate = f"{current}, {part}" + if len(candidate.replace(f"({idx})", "(i)")) > _EOS_WRAP: + lines.append(current + ", &") + current = _EOS_CONT + part + else: + current = candidate + lines.append(current) + return lines + + +def _eos_family_test(families, idx: str) -> str: + """`eoss(i) == eos_a .or. eoss(i) == eos_b`, the membership test for a set of families.""" + return " .or. ".join(f"eoss({idx}) == eos_{family.suffix}" for family in families) + + +def generate_eos_fpp() -> str: + """Fypp macros for the mechanical parts of m_eos: the family predicates and the coefficient init. + + Identical for all targets. The predicates are expression-valued, so `${MACRO('i')}$` splices + them into a larger condition; the init macros are statement-valued, called with `@:MACRO(i)`. + """ + from ..eos_families import EOS_COEFF_DEFAULT, EOS_COEFF_DEFAULTS, EOS_FAMILIES + + idx = "${i}$" + case_fields = _eos_case_fields() + # A family that skips a case-dispatched field would emit an arm leaving it uninitialised + # (case default is not taken when another arm matches) -- worse than a dflt_real sentinel. + for family in EOS_FAMILIES: + if family.state_dependent: + missing = case_fields - set(family.eos_coeffs) + if missing: + raise ValueError(f"EOS family {family.suffix!r} does not assign case-dispatched field(s) {sorted(missing)}") + unknown_defaults = set(EOS_COEFF_DEFAULTS) - case_fields + if unknown_defaults: + raise ValueError(f"EOS_COEFF_DEFAULTS key(s) {sorted(unknown_defaults)} are not case-dispatched fields") + lines = [_HEADER.rstrip()] + lines.append("#! Generated from EOS_FAMILIES in toolchain/mfc/params/eos_families.py.") + lines.append("") + lines.append("#! The families whose coefficients vary with density.") + lines.append("#:def EOS_IS_STATE_DEPENDENT(i)") + lines.append(_eos_family_test([f for f in EOS_FAMILIES if f.state_dependent], idx)) + lines.append("#:enddef") + lines.append("") + lines.append("#! The families whose reference curve is itself an isentrope. Only the family half of the") + lines.append("#! predicate: the runtime gruneisen_a test is not a family property and stays in m_eos.fpp.") + lines.append("#:def EOS_HAS_ISENTROPIC_REFERENCE(i)") + lines.append(_eos_family_test([f for f in EOS_FAMILIES if f.isentropic_reference], idx)) + lines.append("#:enddef") + lines.append("") + lines.append("#! Fields exactly one family writes: one source each, so every fluid gets them.") + lines.append("#:def EOS_INIT_COEFFS(i)") + for family in EOS_FAMILIES: + for field_name, source in family.eos_coeffs.items(): + if field_name not in case_fields: + lines.extend(_eos_assign(_EOS_INDENT, field_name, _eos_rhs(family, source, idx), idx)) + lines.append("#:enddef") + lines.append("") + lines.append("#! Fields several families write: dispatched on the family. case default is the") + lines.append("#! non-state-dependent families' arm.") + lines.append("#:def EOS_INIT_REFERENCE_STATE(i)") + lines.append(f"{_EOS_INDENT}select case (fluid_pp({idx})%eos)") + body = _EOS_INDENT * 2 + for family in EOS_FAMILIES: + written = [(n, s) for n, s in family.eos_coeffs.items() if n in case_fields] + if not written: + continue + lines.append(f"{_EOS_INDENT}case (eos_{family.suffix})") + for field_name, source in written: + lines.extend(_eos_assign(body, field_name, _eos_rhs(family, source, idx), idx)) + lines.append(f"{_EOS_INDENT}case default") + for field_name in _eos_case_field_order(case_fields): + lines.extend(_eos_assign(body, field_name, EOS_COEFF_DEFAULTS.get(field_name, EOS_COEFF_DEFAULT), idx)) + lines.append(f"{_EOS_INDENT}end select") + lines.append("#:enddef") return "\n".join(lines) + "\n" +def _eos_case_field_order(case_fields: set) -> List[str]: + """The case-dispatched fields in registry order, so case default matches the arms above it.""" + from ..eos_families import EOS_FAMILIES + + order = [] + for family in EOS_FAMILIES: + for field_name in family.eos_coeffs: + if field_name in case_fields and field_name not in order: + order.append(field_name) + return order + + # case.py-computed extras that are not CASE_OPT_PARAMS but appear in the # case-optimization declaration block (injected as Fypp #:set variables # by toolchain/mfc/case.py, not from the namelist registry). @@ -694,10 +852,10 @@ def resolve_namelist_content(fpp_path: Path) -> str: def get_generated_files(build_dir: Path) -> List[Tuple[Path, str]]: - """Return (path, content) for all 15 generated .fpp files under build_dir. + """Return (path, content) for all 18 generated .fpp files under build_dir. Paths match the cmake include directory structure: - build_dir/include/{full_target}/generated_{namelist,decls,constants,case_opt_decls,bcast}.fpp + build_dir/include/{full_target}/generated_{namelist,decls,constants,eos,case_opt_decls,bcast}.fpp Every target gets generated_case_opt_decls.fpp: the full case-optimization block for simulation and common computed-scalar declarations for pre/post. Every target gets generated_bcast.fpp with its MPI broadcast statements. @@ -708,6 +866,7 @@ def get_generated_files(build_dir: Path) -> List[Tuple[Path, str]]: result.append((inc / "generated_namelist.fpp", generate_namelist_fpp(short))) result.append((inc / "generated_decls.fpp", generate_decls_fpp(short))) result.append((inc / "generated_constants.fpp", generate_constants_fpp())) + result.append((inc / "generated_eos.fpp", generate_eos_fpp())) sim_gpu_decls = "" for short, full in TARGETS: inc = build_dir / "include" / full diff --git a/toolchain/mfc/params_tests/test_eos_families.py b/toolchain/mfc/params_tests/test_eos_families.py new file mode 100644 index 000000000..b7c33603b --- /dev/null +++ b/toolchain/mfc/params_tests/test_eos_families.py @@ -0,0 +1,240 @@ +"""R1: proves the EOS_FAMILIES registry changed nothing observable. + +R5 (below the R1 tests): proves the registry agrees with the parts of the codebase that are +still hand-written by design — physical_parameters, s_initialize_eos_module, eos.py and +s_reference_curve — so a family that validates cleanly but computes with dflt_real coefficients +is caught here rather than at runtime. + +Expected values are hard-coded (captured from the tree before the registry existed), not derived +from the registry itself — asserting the registry against itself would be circular. +""" + +import re + +from .. import eos as eos_module +from ..case import EOS_STATE_DEPENDENT_VALUES +from ..params import REGISTRY +from ..params.definitions import CONSTRAINTS +from ..params.eos_families import EOS_COEFF_DEFAULTS, EOS_FAMILIES +from ..params.generators.fortran_gen import generate_eos_fpp +from ..params.namelist_parser import get_mfc_root + +_EXPECTED_NAMES = {"stiffened_gas": 1, "ideal_gas": 2, "mie_gruneisen": 3, "jwl": 4, "vinet": 5} +_EXPECTED_LABELS = {1: "stiffened-gas", 2: "ideal-gas", 3: "Mie-Gruneisen", 4: "JWL", 5: "Vinet"} +_EXPECTED_CHOICES = [1, 2, 3, 4, 5] + +# Captured from `sorted(n for n in REGISTRY.all_params if n.startswith("fluid_pp(1)%"))` +# on the unmodified tree (HEAD 949a131c), before eos_families.py existed. +_EXPECTED_FLUID_PP1_PARAMS = [ + "fluid_pp(1)%G", + "fluid_pp(1)%K", + "fluid_pp(1)%Re(1)", + "fluid_pp(1)%Re(2)", + "fluid_pp(1)%cv", + "fluid_pp(1)%eos", + "fluid_pp(1)%gamma", + "fluid_pp(1)%hb_m", + "fluid_pp(1)%jwl_a", + "fluid_pp(1)%jwl_b", + "fluid_pp(1)%jwl_omega", + "fluid_pp(1)%jwl_r1", + "fluid_pp(1)%jwl_r2", + "fluid_pp(1)%jwl_rho0", + "fluid_pp(1)%jwl_t0", + "fluid_pp(1)%mg_c0", + "fluid_pp(1)%mg_gruneisen", + "fluid_pp(1)%mg_gruneisen_a", + "fluid_pp(1)%mg_rho0", + "fluid_pp(1)%mg_s", + "fluid_pp(1)%mg_s2", + "fluid_pp(1)%mg_s3", + "fluid_pp(1)%mg_t0", + "fluid_pp(1)%mu_bulk", + "fluid_pp(1)%mu_max", + "fluid_pp(1)%mu_min", + "fluid_pp(1)%nn", + "fluid_pp(1)%non_newtonian", + "fluid_pp(1)%pi_inf", + "fluid_pp(1)%qv", + "fluid_pp(1)%qvp", + "fluid_pp(1)%tau0", + "fluid_pp(1)%vinet_gruneisen", + "fluid_pp(1)%vinet_gruneisen_a", + "fluid_pp(1)%vinet_k0", + "fluid_pp(1)%vinet_k0p", + "fluid_pp(1)%vinet_rho0", + "fluid_pp(1)%vinet_t0", +] + + +def test_eos_names_unchanged(): + assert CONSTRAINTS["fluid_pp(1)%eos"]["names"] == _EXPECTED_NAMES + + +def test_eos_value_labels_unchanged(): + assert CONSTRAINTS["fluid_pp(1)%eos"]["value_labels"] == _EXPECTED_LABELS + + +def test_eos_choices_unchanged(): + assert CONSTRAINTS["fluid_pp(1)%eos"]["choices"] == _EXPECTED_CHOICES + + +def test_fluid_pp_parameter_set_unchanged(): + names = sorted(n for n in REGISTRY.all_params if n.startswith("fluid_pp(1)%")) + assert names == _EXPECTED_FLUID_PP1_PARAMS + + +# R5: agreement tests. Each closes one leg of the web docs/superpowers/specs/2026-09-12-eos- +# family-registry-design.md maps between the registry and the hand-written code it drives. + +_M_EOS = (get_mfc_root() / "src" / "common" / "m_eos.fpp").read_text() +_M_DERIVED_TYPES = (get_mfc_root() / "src" / "common" / "m_derived_types.fpp").read_text() + +_STATE_DEPENDENT = [f for f in EOS_FAMILIES if f.state_dependent] + + +def _physical_parameters_fields(): + """Field names declared on `type physical_parameters`, splitting comma-joined `::` lines + (e.g. `real(wp) :: mg_s2, mg_s3` is two fields, not one).""" + body = re.search(r"type physical_parameters\b(.*?)end type physical_parameters", _M_DERIVED_TYPES, re.S).group(1) + fields = set() + for line in body.splitlines(): + code = line.split("!", 1)[0] + if "::" not in code: + continue + for raw in code.split("::", 1)[1].split(","): + name = raw.strip() + if name: + fields.add(name) + return fields + + +def test_every_registry_parameter_has_a_physical_parameters_field(): + fields = _physical_parameters_fields() + for family in EOS_FAMILIES: + for suffix, _math in family.required + family.optional: + name = f"{family.prefix}_{suffix}" + assert name in fields, f"{name} (family {family.suffix}) has no physical_parameters field" + + +def _init_module_text(): + """s_initialize_eos_module's body, plus the generated macro bodies it calls (`@:EOS_INIT_COEFFS` + / `@:EOS_INIT_REFERENCE_STATE`) — the parameters are read via those macros, not literally in + m_eos.fpp, since R4 moved them into generated_eos.fpp.""" + body = re.search(r"impure subroutine s_initialize_eos_module\(\)(.*?)end subroutine s_initialize_eos_module", _M_EOS, re.S).group(1) + return body + "\n" + generate_eos_fpp() + + +def test_every_parameter_is_read_by_the_init(): + """Optionals included: the bug this registry exists to prevent (Vinet ignoring vinet_gruneisen_a) + was an optional parameter, so a required-only sweep would still miss it.""" + text = _init_module_text() + for family in _STATE_DEPENDENT: + for suffix, _math in family.required + family.optional: + name = f"{family.prefix}_{suffix}" + assert re.search(rf"\b{re.escape(name)}\b", text), f"{name} (family {family.suffix}) is never read in s_initialize_eos_module" + + +def test_every_coefficients_fn_resolves(): + for family in _STATE_DEPENDENT: + fn = getattr(eos_module, family.coefficients_fn, None) + assert callable(fn), f"eos.py has no callable {family.coefficients_fn!r} (family {family.suffix})" + + +def _reference_curve_text(): + return re.search(r"subroutine s_reference_curve\(.*?\)(.*?)end subroutine s_reference_curve", _M_EOS, re.S).group(1) + + +def test_every_family_has_a_reference_curve_case(): + text = _reference_curve_text() + for family in _STATE_DEPENDENT: + assert f"case (eos_{family.suffix})" in text, f"no case (eos_{family.suffix}) in s_reference_curve" + + +# Verified by hand against the registry (docs/superpowers/plans/2026-09-12-eos-family-registry.md, +# R1 review). Hard-coded so an edit that garbles a LaTeX symbol while keeping the parameter name +# is caught here rather than slipping through unreviewed. +_EXPECTED_MATH = { + ("mie_gruneisen", "rho0"): r"\f$\rho_{0,k}\f$", + ("mie_gruneisen", "c0"): r"\f$c_{0,k}\f$", + ("mie_gruneisen", "s"): r"\f$s_k\f$", + ("mie_gruneisen", "gruneisen"): r"\f$\Gamma_{G,k}\f$", + ("mie_gruneisen", "gruneisen_a"): r"\f$a_k\f$", + ("mie_gruneisen", "t0"): r"\f$T_{0,k}\f$", + ("mie_gruneisen", "s2"): r"\f$s_{2,k}\f$", + ("mie_gruneisen", "s3"): r"\f$s_{3,k}\f$", + ("jwl", "a"): r"\f$A_k\f$", + ("jwl", "b"): r"\f$B_k\f$", + ("jwl", "r1"): r"\f$R_{1,k}\f$", + ("jwl", "r2"): r"\f$R_{2,k}\f$", + ("jwl", "omega"): r"\f$\omega_k\f$", + ("jwl", "rho0"): r"\f$\rho_{0,k}\f$", + ("jwl", "t0"): r"\f$T_{0,k}\f$", + ("vinet", "k0"): r"\f$K_{0,k}\f$", + ("vinet", "k0p"): r"\f$K'_{0,k}\f$", + ("vinet", "rho0"): r"\f$\rho_{0,k}\f$", + ("vinet", "gruneisen"): r"\f$\Gamma_{G,k}\f$", + ("vinet", "gruneisen_a"): r"\f$a_k\f$", + ("vinet", "t0"): r"\f$T_{0,k}\f$", +} + + +def test_doxygen_math_symbols_match_expected(): + actual = {(family.suffix, suffix): math for family in EOS_FAMILIES for suffix, math in family.required + family.optional} + assert actual == _EXPECTED_MATH + + +def _case_dispatched_fields(): + """eos_coeffs fields written by more than one family — the same rule generate_eos_fpp's + _eos_case_fields uses, reimplemented here so this test does not depend on that function.""" + writers = {} + for family in EOS_FAMILIES: + for field_name in family.eos_coeffs: + writers[field_name] = writers.get(field_name, 0) + 1 + return {name for name, count in writers.items() if count > 1} + + +def test_every_state_dependent_family_assigns_every_case_dispatched_field(): + case_fields = _case_dispatched_fields() + assert case_fields # sanity: the set this test exists to check is non-empty + for family in _STATE_DEPENDENT: + missing = case_fields - set(family.eos_coeffs) + assert not missing, f"family {family.suffix} does not assign case-dispatched field(s) {sorted(missing)}" + + +def test_eos_coeff_defaults_keys_are_case_dispatched_fields(): + case_fields = _case_dispatched_fields() + unknown = set(EOS_COEFF_DEFAULTS) - case_fields + assert not unknown, f"EOS_COEFF_DEFAULTS key(s) {sorted(unknown)} are not case-dispatched fields" + + +# The set case.py bakes into any_state_dependent_eos. Hard-coded, not derived: deriving it from +# EOS_FAMILIES would assert the registry against itself and prove nothing. +_EXPECTED_STATE_DEPENDENT_VALUES = {3, 4, 5} + + +def test_case_optimization_state_dependent_values_unchanged(): + assert set(EOS_STATE_DEPENDENT_VALUES) == _EXPECTED_STATE_DEPENDENT_VALUES + + +# Captured from the hand-written if/elif/else the validator used before the dispatch was generated. +# A reordered argument silently changes which initial states are accepted, so pin the order. +_EXPECTED_COEFFICIENT_ARGS = { + "mie_gruneisen": ("mg_rho0", "mg_c0", "mg_s", "mg_gruneisen", "mg_gruneisen_a", "mg_s2", "mg_s3"), + "jwl": ("jwl_rho0", "jwl_a", "jwl_b", "jwl_r1", "jwl_r2", "jwl_omega"), + "vinet": ("vinet_rho0", "vinet_k0", "vinet_k0p", "vinet_gruneisen", "vinet_gruneisen_a"), +} + + +def test_coefficients_call_arguments_unchanged(): + actual = {f.suffix: tuple(f"{f.prefix}_{s}" for s in f.coefficients_args) for f in _STATE_DEPENDENT} + assert actual == _EXPECTED_COEFFICIENT_ARGS + + +def test_coefficient_args_are_declared_parameters(): + """Every call argument must be a registry parameter, or the generated call passes a name the + user cannot set and no `physical_parameters` field backs.""" + for family in _STATE_DEPENDENT: + declared = {suffix for suffix, _math in family.required + family.optional} + unknown = set(family.coefficients_args) - declared + assert not unknown, f"family {family.suffix} passes undeclared parameter(s) {sorted(unknown)}" diff --git a/toolchain/mfc/params_tests/test_eos_selector.py b/toolchain/mfc/params_tests/test_eos_selector.py index 02554b556..49da05ccf 100644 --- a/toolchain/mfc/params_tests/test_eos_selector.py +++ b/toolchain/mfc/params_tests/test_eos_selector.py @@ -1,16 +1,37 @@ """Guards the per-fluid equation-of-state selector. -The eos_* constants are hand-written in m_constants.fpp because generate_constants_fpp skips -compound registry keys (fluid_pp(1)%eos is not a valid Fortran identifier), so nothing but this -test keeps the Fortran and Python enumerations in step. +The eos_* constants are generated by generate_constants_fpp from EOS_FAMILIES (compound registry +keys carry a "fortran_prefix" so fluid_pp(1)%eos's names become standalone eos_ constants). +Asserting the generator agrees with the registry it was generated from would be circular, so the +tests below check the emitter's naming/values directly and cross-check against actual Fortran usage +instead of re-deriving from the same registry. """ +import re + from ..case_validator import CaseValidator from ..params.definitions import CONSTRAINTS -from ..params.namelist_parser import get_fortran_constants +from ..params.generators.fortran_gen import generate_constants_fpp +from ..params.namelist_parser import get_mfc_root EOS = CONSTRAINTS["fluid_pp(1)%eos"]["names"] +# Values the hand-written m_constants.fpp block carried before this generator emitted them. +_EXPECTED_EOS_CONSTANT_LINES = { + "integer, parameter :: eos_stiffened_gas = 1", + "integer, parameter :: eos_ideal_gas = 2", + "integer, parameter :: eos_mie_gruneisen = 3", + "integer, parameter :: eos_jwl = 4", + "integer, parameter :: eos_vinet = 5", +} + +# eos_* identifiers that appear in src/ but are not registry-derived selector constants: +# eos_coeffs/eos_coefficients are a variable/type name, eos_rk4_steps is an unrelated, +# hand-written RK4 step count (not a family selector, so not in EOS_FAMILIES). +_NON_FAMILY_EOS_IDENTIFIERS = {"eos_coeffs", "eos_coefficients", "eos_rk4_steps"} + +_EOS_IDENTIFIER_RE = re.compile(r"\b(eos_[a-zA-Z0-9_]+)\b") + def _errors(**overrides): params = {"num_fluids": 1, "fluid_pp(1)%gamma": 2.5, "fluid_pp(1)%pi_inf": 0.0} @@ -20,10 +41,27 @@ def _errors(**overrides): return validator.errors -def test_fortran_and_python_enums_agree(): - fortran = get_fortran_constants() - for name, value in EOS.items(): - assert fortran[f"eos_{name}"] == value, f"eos_{name} disagrees between m_constants.fpp and definitions.py" +def test_generate_constants_fpp_emits_one_line_per_eos_family(): + out = generate_constants_fpp() + lines = {line for line in out.splitlines() if "eos_" in line} + assert lines == _EXPECTED_EOS_CONSTANT_LINES + + +def _eos_identifiers_referenced_under_src(): + referenced = set() + src = get_mfc_root() / "src" + for path in list(src.rglob("*.fpp")) + list(src.rglob("*.f90")): + for line in path.read_text().splitlines(): + code = line.split("!", 1)[0] # drop comments; identifiers, not doc prose, matter here + referenced.update(_EOS_IDENTIFIER_RE.findall(code)) + return referenced + + +def test_every_eos_identifier_in_src_is_generated_or_known_non_family(): + generated = {f"eos_{name}" for name in EOS} + referenced = _eos_identifiers_referenced_under_src() + unaccounted = referenced - generated - _NON_FAMILY_EOS_IDENTIFIERS + assert not unaccounted, f"eos_* symbol(s) referenced in src/ but absent from EOS_FAMILIES: {sorted(unaccounted)}" def test_stiffened_gas_accepts_stiffness(): diff --git a/toolchain/mfc/params_tests/test_fortran_gen.py b/toolchain/mfc/params_tests/test_fortran_gen.py index 8991ad6ef..0c8235133 100644 --- a/toolchain/mfc/params_tests/test_fortran_gen.py +++ b/toolchain/mfc/params_tests/test_fortran_gen.py @@ -258,13 +258,13 @@ def test_check_target_raises_on_bad_target(): generate_decls_fpp("bad") -def test_get_generated_files_returns_fifteen(): +def test_get_generated_files_returns_eighteen(): from pathlib import Path from mfc.params.generators.fortran_gen import get_generated_files files = get_generated_files(Path("/build")) - assert len(files) == 15 + assert len(files) == 18 paths = [str(p) for p, _ in files] assert any("pre_process/generated_namelist.fpp" in p for p in paths) assert any("simulation/generated_decls.fpp" in p for p in paths) @@ -298,10 +298,11 @@ def test_get_generated_files_includes_bcast(): "generated_namelist.fpp", "generated_decls.fpp", "generated_constants.fpp", + "generated_eos.fpp", "generated_case_opt_decls.fpp", "generated_bcast.fpp", } - assert len(files) == 15 + assert len(files) == 18 def test_generate_case_opt_decls_fpp(): @@ -711,3 +712,79 @@ def test_mpi_proxy_residue_pins_wall_velocity_and_bc_datatypes(): assert "'bc_x%beg', 'bc_x%end', 'bc_y%beg', 'bc_y%end', 'bc_z%beg', 'bc_z%end']" in src seg = src.split("'bc_z%beg', 'bc_z%end']", 1)[1] assert "MPI_INTEGER" in seg.split("#:endfor")[0], f"{target}: BC codes not MPI_INTEGER" + + +def test_eos_coeff_layers_split_by_writer_count(): + """A field several families write must be family-dispatched; a single-writer field need not be.""" + from mfc.params.eos_families import EOS_FAMILIES + from mfc.params.generators.fortran_gen import _eos_case_fields + + all_fields = {name for family in EOS_FAMILIES for name in family.eos_coeffs} + case_fields = _eos_case_fields() + assert case_fields == {"rho0", "t0", "gruneisen0", "gruneisen_a"} + assert all_fields - case_fields == {"c0", "s", "s2", "s3", "a", "b", "r1", "r2", "k0", "k0p", "mu_max"} + + +def test_generate_eos_fpp_predicates(): + from mfc.params.generators.fortran_gen import generate_eos_fpp + + out = generate_eos_fpp() + assert "#:def EOS_IS_STATE_DEPENDENT(i)\neoss(${i}$) == eos_mie_gruneisen .or. eoss(${i}$) == eos_jwl .or. eoss(${i}$) == eos_vinet\n#:enddef" in out + assert "#:def EOS_HAS_ISENTROPIC_REFERENCE(i)\neoss(${i}$) == eos_jwl .or. eoss(${i}$) == eos_vinet\n#:enddef" in out + # The runtime gruneisen_a term is not a family property and must stay in m_eos.fpp. + assert "gruneisen_a == 0._wp" not in out.split("#:def EOS_INIT_COEFFS")[0] + assert out == generate_eos_fpp() + + +def test_generate_eos_fpp_predicates_fit_on_one_line(): + """The predicates expand inline, so the widest call site must stay inside Fortran's 132 columns.""" + from mfc.params.generators.fortran_gen import generate_eos_fpp + + body = {} + lines = generate_eos_fpp().splitlines() + for n, line in enumerate(lines): + if line.startswith("#:def EOS_") and not line.startswith("#:def EOS_INIT_"): + body[line] = lines[n + 1].replace("${i}$", "i") + assert body, "no predicate macros emitted" + # m_eos.fpp's two call sites: a bare assignment, and one wrapped in the runtime gruneisen_a test. + sites = { + "#:def EOS_IS_STATE_DEPENDENT(i)": (" yes = ", ""), + "#:def EOS_HAS_ISENTROPIC_REFERENCE(i)": (" yes = (", ") .and. eos_coeffs(i)%gruneisen_a == 0._wp"), + } + assert set(body) == set(sites) + for macro, (before, after) in sites.items(): + rendered = len(before) + len(body[macro]) + len(after) + assert rendered <= 132, f"{macro} call site would render {rendered} columns" + + +def test_generate_eos_fpp_init_sources(): + """The three eos_coeffs source kinds: a parameter copy, a bare literal, and a computed call.""" + from mfc.params.generators.fortran_gen import generate_eos_fpp + + out = generate_eos_fpp().replace("${i}$", "i") + assert " eos_coeffs(i)%c0 = fluid_pp(i)%mg_c0\n" in out + assert " eos_coeffs(i)%gruneisen_a = 0._wp\n case (eos_vinet)" in out + assert "eos_coeffs(i)%mu_max = f_hugoniot_compression_limit(fluid_pp(i)%mg_c0, fluid_pp(i)%mg_s, fluid_pp(i)%mg_s2, &\n & fluid_pp(i)%mg_s3)" in out + assert max(len(line) for line in out.splitlines()) <= 132 + # case default: dflt_real everywhere but gruneisen_a, which s_reference_curve reads unconditionally. + default_arm = out.split(" case default\n")[1].splitlines() + assert default_arm[:5] == [ + " eos_coeffs(i)%rho0 = dflt_real", + " eos_coeffs(i)%t0 = dflt_real", + " eos_coeffs(i)%gruneisen0 = dflt_real", + " eos_coeffs(i)%gruneisen_a = 0._wp", + " end select", + ] + + +def test_eos_computed_call_rejects_what_it_cannot_qualify(): + import pytest + + from mfc.params.eos_families import EOS_FAMILIES + from mfc.params.generators.fortran_gen import _eos_computed_call + + mg = next(f for f in EOS_FAMILIES if f.suffix == "mie_gruneisen") + with pytest.raises(ValueError, match="not one of the family's parameters"): + _eos_computed_call(mg, "f_x(mg_c0, jwl_a)", "i") + with pytest.raises(ValueError, match="not a plain"): + _eos_computed_call(mg, "mg_c0 + 1._wp", "i") diff --git a/toolchain/mfc/test_case_validator.py b/toolchain/mfc/test_case_validator.py index d6f6d585f..f9741a9ed 100644 --- a/toolchain/mfc/test_case_validator.py +++ b/toolchain/mfc/test_case_validator.py @@ -529,3 +529,15 @@ def test_accepts_and_requires(self): self.assertRejects({**BASE, **self.VINET, "fluid_pp(1)%mg_s2": 0.1}, "fluid_pp(1)%mg_* are only read when") for k in ("gamma", "pi_inf"): self.assertRejects({**BASE, **self.VINET, f"fluid_pp(1)%{k}": 1.0}, f"fluid_pp(1)%{k} is not read with eos = 'vinet'") + + def test_variable_gruneisen_is_read_by_the_initial_state_check(self): + """vinet_gruneisen_a must reach _check_initial_states_inside_eos, not just eos.vinet_coefficients. + + At rho = rho0 (mu = 0, as test_accepts_and_requires uses) gruneisen_a has no effect by + construction, so that test cannot catch a dropped gruneisen_a in the validator's dispatch. + Here mu = 0.2 != 0 and the Gamma_G = gruneisen0 + gruneisen_a*mu term flips the verdict: + with gruneisen_a wired in this state is accepted; with it dropped (gruneisen_a treated as 0, + the pre-fix behaviour) rho e goes negative and the case is rejected. + """ + case = {**BASE, **self.VINET, "patch_icpp(1)%alpha_rho(1)": 1.2 * self.VINET["fluid_pp(1)%vinet_rho0"], "fluid_pp(1)%vinet_gruneisen_a": -20.0} + self.assertAccepts(case) diff --git a/toolchain/mfc/test_eos_mie_gruneisen.py b/toolchain/mfc/test_eos_mie_gruneisen.py index fc210d9b6..63c194b04 100644 --- a/toolchain/mfc/test_eos_mie_gruneisen.py +++ b/toolchain/mfc/test_eos_mie_gruneisen.py @@ -6,7 +6,20 @@ import pytest -from mfc.eos import eos_coefficients, isentrope_rk4, jwl_coefficients, jwl_reference, mg_reference, reference_isentrope, reference_temperature, sound_speed, temperature, vinet_reference +from mfc.eos import ( + coefficients_from_curve, + eos_coefficients, + isentrope_rk4, + jwl_coefficients, + jwl_reference, + mg_reference, + reference_isentrope, + reference_temperature, + sound_speed, + temperature, + vinet_coefficients, + vinet_reference, +) # Copper-like, no calibrated material: order-of-magnitude values only. RHO0, C0, S, G0 = 8930.0, 3940.0, 1.49, 2.0 @@ -218,3 +231,22 @@ def bulk(r): h = RHO0 * 1e-6 assert bulk(RHO0) == pytest.approx(k0, rel=1e-12) assert _fd(bulk, RHO0, h) / vinet_reference(RHO0, RHO0, k0, k0p)[2] == pytest.approx(k0p, rel=1e-6) + + +def test_vinet_coefficients_uses_the_variable_gruneisen(): + """Gamma_G = Gamma_0 + a mu, mirroring m_eos.fpp's G0 assignment for every state-dependent family. + + The validator used to compose coefficients_from_curve(vinet_reference(...), gruneisen) directly, + a constant Gamma_G that silently ignored vinet_gruneisen_a. Away from rho0 (mu != 0) that old + expression must disagree with vinet_coefficients, or this test could pass without the fix. + """ + k0, k0p, g0, a = 1.4e11, 5.0, 2.0, 0.7 + rho = 1.2 * RHO0 + mu = rho / RHO0 - 1.0 + + gamma, pi, dpi, dgamma = vinet_coefficients(rho, RHO0, k0, k0p, g0, a) + assert gamma == pytest.approx(1.0 / (g0 + a * mu), rel=1e-12) + + old = coefficients_from_curve(rho, vinet_reference(rho, RHO0, k0, k0p), g0) + assert gamma != pytest.approx(old[0], rel=1e-6) + assert pi != pytest.approx(old[1], rel=1e-6)