Skip to content

Generate caller-shaped thermochemistry routines #1919

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
sbryngelson wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Generate caller-shaped thermochemistry routines #1919

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c48bf48
Generate caller-shaped thermochemistry routines
sbryngelson Sep 23, 2026
6c93dec
Make get_mixture_caloric_state a leaf routine
sbryngelson Sep 23, 2026
3a47e7e
Move the caller-shaped routines into m_thermochem_state
sbryngelson Sep 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 10 additions & 1 deletion docs/documentation/thermochemistry.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,16 @@ to parse a mechanism does not imply that the generator supports every feature in

The generated module provides species metadata, caloric and ideal-gas properties,
temperature inversion, net production rates, fused creation/destruction rates,
mixture viscosity and thermal conductivity, and species diffusivities. MFC still
mixture viscosity and thermal conductivity, and species diffusivities. Routines
shaped for MFC's callers return what a call site needs from one state:
`get_mixture_caloric_state` (species cp/R and mixture cp, cv and energy from one
pass over the NASA7 fits), `get_mixavg_transport_state` (molecular weight, mole
fractions, mixture-averaged diffusivities and conductivity, sharing the composition
work) and `get_species_enthalpies_mass`. They live in the companion module
`m_thermochem_state`, generated into the same file, so `m_thermochem`, which nearly
every MFC file reads, is unchanged; nvfortran 24.7 and older crash compiling unrelated
files when it grows by these routines. They reuse the separate routines'
arithmetic; optimized builds may still differ from separate calls at roundoff. MFC still
owns reaction time integration, including alpha-QSS, and spatial transport
discretization. The generator emits single- or double-precision routines and the
selected CPU, OpenACC or OpenMP annotations. Nonchemistry builds retain the existing
Expand Down
28 changes: 7 additions & 21 deletions src/common/m_chemistry.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ module m_chemistry
& get_species_mass_diffusivities_mixavg, gas_constant, get_mixture_molecular_weight, get_mixture_energy_mass, &
& get_mixture_thermal_conductivity_mixavg, get_species_enthalpies_rt, get_mixture_viscosity_mixavg, &
& get_mixture_specific_heat_cp_mass, get_mixture_enthalpy_mass

use m_thermochem_state, only: get_mixavg_transport_state, get_species_enthalpies_mass
use m_global_parameters

implicit none
Expand Down Expand Up @@ -385,14 +385,6 @@ contains
& - eqn_idx%species%beg + 1))
end do

! Calculate molecular weights and mole fractions
call get_mixture_molecular_weight(Ys_L, MW_L)
call get_mixture_molecular_weight(Ys_R, MW_R)
MW_cell = 0.5_wp*(MW_L + MW_R)

call get_mole_fractions(MW_L, Ys_L, Xs_L)
call get_mole_fractions(MW_R, Ys_R, Xs_R)

P_L = q_prim_qp(eqn_idx%E)%sf(x, y, z)
P_R = q_prim_qp(eqn_idx%E)%sf(x + offsets(1), y + offsets(2), z + offsets(3))

Expand All @@ -405,23 +397,17 @@ contains
rho_cell = 0.5_wp*(rho_L + rho_R)
dT_dxi = (T_R - T_L)/grid_spacing

! Get transport properties
call get_species_mass_diffusivities_mixavg(P_L, T_L, Ys_L, mass_diffusivities_mixavg1)
call get_species_mass_diffusivities_mixavg(P_R, T_R, Ys_R, mass_diffusivities_mixavg2)

call get_mixture_thermal_conductivity_mixavg(T_L, Ys_L, lambda_L)
call get_mixture_thermal_conductivity_mixavg(T_R, Ys_R, lambda_R)
! Transport state and species enthalpies [J/kg] of each side
call get_mixavg_transport_state(P_L, T_L, Ys_L, MW_L, Xs_L, mass_diffusivities_mixavg1, lambda_L)
call get_mixavg_transport_state(P_R, T_R, Ys_R, MW_R, Xs_R, mass_diffusivities_mixavg2, lambda_R)
MW_cell = 0.5_wp*(MW_L + MW_R)

call get_species_enthalpies_rt(T_L, h_l)
call get_species_enthalpies_rt(T_R, h_r)
call get_species_enthalpies_mass(T_L, h_l)
call get_species_enthalpies_mass(T_R, h_r)

! Calculate species properties and gradients
$:GPU_LOOP(parallelism='[seq]')
do i = eqn_idx%species%beg, eqn_idx%species%end
h_l(i - eqn_idx%species%beg + 1) = h_l(i - eqn_idx%species%beg + 1) &
& *gas_constant*T_L/molecular_weights(i - eqn_idx%species%beg + 1)
h_r(i - eqn_idx%species%beg + 1) = h_r(i - eqn_idx%species%beg + 1) &
& *gas_constant*T_R/molecular_weights(i - eqn_idx%species%beg + 1)
Xs_cell(i - eqn_idx%species%beg + 1) = 0.5_wp*(Xs_L(i - eqn_idx%species%beg + 1) + Xs_R(i &
& - eqn_idx%species%beg + 1))
h_k(i - eqn_idx%species%beg + 1) = 0.5_wp*(h_l(i - eqn_idx%species%beg + 1) + h_r(i &
Expand Down
9 changes: 3 additions & 6 deletions src/simulation/m_cbc.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,9 @@ module m_cbc
use m_compute_cbc
use m_boundary_primitives, only: f_vel_ramp
use m_constants, only: riemann_solver_hll, model_eqns_gamma_law, recon_type_weno, recon_type_muscl
use m_thermochem, only: get_mixture_energy_mass, get_mixture_specific_heat_cv_mass, get_mixture_specific_heat_cp_mass, &
& gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights, get_species_specific_heats_r, &
use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights, &
& get_mole_fractions
use m_thermochem_state, only: get_mixture_caloric_state

implicit none

Expand Down Expand Up @@ -646,18 +646,15 @@ contains
call get_mixture_molecular_weight(Ys, Mw)
R_gas = gas_constant/Mw
T = pres/rho/R_gas
call get_mixture_specific_heat_cp_mass(T, Ys, Cp)
call get_mixture_energy_mass(T, Ys, e_mix)
call get_mixture_caloric_state(T, Ys, Cp_i, Cp, Cv, e_mix)
E = rho*e_mix + 5.e-1_wp*rho*vel_K_sum
if (chem_params%gamma_method == 1) then
!> gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
call get_mole_fractions(Mw, Ys, Xs)
call get_species_specific_heats_r(T, Cp_i)
Gamma_i(1:num_species) = Cp_i(1:num_species)/(Cp_i(1:num_species) - 1.0_wp)
gamma = sum(Xs(1:num_species)/(Gamma_i(1:num_species) - 1.0_wp))
else if (chem_params%gamma_method == 2) then
!> gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
call get_mixture_specific_heat_cv_mass(T, Ys, Cv)
gamma = 1.0_wp/(Cp/Cv - 1.0_wp)
end if
end if
Expand Down
17 changes: 4 additions & 13 deletions src/simulation/m_riemann_solver_hll.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,8 @@ module m_riemann_solver_hll
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
use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_mixture_specific_heat_cv_mass, &
& get_mixture_energy_mass, get_species_specific_heats_r, get_species_enthalpies_rt, get_mixture_specific_heat_cp_mass, &
& molecular_weights
use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights
use m_thermochem_state, only: get_mixture_caloric_state
use m_riemann_state

implicit none
Expand Down Expand Up @@ -224,8 +223,8 @@ contains
T_L = pres_L/rho_L/R_gas_L
T_R = pres_R/rho_R/R_gas_R

call get_species_specific_heats_r(T_L, Cp_iL)
call get_species_specific_heats_r(T_R, Cp_iR)
call get_mixture_caloric_state(T_L, Ys_L, Cp_iL, Cp_L, Cv_L, E_L)
call get_mixture_caloric_state(T_R, Ys_R, Cp_iR, Cp_R, Cv_R, E_R)

if (chem_params%gamma_method == 1) then
! gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
Expand All @@ -236,20 +235,12 @@ contains
gamma_R = sum(Xs_R(1:num_species)/(Gamma_iR(1:num_species) - 1.0_wp))
else if (chem_params%gamma_method == 2) then
! gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
call get_mixture_specific_heat_cp_mass(T_L, Ys_L, Cp_L)
call get_mixture_specific_heat_cp_mass(T_R, Ys_R, Cp_R)
call get_mixture_specific_heat_cv_mass(T_L, Ys_L, Cv_L)
call get_mixture_specific_heat_cv_mass(T_R, Ys_R, Cv_R)

Gamm_L = Cp_L/Cv_L
gamma_L = 1.0_wp/(Gamm_L - 1.0_wp)
Gamm_R = Cp_R/Cv_R
gamma_R = 1.0_wp/(Gamm_R - 1.0_wp)
end if

call get_mixture_energy_mass(T_L, Ys_L, E_L)
call get_mixture_energy_mass(T_R, Ys_R, E_R)

E_L = rho_L*E_L + 5.e-1*rho_L*vel_L_rms
E_R = rho_R*E_R + 5.e-1*rho_R*vel_R_rms
H_L = (E_L + pres_L)/rho_L
Expand Down
17 changes: 4 additions & 13 deletions src/simulation/m_riemann_solver_hllc.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,8 @@ module m_riemann_solver_hllc
use m_bubbles_EE
use m_surface_tension
use m_chemistry
use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_mixture_specific_heat_cv_mass, &
& get_mixture_energy_mass, get_species_specific_heats_r, get_species_enthalpies_rt, get_mixture_specific_heat_cp_mass, &
& molecular_weights
use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights
use m_thermochem_state, only: get_mixture_caloric_state
use m_riemann_state

implicit none
Expand Down Expand Up @@ -1020,8 +1019,8 @@ contains
T_L = pres_L/rho_L/R_gas_L
T_R = pres_R/rho_R/R_gas_R

call get_species_specific_heats_r(T_L, Cp_iL)
call get_species_specific_heats_r(T_R, Cp_iR)
call get_mixture_caloric_state(T_L, Ys_L, Cp_iL, Cp_L, Cv_L, E_L)
call get_mixture_caloric_state(T_R, Ys_R, Cp_iR, Cp_R, Cv_R, E_R)

if (chem_params%gamma_method == 1) then
!> gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
Expand All @@ -1032,18 +1031,10 @@ contains
gamma_R = sum(Xs_R(1:num_species)/(Gamma_iR(1:num_species) - 1.0_wp))
else if (chem_params%gamma_method == 2) then
!> gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
call get_mixture_specific_heat_cp_mass(T_L, Ys_L, Cp_L)
call get_mixture_specific_heat_cp_mass(T_R, Ys_R, Cp_R)
call get_mixture_specific_heat_cv_mass(T_L, Ys_L, Cv_L)
call get_mixture_specific_heat_cv_mass(T_R, Ys_R, Cv_R)

Gamm_L = Cp_L/Cv_L; Gamm_R = Cp_R/Cv_R
gamma_L = 1.0_wp/(Gamm_L - 1.0_wp); gamma_R = 1.0_wp/(Gamm_R - 1.0_wp)
end if

call get_mixture_energy_mass(T_L, Ys_L, E_L)
call get_mixture_energy_mass(T_R, Ys_R, E_R)

E_L = rho_L*E_L + 5.e-1*rho_L*vel_L_rms
E_R = rho_R*E_R + 5.e-1*rho_R*vel_R_rms
H_L = (E_L + pres_L)/rho_L
Expand Down
16 changes: 4 additions & 12 deletions src/simulation/m_riemann_solver_lf.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@ module m_riemann_solver_lf
use m_global_parameters
use m_variables_conversion
use m_constants, only: riemann_solver_hll, riemann_solver_hllc, riemann_solver_lax_friedrichs
use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_mixture_specific_heat_cv_mass, &
& get_mixture_energy_mass, get_species_specific_heats_r, get_mixture_specific_heat_cp_mass, molecular_weights
use m_thermochem, only: gas_constant, get_mixture_molecular_weight, molecular_weights
use m_thermochem_state, only: get_mixture_caloric_state
use m_riemann_state

implicit none
Expand Down Expand Up @@ -178,8 +178,8 @@ contains
T_L = pres_L/rho_L/R_gas_L
T_R = pres_R/rho_R/R_gas_R

call get_species_specific_heats_r(T_L, Cp_iL)
call get_species_specific_heats_r(T_R, Cp_iR)
call get_mixture_caloric_state(T_L, Ys_L, Cp_iL, Cp_L, Cv_L, E_L)
call get_mixture_caloric_state(T_R, Ys_R, Cp_iR, Cp_R, Cv_R, E_R)

if (chem_params%gamma_method == 1) then
! gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
Expand All @@ -190,20 +190,12 @@ contains
gamma_R = sum(Xs_R(1:num_species)/(Gamma_iR(1:num_species) - 1.0_wp))
else if (chem_params%gamma_method == 2) then
! gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
call get_mixture_specific_heat_cp_mass(T_L, Ys_L, Cp_L)
call get_mixture_specific_heat_cp_mass(T_R, Ys_R, Cp_R)
call get_mixture_specific_heat_cv_mass(T_L, Ys_L, Cv_L)
call get_mixture_specific_heat_cv_mass(T_R, Ys_R, Cv_R)

Gamm_L = Cp_L/Cv_L
gamma_L = 1.0_wp/(Gamm_L - 1.0_wp)
Gamm_R = Cp_R/Cv_R
gamma_R = 1.0_wp/(Gamm_R - 1.0_wp)
end if

call get_mixture_energy_mass(T_L, Ys_L, E_L)
call get_mixture_energy_mass(T_R, Ys_R, E_R)

E_L = rho_L*E_L + 5.e-1*rho_L*vel_L_rms
E_R = rho_R*E_R + 5.e-1*rho_R*vel_R_rms
else
Expand Down
53 changes: 53 additions & 0 deletions toolchain/mfc/test_thermochem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -223,3 +223,56 @@ def test_zero_concentration_falloff(tmp_path, precision):
for i, reaction in enumerate(gas.reactions()):
if isinstance(reaction.rate, ct.FalloffRate):
assert rates[i] == 0


FUSED_DRIVER = """
program fused
use m_thermochem
use m_thermochem_state
implicit none
integer :: ierr
real(KIND) :: t, pressure, cp, cv, e, lambda, mw, cp_f, cv_f, e_f, lambda_f, mw_f
real(KIND), dimension(num_species) :: y, cp0_r, cp0_r_f, x, x_f, diffusion, diffusion_f, h_rt, h
do
read(*,*,iostat=ierr) t, pressure, y
if (ierr /= 0) exit
call get_species_specific_heats_r(t, cp0_r)
call get_mixture_specific_heat_cp_mass(t, y, cp)
call get_mixture_specific_heat_cv_mass(t, y, cv)
call get_mixture_energy_mass(t, y, e)
call get_mixture_caloric_state(t, y, cp0_r_f, cp_f, cv_f, e_f)
call get_mixture_molecular_weight(y, mw)
call get_mole_fractions(mw, y, x)
call get_species_mass_diffusivities_mixavg(pressure, t, y, diffusion)
call get_mixture_thermal_conductivity_mixavg(t, y, lambda)
call get_mixavg_transport_state(pressure, t, y, mw_f, x_f, diffusion_f, lambda_f)
call get_species_enthalpies_rt(t, h_rt)
call get_species_enthalpies_mass(t, h)
! Fused routines share the separate routines' arithmetic, so they agree exactly.
if (any(cp0_r /= cp0_r_f) .or. cp /= cp_f .or. cv /= cv_f .or. e /= e_f) stop 1
if (mw /= mw_f .or. any(x /= x_f) .or. any(diffusion /= diffusion_f) .or. lambda /= lambda_f) stop 2
if (any(h /= h_rt*gas_constant*t/molecular_weights)) stop 3
write(*,'(*(ES25.16E3,1X))') cp_f, cv_f, e_f, lambda_f, diffusion_f, h
end do
end program
"""


@pytest.mark.parametrize("mechanism", MECHANISMS[:2])
def test_fused_routines(tmp_path, mechanism):
"""Caller-shaped routines equal the separate calls bitwise and agree with Cantera."""
gas = ct.Solution(mechanism)
executable = compile_kernel(tmp_path, gas, driver_source=FUSED_DRIVER)
states = list(reference_states(gas))
inputs = "\n".join(" ".join(map(str, [t, p, *y])) for t, p, y in states) + "\n"
result = subprocess.run([str(executable)], input=inputs, capture_output=True, text=True, check=True)
rows = np.array([np.fromstring(line, sep=" ") for line in result.stdout.splitlines()])
assert len(rows) == len(states)
for actual, (t, p, y) in zip(rows, states):
gas.TPY = t, p, y
diffusion = gas.mix_diff_coeffs.copy()
for k in np.flatnonzero(y == 1):
diffusion[k] = gas.binary_diff_coeffs[k, k]
h = gas.standard_enthalpies_RT * ct.gas_constant * t / gas.molecular_weights
expected = np.concatenate(([gas.cp_mass, gas.cv_mass, gas.int_energy_mass, gas.thermal_conductivity], diffusion, h))
np.testing.assert_allclose(actual, expected, rtol=2e-11, atol=1e-10)
Loading
Loading