From c48bf484ebfd019801b469b2bd08e89e09f637a5 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 22 Sep 2026 22:52:35 -0500 Subject: [PATCH 1/3] Generate caller-shaped thermochemistry routines MFC's hot loops chain several generated calls per gas state, each re-evaluating the same NASA7 fits or composition. The generator now also emits routines shaped like those call sites: - get_mixture_caloric_state: species cp/R and mixture cp, cv and internal energy from one pass over the NASA7 fits. - get_mixavg_transport_state: molecular weight, mole fractions, mixture-averaged diffusivities and thermal conductivity, computing the composition once. - get_species_enthalpies_mass: species enthalpies in J/kg. They reuse the separate routines' arithmetic: a new test finds them bitwise equal to the separate calls (gfortran -O0) and in agreement with Cantera. Optimized builds may contract or inline them differently, so results can differ at roundoff. The shared diffusivity and conductivity bodies are expanded in place rather than called, so no device routine is nested deeper than before, which CCE is sensitive to. The HLL, HLLC and LF solvers and CBC now make one caloric call per state instead of up to four, and the mixture-averaged diffusion flux one transport call per side instead of five. The calls stay in the loop bodies. Done with Claude Code. --- docs/documentation/thermochemistry.md | 8 +- src/common/m_chemistry.fpp | 28 ++--- src/simulation/m_cbc.fpp | 10 +- src/simulation/m_riemann_solver_hll.fpp | 15 +-- src/simulation/m_riemann_solver_hllc.fpp | 15 +-- src/simulation/m_riemann_solver_lf.fpp | 15 +-- toolchain/mfc/test_thermochem.py | 52 ++++++++++ toolchain/mfc/thermochem/module.f90.mako | 125 ++++++++++++++++++----- 8 files changed, 176 insertions(+), 92 deletions(-) diff --git a/docs/documentation/thermochemistry.md b/docs/documentation/thermochemistry.md index 7398d0593..d2ceec540 100644 --- a/docs/documentation/thermochemistry.md +++ b/docs/documentation/thermochemistry.md @@ -31,7 +31,13 @@ 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 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 diff --git a/src/common/m_chemistry.fpp b/src/common/m_chemistry.fpp index 687b2a118..f34775f85 100644 --- a/src/common/m_chemistry.fpp +++ b/src/common/m_chemistry.fpp @@ -13,7 +13,7 @@ module m_chemistry & get_creation_destruction_rates, get_mole_fractions, get_species_binary_mass_diffusivities, & & 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 + & get_mixture_specific_heat_cp_mass, get_mixture_enthalpy_mass, get_mixavg_transport_state, get_species_enthalpies_mass use m_global_parameters @@ -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)) @@ -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 & diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index 7bddbdc31..d49ea2957 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -15,9 +15,8 @@ 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, & - & get_mole_fractions + use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_mixture_caloric_state, get_species_enthalpies_rt, & + & molecular_weights, get_mole_fractions implicit none @@ -646,18 +645,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 diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index 969d6996a..39dd8e0f8 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -15,8 +15,7 @@ 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, & + use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_mixture_caloric_state, get_species_enthalpies_rt, & & molecular_weights use m_riemann_state @@ -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. @@ -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 diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index d96c8da7e..b1cadf453 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -18,8 +18,7 @@ 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, & + use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_mixture_caloric_state, get_species_enthalpies_rt, & & molecular_weights use m_riemann_state @@ -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. @@ -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 diff --git a/src/simulation/m_riemann_solver_lf.fpp b/src/simulation/m_riemann_solver_lf.fpp index 82dd93059..d5d78ec1c 100644 --- a/src/simulation/m_riemann_solver_lf.fpp +++ b/src/simulation/m_riemann_solver_lf.fpp @@ -12,8 +12,7 @@ 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, get_mixture_caloric_state, molecular_weights use m_riemann_state implicit none @@ -178,8 +177,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. @@ -190,20 +189,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 diff --git a/toolchain/mfc/test_thermochem.py b/toolchain/mfc/test_thermochem.py index 17a0f5f34..e0b562f66 100644 --- a/toolchain/mfc/test_thermochem.py +++ b/toolchain/mfc/test_thermochem.py @@ -223,3 +223,55 @@ 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 + 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) diff --git a/toolchain/mfc/thermochem/module.f90.mako b/toolchain/mfc/thermochem/module.f90.mako index 0a364e4f6..15856d223 100644 --- a/toolchain/mfc/thermochem/module.f90.mako +++ b/toolchain/mfc/thermochem/module.f90.mako @@ -256,6 +256,33 @@ contains end subroutine get_mixture_energy_mass + !> Species cp/R and mixture cp, cv, e [J/kg] from one NASA7 pass. + subroutine get_mixture_caloric_state(temperature, mass_fractions, cp0_r, cp_mix, cv_mix, e_mix) + + GPU_ROUTINE(get_mixture_caloric_state) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: cp0_r + ${real_type}, intent(out) :: cp_mix, cv_mix, e_mix + + ${real_type}, dimension(${sol.n_species}) :: shifted + + call get_species_specific_heats_r(temperature, cp0_r) + call get_mass_averaged_property(mass_fractions, cp0_r, cp_mix) + cp_mix = cp_mix * gas_constant + + shifted = cp0_r - 1.e0_${kind} + call get_mass_averaged_property(mass_fractions, shifted, cv_mix) + cv_mix = cv_mix * gas_constant + + call get_species_enthalpies_rt(temperature, shifted) + shifted = shifted - 1.e0_${kind} + call get_mass_averaged_property(mass_fractions, shifted, e_mix) + e_mix = e_mix * gas_constant * temperature + + end subroutine get_mixture_caloric_state + subroutine get_species_specific_heats_r(temperature, cp0_r) GPU_ROUTINE(get_species_specific_heats_r) @@ -282,6 +309,21 @@ contains end subroutine get_species_enthalpies_rt + !> Species enthalpies per unit mass [J/kg]. + subroutine get_species_enthalpies_mass(temperature, enthalpies) + + GPU_ROUTINE(get_species_enthalpies_mass) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(out), dimension(${sol.n_species}) :: enthalpies + + call get_species_enthalpies_rt(temperature, enthalpies) + %for i in range(sol.n_species): + enthalpies(${i+1}) = enthalpies(${i+1})*gas_constant*temperature/molecular_weights(${i+1}) + %endfor + + end subroutine get_species_enthalpies_mass + subroutine get_species_entropies_r(temperature, s0_r) GPU_ROUTINE(get_species_entropies_r) @@ -709,6 +751,38 @@ contains end subroutine get_mixture_viscosity_mixavg +<%def name="conductivity_body(result)">\ + call get_species_thermal_conductivities(temperature, conductivities) + + ${result} = 0.5_${kind}*(& + sum(mole_fractions*conductivities) + & + 1/sum(mole_fractions/conductivities)) +\ +<%def name="diffusivities_body(result)">\ + call get_species_binary_mass_diffusivities(temperature, bdiff_ij) + + %for sp in range(sol.n_species): + x_sum(${sp + 1}) = ${cgm(ce.diffusivity_mixture_rule_denom_expr( + sol, sp, Variable("mole_fractions"), Variable("bdiff_ij")))} + %endfor + + %for sp in range(sol.n_species): + denom(${sp + 1}) = x_sum(${sp + 1}) - & + mole_fractions(${sp + 1})/bdiff_ij(${sp + 1}, ${sp + 1}) + %endfor + + %for sp in range(sol.n_species): + if (denom(${sp + 1}) .gt. 0e0_${kind}) then + ${result}(${sp + 1}) = & + (mix_mol_weight - & + mole_fractions(${sp + 1})*molecular_weights(${sp + 1}))& + /(pressure * mix_mol_weight * denom(${sp + 1})) + else + ${result}(${sp + 1}) = & + bdiff_ij(${sp + 1}, ${sp + 1}) / pressure + end if + %endfor +\ subroutine get_mixture_thermal_conductivity_mixavg(temperature, & mass_fractions, mixture_thermal_conductivity_mixavg) @@ -723,12 +797,7 @@ contains call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) call get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) - call get_species_thermal_conductivities(temperature, conductivities) - - mixture_thermal_conductivity_mixavg = 0.5_${kind}*(& - sum(mole_fractions*conductivities) + & - 1/sum(mole_fractions/conductivities)) - +${conductivity_body("mixture_thermal_conductivity_mixavg")} end subroutine get_mixture_thermal_conductivity_mixavg subroutine get_species_mass_diffusivities_mixavg(& @@ -747,30 +816,32 @@ contains call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) call get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) - call get_species_binary_mass_diffusivities(temperature, bdiff_ij) +${diffusivities_body("mass_diffusivities_mixavg")} + end subroutine get_species_mass_diffusivities_mixavg - %for sp in range(sol.n_species): - x_sum(${sp + 1}) = ${cgm(ce.diffusivity_mixture_rule_denom_expr( - sol, sp, Variable("mole_fractions"), Variable("bdiff_ij")))} - %endfor + !> Mixture-averaged transport of one state; shares the composition work. + !> The shared bodies are expanded, not called, to keep the call depth. + subroutine get_mixavg_transport_state(pressure, temperature, mass_fractions, & + mix_mol_weight, mole_fractions, mass_diffusivities_mixavg, & + mixture_thermal_conductivity_mixavg) - %for sp in range(sol.n_species): - denom(${sp + 1}) = x_sum(${sp + 1}) - & - mole_fractions(${sp + 1})/bdiff_ij(${sp + 1}, ${sp + 1}) - %endfor + GPU_ROUTINE(get_mixavg_transport_state) - %for sp in range(sol.n_species): - if (denom(${sp + 1}) .gt. 0e0_${kind}) then - mass_diffusivities_mixavg(${sp + 1}) = & - (mix_mol_weight - & - mole_fractions(${sp + 1})*molecular_weights(${sp + 1}))& - /(pressure * mix_mol_weight * denom(${sp + 1})) - else - mass_diffusivities_mixavg(${sp + 1}) = & - bdiff_ij(${sp + 1}, ${sp + 1}) / pressure - end if - %endfor + ${real_type}, intent(in) :: pressure, temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: mix_mol_weight + ${real_type}, intent(out), dimension(${sol.n_species}) :: mole_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: & + mass_diffusivities_mixavg + ${real_type}, intent(out) :: mixture_thermal_conductivity_mixavg - end subroutine get_species_mass_diffusivities_mixavg + ${real_type}, dimension(${sol.n_species}) :: conductivities, x_sum, denom + ${real_type}, dimension(${sol.n_species}, ${sol.n_species}) :: bdiff_ij + + call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) + call get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) +${diffusivities_body("mass_diffusivities_mixavg")} +${conductivity_body("mixture_thermal_conductivity_mixavg")} + end subroutine get_mixavg_transport_state end module ${module_name} From 6c93dec353e0422e126ee8821078377264524ad2 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 23 Sep 2026 10:33:36 -0500 Subject: [PATCH 2/3] Make get_mixture_caloric_state a leaf routine Frontier CCE (OpenACC and OpenMP) faulted in the HLLC kernel with a GPU memory access violation: get_mixture_caloric_state passed the species cp array on to get_species_specific_heats_r from inside the kernel, the nested-call pattern CCE is known to fault on. The routine now expands the NASA7 cp/R and h/RT expressions and the mass averaging in place and calls nothing, so HLLC's call is shallower than master's get_mixture_energy_mass. The expressions are the separate routines', so the fused-versus-separate bitwise test still passes. Checked: test_thermochem (17 passed), amdflang 23.2.1 OpenMP offload on MI210 --only Chemistry (16 passed), nvfortran 24.1 --gpu mp --mpi build. Done with Claude Code. --- toolchain/mfc/thermochem/module.f90.mako | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/toolchain/mfc/thermochem/module.f90.mako b/toolchain/mfc/thermochem/module.f90.mako index 15856d223..607ddfdf3 100644 --- a/toolchain/mfc/thermochem/module.f90.mako +++ b/toolchain/mfc/thermochem/module.f90.mako @@ -256,7 +256,16 @@ contains end subroutine get_mixture_energy_mass +<%def name="mass_average(result, prop)">\ + ${result} = ( & + %for i in range(sol.n_species): + + inv_molecular_weights(${i+1})*mass_fractions(${i+1}) & + *${prop}(${i+1}) & + %endfor + ) +\ !> Species cp/R and mixture cp, cv, e [J/kg] from one NASA7 pass. + !> A leaf routine: CCE faults on thermochem calls nested below a kernel. subroutine get_mixture_caloric_state(temperature, mass_fractions, cp0_r, cp_mix, cv_mix, e_mix) GPU_ROUTINE(get_mixture_caloric_state) @@ -268,17 +277,21 @@ contains ${real_type}, dimension(${sol.n_species}) :: shifted - call get_species_specific_heats_r(temperature, cp0_r) - call get_mass_averaged_property(mass_fractions, cp0_r, cp_mix) + %for i, sp in enumerate(sol.species()): + cp0_r(${i+1}) = ${cgm(ce.poly_to_expr(sp.thermo, "temperature"))} + %endfor +${mass_average("cp_mix", "cp0_r")} cp_mix = cp_mix * gas_constant shifted = cp0_r - 1.e0_${kind} - call get_mass_averaged_property(mass_fractions, shifted, cv_mix) +${mass_average("cv_mix", "shifted")} cv_mix = cv_mix * gas_constant - call get_species_enthalpies_rt(temperature, shifted) + %for i, sp in enumerate(sol.species()): + shifted(${i+1}) = ${cgm(ce.poly_to_enthalpy_expr(sp.thermo, "temperature"))} + %endfor shifted = shifted - 1.e0_${kind} - call get_mass_averaged_property(mass_fractions, shifted, e_mix) +${mass_average("e_mix", "shifted")} e_mix = e_mix * gas_constant * temperature end subroutine get_mixture_caloric_state From 3a47e7ec7539b3ab62f76fe763b9ab85771e9d25 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 23 Sep 2026 13:38:08 -0500 Subject: [PATCH 3/3] Move the caller-shaped routines into m_thermochem_state NVHPC 23.11 through 24.7 (OpenACC) crashed in the front end (fort1, signal 11) compiling m_compute_levelset, a file this PR does not touch. Reproduced with nvfortran 24.1 and bisected: master's module files plus this PR's m_thermochem.mod alone crash it, and removing any one of the three new routines avoids it. m_global_parameters imports m_thermochem, so nearly every file reads that module; the crash depends on its contents, not simply its size (gri30's much larger module compiles). The three routines now live in a companion module, m_thermochem_state, generated into the same file and imported only by their callers. m_thermochem is byte-identical to master's for h2o2 and gri30. Checked: test_thermochem (17 passed); nvfortran 24.1 --gpu acc --mpi Chemistry build (previously crashed); gfortran CPU and amdflang OpenMP offload on MI210 --only Chemistry (16 passed each). Done with Claude Code. --- docs/documentation/thermochemistry.md | 5 +- src/common/m_chemistry.fpp | 4 +- src/simulation/m_cbc.fpp | 5 +- src/simulation/m_riemann_solver_hll.fpp | 4 +- src/simulation/m_riemann_solver_hllc.fpp | 4 +- src/simulation/m_riemann_solver_lf.fpp | 3 +- toolchain/mfc/test_thermochem.py | 1 + toolchain/mfc/thermochem/module.f90.mako | 218 ++++++++++++++--------- 8 files changed, 147 insertions(+), 97 deletions(-) diff --git a/docs/documentation/thermochemistry.md b/docs/documentation/thermochemistry.md index d2ceec540..12a3bac0f 100644 --- a/docs/documentation/thermochemistry.md +++ b/docs/documentation/thermochemistry.md @@ -36,7 +36,10 @@ 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 reuse the separate routines' +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 diff --git a/src/common/m_chemistry.fpp b/src/common/m_chemistry.fpp index f34775f85..2735c712a 100644 --- a/src/common/m_chemistry.fpp +++ b/src/common/m_chemistry.fpp @@ -13,8 +13,8 @@ module m_chemistry & get_creation_destruction_rates, get_mole_fractions, get_species_binary_mass_diffusivities, & & 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, get_mixavg_transport_state, get_species_enthalpies_mass - + & 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 diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index d49ea2957..705cca23c 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -15,8 +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: gas_constant, get_mixture_molecular_weight, get_mixture_caloric_state, get_species_enthalpies_rt, & - & molecular_weights, get_mole_fractions + 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 diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index 39dd8e0f8..b9ea17387 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -15,8 +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_caloric_state, get_species_enthalpies_rt, & - & 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 diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index b1cadf453..6d49f2931 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -18,8 +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_caloric_state, get_species_enthalpies_rt, & - & 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 diff --git a/src/simulation/m_riemann_solver_lf.fpp b/src/simulation/m_riemann_solver_lf.fpp index d5d78ec1c..5672d955a 100644 --- a/src/simulation/m_riemann_solver_lf.fpp +++ b/src/simulation/m_riemann_solver_lf.fpp @@ -12,7 +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_caloric_state, 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 diff --git a/toolchain/mfc/test_thermochem.py b/toolchain/mfc/test_thermochem.py index e0b562f66..d31b50960 100644 --- a/toolchain/mfc/test_thermochem.py +++ b/toolchain/mfc/test_thermochem.py @@ -228,6 +228,7 @@ def test_zero_concentration_falloff(tmp_path, precision): 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 diff --git a/toolchain/mfc/thermochem/module.f90.mako b/toolchain/mfc/thermochem/module.f90.mako index 607ddfdf3..392366397 100644 --- a/toolchain/mfc/thermochem/module.f90.mako +++ b/toolchain/mfc/thermochem/module.f90.mako @@ -256,46 +256,6 @@ contains end subroutine get_mixture_energy_mass -<%def name="mass_average(result, prop)">\ - ${result} = ( & - %for i in range(sol.n_species): - + inv_molecular_weights(${i+1})*mass_fractions(${i+1}) & - *${prop}(${i+1}) & - %endfor - ) -\ - !> Species cp/R and mixture cp, cv, e [J/kg] from one NASA7 pass. - !> A leaf routine: CCE faults on thermochem calls nested below a kernel. - subroutine get_mixture_caloric_state(temperature, mass_fractions, cp0_r, cp_mix, cv_mix, e_mix) - - GPU_ROUTINE(get_mixture_caloric_state) - - ${real_type}, intent(in) :: temperature - ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions - ${real_type}, intent(out), dimension(${sol.n_species}) :: cp0_r - ${real_type}, intent(out) :: cp_mix, cv_mix, e_mix - - ${real_type}, dimension(${sol.n_species}) :: shifted - - %for i, sp in enumerate(sol.species()): - cp0_r(${i+1}) = ${cgm(ce.poly_to_expr(sp.thermo, "temperature"))} - %endfor -${mass_average("cp_mix", "cp0_r")} - cp_mix = cp_mix * gas_constant - - shifted = cp0_r - 1.e0_${kind} -${mass_average("cv_mix", "shifted")} - cv_mix = cv_mix * gas_constant - - %for i, sp in enumerate(sol.species()): - shifted(${i+1}) = ${cgm(ce.poly_to_enthalpy_expr(sp.thermo, "temperature"))} - %endfor - shifted = shifted - 1.e0_${kind} -${mass_average("e_mix", "shifted")} - e_mix = e_mix * gas_constant * temperature - - end subroutine get_mixture_caloric_state - subroutine get_species_specific_heats_r(temperature, cp0_r) GPU_ROUTINE(get_species_specific_heats_r) @@ -322,21 +282,6 @@ ${mass_average("e_mix", "shifted")} end subroutine get_species_enthalpies_rt - !> Species enthalpies per unit mass [J/kg]. - subroutine get_species_enthalpies_mass(temperature, enthalpies) - - GPU_ROUTINE(get_species_enthalpies_mass) - - ${real_type}, intent(in) :: temperature - ${real_type}, intent(out), dimension(${sol.n_species}) :: enthalpies - - call get_species_enthalpies_rt(temperature, enthalpies) - %for i in range(sol.n_species): - enthalpies(${i+1}) = enthalpies(${i+1})*gas_constant*temperature/molecular_weights(${i+1}) - %endfor - - end subroutine get_species_enthalpies_mass - subroutine get_species_entropies_r(temperature, s0_r) GPU_ROUTINE(get_species_entropies_r) @@ -764,14 +709,44 @@ ${mass_average("e_mix", "shifted")} end subroutine get_mixture_viscosity_mixavg -<%def name="conductivity_body(result)">\ + subroutine get_mixture_thermal_conductivity_mixavg(temperature, & + mass_fractions, mixture_thermal_conductivity_mixavg) + + GPU_ROUTINE(get_mixture_thermal_conductivity_mixavg) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: mixture_thermal_conductivity_mixavg + + ${real_type} :: mix_mol_weight + ${real_type}, dimension(${sol.n_species}) :: mole_fractions, conductivities + + call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) + call get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) call get_species_thermal_conductivities(temperature, conductivities) - ${result} = 0.5_${kind}*(& + mixture_thermal_conductivity_mixavg = 0.5_${kind}*(& sum(mole_fractions*conductivities) + & 1/sum(mole_fractions/conductivities)) -\ -<%def name="diffusivities_body(result)">\ + + end subroutine get_mixture_thermal_conductivity_mixavg + + subroutine get_species_mass_diffusivities_mixavg(& + pressure, temperature, mass_fractions, mass_diffusivities_mixavg) + + GPU_ROUTINE(get_species_mass_diffusivities_mixavg) + + ${real_type}, intent(in) :: pressure, temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: & + mass_diffusivities_mixavg + + ${real_type} :: mix_mol_weight + ${real_type}, dimension(${sol.n_species}) :: mole_fractions, x_sum, denom + ${real_type}, dimension(${sol.n_species}, ${sol.n_species}) :: bdiff_ij + + call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) + call get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) call get_species_binary_mass_diffusivities(temperature, bdiff_ij) %for sp in range(sol.n_species): @@ -786,51 +761,80 @@ ${mass_average("e_mix", "shifted")} %for sp in range(sol.n_species): if (denom(${sp + 1}) .gt. 0e0_${kind}) then - ${result}(${sp + 1}) = & + mass_diffusivities_mixavg(${sp + 1}) = & (mix_mol_weight - & mole_fractions(${sp + 1})*molecular_weights(${sp + 1}))& /(pressure * mix_mol_weight * denom(${sp + 1})) else - ${result}(${sp + 1}) = & + mass_diffusivities_mixavg(${sp + 1}) = & bdiff_ij(${sp + 1}, ${sp + 1}) / pressure end if %endfor -\ - subroutine get_mixture_thermal_conductivity_mixavg(temperature, & - mass_fractions, mixture_thermal_conductivity_mixavg) - GPU_ROUTINE(get_mixture_thermal_conductivity_mixavg) + end subroutine get_species_mass_diffusivities_mixavg + +end module ${module_name} + +!> Routines shaped like MFC's call sites, kept out of the module above: +!> nearly every MFC file reads it, and nvfortran <= 24.7 crashes (fort1) +!> compiling unrelated files once it grows by these routines. +module ${module_name}_state + + use ${module_name} + + implicit none + + private + public :: get_mixture_caloric_state, get_species_enthalpies_mass, get_mixavg_transport_state + +contains + + !> Species cp/R and mixture cp, cv, e [J/kg] from one NASA7 pass. + !> A leaf routine: CCE faults on thermochem calls nested below a kernel. + subroutine get_mixture_caloric_state(temperature, mass_fractions, cp0_r, cp_mix, cv_mix, e_mix) + + GPU_ROUTINE(get_mixture_caloric_state) ${real_type}, intent(in) :: temperature ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions - ${real_type}, intent(out) :: mixture_thermal_conductivity_mixavg + ${real_type}, intent(out), dimension(${sol.n_species}) :: cp0_r + ${real_type}, intent(out) :: cp_mix, cv_mix, e_mix - ${real_type} :: mix_mol_weight - ${real_type}, dimension(${sol.n_species}) :: mole_fractions, conductivities + ${real_type}, dimension(${sol.n_species}) :: shifted - call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) - call get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) -${conductivity_body("mixture_thermal_conductivity_mixavg")} - end subroutine get_mixture_thermal_conductivity_mixavg + %for i, sp in enumerate(sol.species()): + cp0_r(${i+1}) = ${cgm(ce.poly_to_expr(sp.thermo, "temperature"))} + %endfor +${mass_average("cp_mix", "cp0_r")} + cp_mix = cp_mix * gas_constant - subroutine get_species_mass_diffusivities_mixavg(& - pressure, temperature, mass_fractions, mass_diffusivities_mixavg) + shifted = cp0_r - 1.e0_${kind} +${mass_average("cv_mix", "shifted")} + cv_mix = cv_mix * gas_constant - GPU_ROUTINE(get_species_mass_diffusivities_mixavg) + %for i, sp in enumerate(sol.species()): + shifted(${i+1}) = ${cgm(ce.poly_to_enthalpy_expr(sp.thermo, "temperature"))} + %endfor + shifted = shifted - 1.e0_${kind} +${mass_average("e_mix", "shifted")} + e_mix = e_mix * gas_constant * temperature - ${real_type}, intent(in) :: pressure, temperature - ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions - ${real_type}, intent(out), dimension(${sol.n_species}) :: & - mass_diffusivities_mixavg + end subroutine get_mixture_caloric_state - ${real_type} :: mix_mol_weight - ${real_type}, dimension(${sol.n_species}) :: mole_fractions, x_sum, denom - ${real_type}, dimension(${sol.n_species}, ${sol.n_species}) :: bdiff_ij + !> Species enthalpies per unit mass [J/kg]. + subroutine get_species_enthalpies_mass(temperature, enthalpies) - call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) - call get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) -${diffusivities_body("mass_diffusivities_mixavg")} - end subroutine get_species_mass_diffusivities_mixavg + GPU_ROUTINE(get_species_enthalpies_mass) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(out), dimension(${sol.n_species}) :: enthalpies + + call get_species_enthalpies_rt(temperature, enthalpies) + %for i in range(sol.n_species): + enthalpies(${i+1}) = enthalpies(${i+1})*gas_constant*temperature/molecular_weights(${i+1}) + %endfor + + end subroutine get_species_enthalpies_mass !> Mixture-averaged transport of one state; shares the composition work. !> The shared bodies are expanded, not called, to keep the call depth. @@ -857,4 +861,44 @@ ${diffusivities_body("mass_diffusivities_mixavg")} ${conductivity_body("mixture_thermal_conductivity_mixavg")} end subroutine get_mixavg_transport_state -end module ${module_name} +end module ${module_name}_state +<%def name="mass_average(result, prop)">\ + ${result} = ( & + %for i in range(sol.n_species): + + inv_molecular_weights(${i+1})*mass_fractions(${i+1}) & + *${prop}(${i+1}) & + %endfor + ) +\ +<%def name="conductivity_body(result)">\ + call get_species_thermal_conductivities(temperature, conductivities) + + ${result} = 0.5_${kind}*(& + sum(mole_fractions*conductivities) + & + 1/sum(mole_fractions/conductivities)) +\ +<%def name="diffusivities_body(result)">\ + call get_species_binary_mass_diffusivities(temperature, bdiff_ij) + + %for sp in range(sol.n_species): + x_sum(${sp + 1}) = ${cgm(ce.diffusivity_mixture_rule_denom_expr( + sol, sp, Variable("mole_fractions"), Variable("bdiff_ij")))} + %endfor + + %for sp in range(sol.n_species): + denom(${sp + 1}) = x_sum(${sp + 1}) - & + mole_fractions(${sp + 1})/bdiff_ij(${sp + 1}, ${sp + 1}) + %endfor + + %for sp in range(sol.n_species): + if (denom(${sp + 1}) .gt. 0e0_${kind}) then + ${result}(${sp + 1}) = & + (mix_mol_weight - & + mole_fractions(${sp + 1})*molecular_weights(${sp + 1}))& + /(pressure * mix_mol_weight * denom(${sp + 1})) + else + ${result}(${sp + 1}) = & + bdiff_ij(${sp + 1}, ${sp + 1}) / pressure + end if + %endfor +\