diff --git a/src/common/include/inline_eos.fpp b/src/common/include/inline_eos.fpp new file mode 100644 index 000000000..aed463362 --- /dev/null +++ b/src/common/include/inline_eos.fpp @@ -0,0 +1,38 @@ +!> +!! @file inline_eos.fpp +!! @brief Helpers for modules that compile a private copy of the hot EOS routines (see inline_eos_procs.fpp). +!! +!! A kernel module does three things: +!! 1. $:EOS_INLINE_USE('m_eos') in place of its "use m_eos": the public routines are renamed to _ext so they +!! do not clash with the local copies. +!! 2. $:EOS_INLINE_PRIVATE() after "implicit none": the copies stay private, so they are never re-exported. +!! 3. $:EOS_INLINE_PROCS() right after "contains" (the macro comes from inline_eos_procs.fpp, included below). +!! Call sites stay as they are and resolve to the local copy. + +#:include 'inline_eos_procs.fpp' + +#:set EOS_INLINE_NAMES = ['f_isentrope_exponent', 'f_isentrope_pressure', 's_phase_coefficients', 'f_bulk_modulus', & + & 'f_relativistic_enthalpy', 's_compute_mixture_coefficients', 's_compute_speed_of_sound', & + & 's_compute_speed_of_sound_avg', 's_compute_energy', 'f_is_state_dependent', 'f_pressure', 's_phase_internal_energy', & + & 's_phase_bulk_modulus', 's_phase_pressure_on_isentrope', 's_phase_temperature', 's_compute_mixture_coefficients_dt'] + +#:def EOS_INLINE_USE(module, only=None) + #:if only is None + use ${module}$, & + #:else + #! the copies' own dependencies, which a plain "use m_eos" would bring in + use ${module}$, only: ${''.join(x + ', ' for x in only)}$s_eos_coefficients, f_has_isentropic_reference, s_reference_curve, & + & s_rk4, & + #:endif + #:for n in EOS_INLINE_NAMES[:-1] + & ${n}$_ext => ${n}$, & + #:endfor + & ${EOS_INLINE_NAMES[-1]}$_ext => ${EOS_INLINE_NAMES[-1]}$ + use m_constants, only: ode_isentrope, ode_reference_temperature +#:enddef + +#:def EOS_INLINE_PRIVATE() + #:for n in EOS_INLINE_NAMES + private :: ${n}$ + #:endfor +#:enddef diff --git a/src/common/include/inline_eos_procs.fpp b/src/common/include/inline_eos_procs.fpp new file mode 100644 index 000000000..df7234322 --- /dev/null +++ b/src/common/include/inline_eos_procs.fpp @@ -0,0 +1,440 @@ +!> +!! @file inline_eos_procs.fpp +!! @brief Hot equation-of-state routines, one source for every module that calls them from a GPU kernel. +!! +!! nvfortran 26.5 inlines a routine into an OpenACC kernel only when it is defined in the same file; a routine from +!! another module stays a real device CALL (-Minline=lib inlines host-side only). These routines are called per face +!! and per cell, so they are kept in the Fypp macro EOS_INLINE_PROCS below: m_eos expands it for the public definitions +!! and each kernel module expands it again inside its own contains section for a private copy. The kernel module must rename the public ones on use +!! (@:EOS_INLINE_USE) and mark its copies private ($:EOS_INLINE_PRIVATE()), both from inline_eos.fpp. +!! Edit the routines here only. + +#:def EOS_INLINE_PROCS() + !> 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 + + !> 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 + + !> 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 + + !> Mixture coefficients of one state. Under bubbles_euler with num_fluids == 1 the sole advection slot aliases the void fraction + !! (eqn_idx%alf == eqn_idx%adv%end), so alpha is not a composition there and the coefficients are the liquid's. Clipping stays + !! with callers; it differs between solvers and cannot coincide with that case, as mpp_lim requires num_fluids > 1. + subroutine s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K) + + $:GPU_ROUTINE(function_name='s_compute_mixture_coefficients', parallelism='[seq]', cray_inline=True) + + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in) :: alpha_rho_K, alpha_K + #:else + real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K + #:endif + real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K + real(wp) :: gamma_i, pi_inf_i, dpi_i, dgamma_i + real(wp) :: rho_i, alpha_i, alpha_rho_i + integer :: i !< Loop iterator over fluids + + ! The bubbly closure is written for one carrier liquid, which keeps its own coefficients + ! undiluted: Gamma_l*p_l = (E - rho|u|^2/2)/(1 - alf) - Pi_inf_l, the void entering only through + ! the (1 - alf) that s_compute_pressure applies. There is nothing to sum - the last advection + ! slot is the void, not a material - and the checker holds num_fluids <= 2 here. + if (bubbles_euler) then + rho_K = alpha_rho_K(1) + gamma_K = gammas(1) + pi_inf_K = pi_infs(1) + ! Energy per unit volume, as below: alpha_rho_K(1) is the liquid partial density + qv_K = alpha_rho_K(1)*qvs(1) + else + rho_K = 0._wp + gamma_K = 0._wp + pi_inf_K = 0._wp + qv_K = 0._wp + + $:GPU_LOOP(parallelism='[seq]') + do i = 1, num_fluids + rho_K = rho_K + alpha_rho_K(i) + alpha_rho_i = alpha_rho_K(i) + alpha_i = alpha_K(i) + call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i) + gamma_K = gamma_K + alpha_K(i)*gamma_i + pi_inf_K = pi_inf_K + alpha_K(i)*pi_inf_i + qv_K = qv_K + alpha_rho_K(i)*qvs(i) + end do + end if + + end subroutine s_compute_mixture_coefficients + + !> 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) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: pres, rho, gamma, pi_inf + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in) :: adv + #:else + real(wp), dimension(num_fluids), intent(in) :: adv + #:endif + real(wp), intent(out) :: c + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in), optional :: alpha_rho + #:else + real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho + #:endif + real(wp) :: alf !< Subgrid void fraction; dilute by construction + real(wp) :: blkmod_q, alpha_q, alpha_rho_q, gamma_q, pi_inf_q + integer :: q + + if (chemistry) then ! Reacting mixture sound speed + c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho) + else if (relativity) then ! Relativistic sound speed, whose enthalpy is 1 + (Gamma + 1)p/rho + c = sqrt((1._wp + 1._wp/gamma)*pres/rho/f_relativistic_enthalpy(pres, rho, gamma)) + else + ! Every case below is a bulk modulus over a density. The equation of state enters + ! only through f_bulk_modulus; the cases differ in how the phases are mixed. + if (any_state_dependent_eos .and. present(alpha_rho)) then ! frozen mixing: each phase's modulus at its own density + c = 0._wp + $:GPU_LOOP(parallelism='[seq]') + do q = 1, num_fluids + alpha_q = adv(q) + alpha_rho_q = alpha_rho(q) + call s_phase_bulk_modulus(pres, alpha_q, alpha_rho_q, q, blkmod_q) + if (alt_soundspeed) then + c = c + adv(q)/blkmod_q + else + c = c + adv(q)*blkmod_q + end if + end do + if (alt_soundspeed) then + c = 1._wp/(rho*c) + else + c = c/rho + end if + else if (alt_soundspeed) then ! Wood's law: volume-weighted harmonic mean + c = 0._wp + $:GPU_LOOP(parallelism='[seq]') + do q = 1, num_fluids + gamma_q = gammas(q) + pi_inf_q = pi_infs(q) + c = c + adv(q)/f_bulk_modulus(pres, gamma_q, pi_inf_q) + end do + c = 1._wp/(rho*c) + else if (model_eqns == model_eqns_6eq) then ! volume-weighted arithmetic mean + c = 0._wp + $:GPU_LOOP(parallelism='[seq]') + do q = 1, num_fluids + gamma_q = gammas(q) + pi_inf_q = pi_infs(q) + c = c + adv(q)*f_bulk_modulus(pres, gamma_q, pi_inf_q) + end do + c = c/rho + else ! the mixture coefficients already carry the mixing + c = f_bulk_modulus(pres, gamma, pi_inf)/rho + + ! Subgrid bubbles: c = c_l/(1 - alf), the carrier-liquid speed with an O(alf) void + ! correction. alf is dilute by construction; near one means a wrong index or an + ! out-of-regime case, which the toolchain warns about at case load (#1793). + if (model_eqns == model_eqns_5eq .and. bubbles_euler .and. .not. (mpp_lim .and. num_fluids > 1)) then + alf = adv(num_fluids) + c = c/(1._wp - alf) + end if + end if + + if (mixture_err .and. c < 0._wp) then + c = 100._wp*sgm_eps + else + c = sqrt(c) + end if + end if + + end subroutine s_compute_speed_of_sound + + !> Speed of sound of an interface-averaged state. An average of two states is not a state - its enthalpy is not the one its + !! pressure and density imply - so the caller supplies H, |u|^2 and qv. Only the enthalpy-reading branches differ from + !! s_compute_speed_of_sound; keep the condition below in step with the branch list there. + subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c, alpha_rho) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in) :: adv + #:else + real(wp), dimension(num_fluids), intent(in) :: adv + #:endif + real(wp), intent(out) :: c + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in), optional :: alpha_rho + #:else + real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho + #:endif + + if (chemistry) then ! Reacting mixture sound speed + if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then + c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - H)) + else + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho) + end if + else if (relativity) then ! Relativistic sound speed + c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H) + else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler) & + & .or. any_state_dependent_eos) then + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho) + else ! Stiffened-gas mixture, the one branch where the averaged enthalpy survives + c = (H - 5.e-1*vel_sum - qv/rho)/gamma + + if (mixture_err .and. c < 0._wp) then + c = 100._wp*sgm_eps + else + c = sqrt(c) + end if + end if + + end subroutine s_compute_speed_of_sound_avg + + !> Total energy per unit volume, thermodynamic terms only. Callers add magnetic and elastic energy, which are not + !! equation-of-state terms. The chemistry and relativistic branches use a different relation and stay open-coded. + subroutine s_compute_energy(pres, alpha_rho_K, alpha_K, vel_sum, E) + + $:GPU_ROUTINE(function_name='s_compute_energy', parallelism='[seq]', cray_inline=True) + + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in) :: alpha_rho_K, alpha_K + #:else + real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K + #:endif + real(wp), intent(in) :: pres, vel_sum + real(wp), intent(out) :: E + real(wp) :: rho, gamma, pi_inf, qv + + call s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho, gamma, pi_inf, qv) + + ! E = dyn_p + (1 - alf)(gamma p + pi_inf) + qv. Only the liquid's internal energy is diluted; qv is already an energy + ! density. + E = gamma*pres + pi_inf + if (bubbles_euler) E = E*(1._wp - alpha_K(num_fluids)) + E = E + qv + 5.e-1_wp*rho*vel_sum + + end subroutine s_compute_energy + + !> 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 + + #:if defined('EOS_INLINE_NAMES') + #:if MFC_CASE_OPTIMIZATION and not any_state_dependent_eos + ! Module-local copy in a case-optimized build without a state-dependent EOS: s_initialize_eos_module + ! aborts unless every fluid is state-independent, so this is a constant and the state-dependent + ! branches of the phase routines drop out of the kernels. + yes = .false. + #:else + ! Module-local copy: the family table is generated for m_eos only, so ask the public routine. + yes = f_is_state_dependent_ext(i) + #:endif + #:else + yes = ${EOS_IS_STATE_DEPENDENT('i')}$ + #:endif + + end function f_is_state_dependent + + !> 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 + + !> 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 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 + + !> Time derivative of the mixture coefficients, mirroring s_compute_mixture_coefficients. + subroutine s_compute_mixture_coefficients_dt(dalpha_rho_dt, dadv_dt, alpha_rho, adv, drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt) + + $:GPU_ROUTINE(function_name='s_compute_mixture_coefficients_dt', parallelism='[seq]', cray_inline=True) + + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv + #:else + real(wp), dimension(num_fluids), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv + #:endif + real(wp), intent(out) :: drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt + real(wp) :: rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i, alpha_i, alpha_rho_i + integer :: i !< Loop iterator over fluids + + dgamma_dt = 0._wp + dpi_inf_dt = 0._wp + dqv_dt = 0._wp + + if (num_fluids == 1 .and. bubbles_euler) then + ! Fluid 1's coefficients are constants here, so only rho varies. + drho_dt = dalpha_rho_dt(1) + else + drho_dt = 0._wp + + $:GPU_LOOP(parallelism='[seq]') + do i = 1, num_fluids + drho_dt = drho_dt + dalpha_rho_dt(i) + alpha_rho_i = alpha_rho(i) + alpha_i = adv(i) + call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i) + ! d(alpha X(rho_i))/dt with rho_i = alpha_rho/alpha; the alpha in dX/dt cancels + dgamma_dt = dgamma_dt + dadv_dt(i)*gamma_i + dgamma_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i)) + dpi_inf_dt = dpi_inf_dt + dadv_dt(i)*pi_inf_i + dpi_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i)) + dqv_dt = dqv_dt + dalpha_rho_dt(i)*qvs(i) + end do + end if + + end subroutine s_compute_mixture_coefficients_dt +#:enddef diff --git a/src/common/m_eos.fpp b/src/common/m_eos.fpp index 16e126297..7e73e6193 100644 --- a/src/common/m_eos.fpp +++ b/src/common/m_eos.fpp @@ -3,6 +3,7 @@ !! @brief Contains module m_eos #:include 'macros.fpp' +#:include 'inline_eos_procs.fpp' #:include 'case.fpp' #:include 'generated_eos.fpp' @@ -41,7 +42,8 @@ module m_eos & s_compute_speed_of_sound_avg, 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_mixture_temperature, & & 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 + & s_phase_density_on_isentrope, s_phase_internal_energy, s_phase_bulk_modulus, s_eos_coefficients, s_compute_energy, & + & f_has_isentropic_reference, s_reference_curve, s_rk4 contains @@ -183,18 +185,6 @@ contains 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. @@ -262,30 +252,10 @@ contains 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 - + ! Hot EOS routines: the source is the Fypp macro EOS_INLINE_PROCS (include/inline_eos_procs.fpp), so modules that call them from + ! GPU kernels can compile a module-local copy (nvfortran 26.5 inlines same-file routines into a kernel, never + ! routines from another file). This is the public definition. + $:EOS_INLINE_PROCS() !> 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) @@ -325,28 +295,6 @@ contains end function f_mixture_temperature - !> 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) @@ -421,53 +369,6 @@ contains 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) @@ -491,290 +392,4 @@ contains 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 - - !> Mixture coefficients of one state. Under bubbles_euler with num_fluids == 1 the sole advection slot aliases the void fraction - !! (eqn_idx%alf == eqn_idx%adv%end), so alpha is not a composition there and the coefficients are the liquid's. Clipping stays - !! with callers; it differs between solvers and cannot coincide with that case, as mpp_lim requires num_fluids > 1. - subroutine s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K) - - $:GPU_ROUTINE(function_name='s_compute_mixture_coefficients', parallelism='[seq]', cray_inline=True) - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: alpha_rho_K, alpha_K - #:else - real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K - #:endif - real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K - real(wp) :: gamma_i, pi_inf_i, dpi_i, dgamma_i - real(wp) :: rho_i, alpha_i, alpha_rho_i - integer :: i !< Loop iterator over fluids - - ! The bubbly closure is written for one carrier liquid, which keeps its own coefficients - ! undiluted: Gamma_l*p_l = (E - rho|u|^2/2)/(1 - alf) - Pi_inf_l, the void entering only through - ! the (1 - alf) that s_compute_pressure applies. There is nothing to sum - the last advection - ! slot is the void, not a material - and the checker holds num_fluids <= 2 here. - if (bubbles_euler) then - rho_K = alpha_rho_K(1) - gamma_K = gammas(1) - pi_inf_K = pi_infs(1) - ! Energy per unit volume, as below: alpha_rho_K(1) is the liquid partial density - qv_K = alpha_rho_K(1)*qvs(1) - else - rho_K = 0._wp - gamma_K = 0._wp - pi_inf_K = 0._wp - qv_K = 0._wp - - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - rho_K = rho_K + alpha_rho_K(i) - alpha_rho_i = alpha_rho_K(i) - alpha_i = alpha_K(i) - call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i) - gamma_K = gamma_K + alpha_K(i)*gamma_i - pi_inf_K = pi_inf_K + alpha_K(i)*pi_inf_i - qv_K = qv_K + alpha_rho_K(i)*qvs(i) - end do - end if - - end subroutine s_compute_mixture_coefficients - - !> Time derivative of the mixture coefficients, mirroring s_compute_mixture_coefficients. - subroutine s_compute_mixture_coefficients_dt(dalpha_rho_dt, dadv_dt, alpha_rho, adv, drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt) - - $:GPU_ROUTINE(function_name='s_compute_mixture_coefficients_dt', parallelism='[seq]', cray_inline=True) - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv - #:else - real(wp), dimension(num_fluids), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv - #:endif - real(wp), intent(out) :: drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt - real(wp) :: rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i, alpha_i, alpha_rho_i - integer :: i !< Loop iterator over fluids - - dgamma_dt = 0._wp - dpi_inf_dt = 0._wp - dqv_dt = 0._wp - - if (num_fluids == 1 .and. bubbles_euler) then - ! Fluid 1's coefficients are constants here, so only rho varies. - drho_dt = dalpha_rho_dt(1) - else - drho_dt = 0._wp - - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - drho_dt = drho_dt + dalpha_rho_dt(i) - alpha_rho_i = alpha_rho(i) - alpha_i = adv(i) - call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i) - ! d(alpha X(rho_i))/dt with rho_i = alpha_rho/alpha; the alpha in dX/dt cancels - dgamma_dt = dgamma_dt + dadv_dt(i)*gamma_i + dgamma_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i)) - dpi_inf_dt = dpi_inf_dt + dadv_dt(i)*pi_inf_i + dpi_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i)) - dqv_dt = dqv_dt + dalpha_rho_dt(i)*qvs(i) - end do - end if - - end subroutine s_compute_mixture_coefficients_dt - - !> 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) - - $:GPU_ROUTINE(parallelism='[seq]') - - real(wp), intent(in) :: pres, rho, gamma, pi_inf - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: adv - #:else - real(wp), dimension(num_fluids), intent(in) :: adv - #:endif - real(wp), intent(out) :: c - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in), optional :: alpha_rho - #:else - real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho - #:endif - real(wp) :: alf !< Subgrid void fraction; dilute by construction - real(wp) :: blkmod_q, alpha_q, alpha_rho_q, gamma_q, pi_inf_q - integer :: q - - if (chemistry) then ! Reacting mixture sound speed - c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho) - else if (relativity) then ! Relativistic sound speed, whose enthalpy is 1 + (Gamma + 1)p/rho - c = sqrt((1._wp + 1._wp/gamma)*pres/rho/f_relativistic_enthalpy(pres, rho, gamma)) - else - ! Every case below is a bulk modulus over a density. The equation of state enters - ! only through f_bulk_modulus; the cases differ in how the phases are mixed. - if (any_state_dependent_eos .and. present(alpha_rho)) then ! frozen mixing: each phase's modulus at its own density - c = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do q = 1, num_fluids - alpha_q = adv(q) - alpha_rho_q = alpha_rho(q) - call s_phase_bulk_modulus(pres, alpha_q, alpha_rho_q, q, blkmod_q) - if (alt_soundspeed) then - c = c + adv(q)/blkmod_q - else - c = c + adv(q)*blkmod_q - end if - end do - if (alt_soundspeed) then - c = 1._wp/(rho*c) - else - c = c/rho - end if - else if (alt_soundspeed) then ! Wood's law: volume-weighted harmonic mean - c = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do q = 1, num_fluids - gamma_q = gammas(q) - pi_inf_q = pi_infs(q) - c = c + adv(q)/f_bulk_modulus(pres, gamma_q, pi_inf_q) - end do - c = 1._wp/(rho*c) - else if (model_eqns == model_eqns_6eq) then ! volume-weighted arithmetic mean - c = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do q = 1, num_fluids - gamma_q = gammas(q) - pi_inf_q = pi_infs(q) - c = c + adv(q)*f_bulk_modulus(pres, gamma_q, pi_inf_q) - end do - c = c/rho - else ! the mixture coefficients already carry the mixing - c = f_bulk_modulus(pres, gamma, pi_inf)/rho - - ! Subgrid bubbles: c = c_l/(1 - alf), the carrier-liquid speed with an O(alf) void - ! correction. alf is dilute by construction; near one means a wrong index or an - ! out-of-regime case, which the toolchain warns about at case load (#1793). - if (model_eqns == model_eqns_5eq .and. bubbles_euler .and. .not. (mpp_lim .and. num_fluids > 1)) then - alf = adv(num_fluids) - c = c/(1._wp - alf) - end if - end if - - if (mixture_err .and. c < 0._wp) then - c = 100._wp*sgm_eps - else - c = sqrt(c) - end if - end if - - end subroutine s_compute_speed_of_sound - - !> Speed of sound of an interface-averaged state. An average of two states is not a state - its enthalpy is not the one its - !! pressure and density imply - so the caller supplies H, |u|^2 and qv. Only the enthalpy-reading branches differ from - !! s_compute_speed_of_sound; keep the condition below in step with the branch list there. - subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c, alpha_rho) - - $:GPU_ROUTINE(parallelism='[seq]') - - real(wp), intent(in) :: pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: adv - #:else - real(wp), dimension(num_fluids), intent(in) :: adv - #:endif - real(wp), intent(out) :: c - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in), optional :: alpha_rho - #:else - real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho - #:endif - - if (chemistry) then ! Reacting mixture sound speed - if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then - c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - H)) - else - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho) - end if - else if (relativity) then ! Relativistic sound speed - c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H) - else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler) & - & .or. any_state_dependent_eos) then - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho) - else ! Stiffened-gas mixture, the one branch where the averaged enthalpy survives - c = (H - 5.e-1*vel_sum - qv/rho)/gamma - - if (mixture_err .and. c < 0._wp) then - c = 100._wp*sgm_eps - else - c = sqrt(c) - end if - end if - - end subroutine s_compute_speed_of_sound_avg - end module m_eos diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 87c3d2654..ef6fd1f14 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -3,13 +3,14 @@ !! @brief Contains module m_variables_conversion #:include 'macros.fpp' +#:include 'inline_eos.fpp' #:include 'case.fpp' !> @brief Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation module m_variables_conversion use m_derived_types - use m_eos + $:EOS_INLINE_USE('m_eos') use m_global_parameters use m_mpi_proxy use m_helper_basic @@ -21,14 +22,15 @@ module m_variables_conversion implicit none + $:EOS_INLINE_PRIVATE() + private public :: s_initialize_variables_conversion_module, s_initialize_pb, s_initialize_mv, s_convert_to_mixture_variables, & & 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, & - & s_compute_mixture_coefficients_dt, s_compute_speed_of_sound_avg, s_compute_fast_magnetosonic_speed, f_elastic_energy, & - & f_hypoelastic_energy, s_finalize_variables_conversion_module, gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps + & s_compute_species_fraction, s_compute_fast_magnetosonic_speed, f_elastic_energy, 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 @@ -49,6 +51,8 @@ module m_variables_conversion contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Dispatch to the s_convert_mixture_to_mixture_variables and s_convert_species_to_mixture_variables subroutines. Replaces a !! procedure pointer. subroutine s_convert_to_mixture_variables(q_vf, i, j, k, rho, gamma, pi_inf, qv, Re_K, G_K, G) @@ -1181,31 +1185,6 @@ contains end subroutine s_finalize_variables_conversion_module - !> Total energy per unit volume, thermodynamic terms only. Callers add magnetic and elastic energy, which are not - !! equation-of-state terms. The chemistry and relativistic branches use a different relation and stay open-coded. - subroutine s_compute_energy(pres, alpha_rho_K, alpha_K, vel_sum, E) - - $:GPU_ROUTINE(function_name='s_compute_energy', parallelism='[seq]', cray_inline=True) - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: alpha_rho_K, alpha_K - #:else - real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K - #:endif - real(wp), intent(in) :: pres, vel_sum - real(wp), intent(out) :: E - real(wp) :: rho, gamma, pi_inf, qv - - call s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho, gamma, pi_inf, qv) - - ! E = dyn_p + (1 - alf)(gamma p + pi_inf) + qv. Only the liquid's internal energy is diluted; - ! qv is already an energy density. - E = gamma*pres + pi_inf - if (bubbles_euler) E = E*(1._wp - alpha_K(num_fluids)) - E = E + qv + 5.e-1_wp*rho*vel_sum - - end subroutine s_compute_energy - !> 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) diff --git a/src/simulation/include/riemann_state_inline.fpp b/src/simulation/include/riemann_state_inline.fpp new file mode 100644 index 000000000..73dd6034c --- /dev/null +++ b/src/simulation/include/riemann_state_inline.fpp @@ -0,0 +1,23 @@ +!> +!! @file riemann_state_inline.fpp +!! @brief Rename / private helpers for modules that compile a private copy of riemann_state_inline_procs.fpp. +!! Usage as in common/include/inline_eos.fpp. + +#:include 'riemann_state_inline_procs.fpp' + +#:set RS_INLINE_NAMES = ['f_low_Mach_zcoef', 'f_low_Mach_pcorr_hllc', 's_apply_low_Mach_velocity', 's_compute_average_state', & + & 's_compute_hypoelastic_interface_energy', 'f_compute_hllc_star_momentum_flux'] + +#:def RS_INLINE_USE(module) + use ${module}$, & + #:for n in RS_INLINE_NAMES[:-1] + & ${n}$_ext => ${n}$, & + #:endfor + & ${RS_INLINE_NAMES[-1]}$_ext => ${RS_INLINE_NAMES[-1]}$ +#:enddef + +#:def RS_INLINE_PRIVATE() + #:for n in RS_INLINE_NAMES + private :: ${n}$ + #:endfor +#:enddef diff --git a/src/simulation/include/riemann_state_inline_procs.fpp b/src/simulation/include/riemann_state_inline_procs.fpp new file mode 100644 index 000000000..56ebc3560 --- /dev/null +++ b/src/simulation/include/riemann_state_inline_procs.fpp @@ -0,0 +1,191 @@ +!> +!! @file riemann_state_inline_procs.fpp +!! @brief Per-face Riemann helpers, one source for m_riemann_state and the Riemann solver modules that call them from +!! GPU kernels, kept in the Fypp macro RS_INLINE_PROCS. Same reason and usage as common/include/inline_eos_procs.fpp; +!! the rename / private helpers are RS_INLINE_USE / RS_INLINE_PRIVATE in riemann_state_inline.fpp. Edit the routines +!! here only. + +#:def RS_INLINE_PROCS() + !> Low-Mach parameter of Thornber et al. JCP (2008): the larger of the two face Mach numbers, capped at one so the correction + !! switches itself off once the flow is no longer low speed. + function f_low_Mach_zcoef(vel_L_rms, vel_R_rms, c_L, c_R) result(zcoef) + + $:GPU_ROUTINE(function_name='f_low_Mach_zcoef', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: vel_L_rms, vel_R_rms !< Left and right squared velocity magnitudes + real(wp), intent(in) :: c_L, c_R !< Left and right sound speeds + real(wp) :: zcoef + + zcoef = min(1._wp, max(vel_L_rms**5.e-1_wp/c_L, vel_R_rms**5.e-1_wp/c_R)) + + end function f_low_Mach_zcoef + + !> The same correction for the HLLC flux, where the star state supplies the pressure jump directly and the correction scales + !! with the mass flux through the acoustic waves instead. Zero unless low_Mach == 1. + function f_low_Mach_pcorr_hllc(vel_L_rms, vel_R_rms, c_L, c_R, rho_L, rho_R, s_L, s_R, vel_L_norm, vel_R_norm) result(pcorr) + + $:GPU_ROUTINE(function_name='f_low_Mach_pcorr_hllc', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: vel_L_rms, vel_R_rms !< Left and right squared velocity magnitudes + real(wp), intent(in) :: c_L, c_R !< Left and right sound speeds + real(wp), intent(in) :: rho_L, rho_R !< Left and right densities + real(wp), intent(in) :: s_L, s_R !< Left and right wave speeds + real(wp), intent(in) :: vel_L_norm, vel_R_norm !< Left and right wave-normal velocities + real(wp) :: pcorr + + pcorr = 0._wp + if (low_Mach == 1) then + pcorr = rho_L*rho_R*(s_L - vel_L_norm)*(s_R - vel_R_norm)*(vel_R_norm - vel_L_norm)/(rho_R*(s_R - vel_R_norm) & + & - rho_L*(s_L - vel_L_norm))*(f_low_Mach_zcoef(vel_L_rms, vel_R_rms, c_L, c_R) - 1._wp) + end if + + end function f_low_Mach_pcorr_hllc + + !> The alternative low-Mach treatment of Thornber et al. JCP (2008) selected by low_Mach == 2: rather than correct the flux, + !! blend the wave-normal velocities towards their mean before the wave speeds are computed, which is why this mutates its + !! arguments and must be called ahead of s_L, s_R and s_S. The tangential velocities and vel_L/R_rms are deliberately left + !! untouched. + subroutine s_apply_low_Mach_velocity(vel_L_rms, vel_R_rms, c_L, c_R, vel_L_norm, vel_R_norm) + + $:GPU_ROUTINE(function_name='s_apply_low_Mach_velocity', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: vel_L_rms, vel_R_rms !< Left and right squared velocity magnitudes + real(wp), intent(in) :: c_L, c_R !< Left and right sound speeds + real(wp), intent(inout) :: vel_L_norm, vel_R_norm !< Left and right wave-normal velocities, blended in place + real(wp) :: zcoef, vel_L_tmp, vel_R_tmp + + zcoef = f_low_Mach_zcoef(vel_L_rms, vel_R_rms, c_L, c_R) + + vel_L_tmp = 5.e-1_wp*((vel_L_norm + vel_R_norm) + zcoef*(vel_L_norm - vel_R_norm)) + vel_R_tmp = 5.e-1_wp*((vel_L_norm + vel_R_norm) + zcoef*(vel_R_norm - vel_L_norm)) + + vel_L_norm = vel_L_tmp + vel_R_norm = vel_R_tmp + + end subroutine s_apply_low_Mach_velocity + + !> Interface-averaged state that the pressure-based wave-speed estimate reads. avg_state selects between the density-weighted + !! Roe average, which costs eight square roots per face, and the plain arithmetic mean; unlike the other solver switches this + !! one is not implied by the call site, so the dispatch stays here. + subroutine s_compute_average_state(rho_L, rho_R, vel_L, vel_R, H_L, H_R, gamma_L, gamma_R, qv_L, qv_R, rho_avg, vel_avg_rms, & + & H_avg, gamma_avg, qv_avg) + + $:GPU_ROUTINE(function_name='s_compute_average_state', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: rho_L, rho_R !< Left and right densities + real(wp), intent(in) :: H_L, H_R !< Left and right total enthalpies + real(wp), intent(in) :: gamma_L, gamma_R !< Left and right specific heat ratio functions + real(wp), intent(in) :: qv_L, qv_R !< Left and right reference energies + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in) :: vel_L, vel_R + #:else + real(wp), dimension(num_vels), intent(in) :: vel_L, vel_R + #:endif + real(wp), intent(out) :: rho_avg, H_avg, gamma_avg, qv_avg + real(wp), intent(out) :: vel_avg_rms !< Squared magnitude of the averaged velocity, summed over all components + integer :: i + + vel_avg_rms = 0._wp + + if (avg_state == avg_state_roe) then + rho_avg = sqrt(rho_L*rho_R) + + $:GPU_LOOP(parallelism='[seq]') + do i = 1, num_vels + vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))**2._wp/(sqrt(rho_L) + sqrt(rho_R))**2._wp + end do + + H_avg = (sqrt(rho_L)*H_L + sqrt(rho_R)*H_R)/(sqrt(rho_L) + sqrt(rho_R)) + gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/(sqrt(rho_L) + sqrt(rho_R)) + qv_avg = (sqrt(rho_L)*qv_L + sqrt(rho_R)*qv_R)/(sqrt(rho_L) + sqrt(rho_R)) + else + rho_avg = 5.e-1_wp*(rho_L + rho_R) + + $:GPU_LOOP(parallelism='[seq]') + do i = 1, num_vels + vel_avg_rms = vel_avg_rms + (5.e-1_wp*(vel_L(i) + vel_R(i)))**2._wp + end do + + H_avg = 5.e-1_wp*(H_L + H_R) + gamma_avg = 5.e-1_wp*(gamma_L + gamma_R) + qv_avg = 5.e-1_wp*(qv_L + qv_R) + end if + + end subroutine s_compute_average_state + + !> Accumulate the hypoelastic stress contribution to the energies of the left and right Riemann states: mix the shear modulus + !! over the fluids, add the elastic energy of each stress component (doubled for the shear components) on each side whose + !! mixture modulus is non-negligible, then scale the returned moduli by the continuum damage state when damage is modeled + !! (energy uses the undamaged modulus; the damaged moduli feed the callers' wave speeds). The elastic shear stresses are loaded + !! from the state buffers by the caller, which reuses them for the stress fluxes and elastic wave speeds. The G > verysmall + !! per-side gate is a deliberate maintainer ruling that replaces HLL's former hard-coded G > 1000 stability floor, retiring its + !! "TODO take out if statement if stable without". + subroutine s_compute_hypoelastic_interface_energy(nf, alpha_L, alpha_R, damage_L, damage_R, tau_e_L, tau_e_R, G_L, G_R, E_L, & + & E_R) + + $:GPU_ROUTINE(function_name='s_compute_hypoelastic_interface_energy', parallelism='[seq]', cray_inline=True) + + integer, intent(in) :: nf !< Number of fluids to mix the shear modulus over + real(wp), dimension(nf), intent(in) :: alpha_L, alpha_R !< Left and right volume fractions + real(wp), intent(in) :: damage_L, damage_R !< Continuum damage states (referenced only when cont_damage) + real(wp), dimension(6), intent(in) :: tau_e_L, tau_e_R !< Left and right elastic shear stresses + real(wp), intent(out) :: G_L, G_R !< Left and right mixture shear moduli + real(wp), intent(inout) :: E_L, E_R !< Left and right state energies + integer :: i !< Loop iterator + + G_L = 0._wp; G_R = 0._wp + + $:GPU_LOOP(parallelism='[seq]') + do i = 1, nf + G_L = G_L + alpha_L(i)*Gs_rs(i) + G_R = G_R + alpha_R(i)*Gs_rs(i) + end do + + ! Elastic energy uses the undamaged modulus, so this loop precedes the damage scaling + $:GPU_LOOP(parallelism='[seq]') + do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1 + ! Elastic contribution to energy if G large enough + if (G_L > verysmall) then + E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) + ! Double for shear stresses + if (any(eqn_idx%stress%beg - 1 + i == shear_indices)) then + E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) + end if + end if + if (G_R > verysmall) then + E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) + ! Double for shear stresses + if (any(eqn_idx%stress%beg - 1 + i == shear_indices)) then + E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) + end if + end if + end do + + ! Damage scaling applies only to the returned moduli (used for wave speeds) + if (cont_damage) then + G_L = G_L*max((1._wp - damage_L), 0._wp) + G_R = G_R*max((1._wp - damage_R), 0._wp) + end if + + end subroutine s_compute_hypoelastic_interface_energy + + !> Compute the advective part of the HLLC star-state momentum flux in the wave-normal direction (pressure excluded), used to + !! assemble the geometrical source flux of the cylindrical and azimuthal sweeps. + function f_compute_hllc_star_momentum_flux(rho_L, rho_R, vel_L_norm, vel_R_norm, s_M, s_P, s_S, xi_L, xi_R, xi_M, xi_P, & + & dir_flg_norm) result(flux_mom) + + $:GPU_ROUTINE(function_name='f_compute_hllc_star_momentum_flux', parallelism='[seq]', cray_inline=True) + + real(wp), intent(in) :: rho_L, rho_R !< Left and right densities + real(wp), intent(in) :: vel_L_norm, vel_R_norm !< Left and right wave-normal velocities + real(wp), intent(in) :: s_M, s_P, s_S !< Clamped left/right and contact wave speeds + real(wp), intent(in) :: xi_L, xi_R, xi_M, xi_P !< Star-state compression factors and upwind selectors + real(wp), intent(in) :: dir_flg_norm !< Direction flag of the wave-normal direction + real(wp) :: flux_mom + + flux_mom = xi_M*(rho_L*(vel_L_norm*vel_L_norm + s_M*(xi_L*(dir_flg_norm*s_S + (1._wp - dir_flg_norm)*vel_L_norm) & + & - vel_L_norm))) + xi_P*(rho_R*(vel_R_norm*vel_R_norm + s_P*(xi_R*(dir_flg_norm*s_S + (1._wp & + & - dir_flg_norm)*vel_R_norm) - vel_R_norm))) + + end function f_compute_hllc_star_momentum_flux +#:enddef diff --git a/src/simulation/m_acoustic_src.fpp b/src/simulation/m_acoustic_src.fpp index 5ca556ede..9f9fcee73 100644 --- a/src/simulation/m_acoustic_src.fpp +++ b/src/simulation/m_acoustic_src.fpp @@ -3,6 +3,7 @@ !! @brief Contains module m_acoustic_src #:include 'macros.fpp' +#:include 'inline_eos.fpp' !> @brief One-way acoustic source injection, Maeda and Colonius JCP (2017) module m_acoustic_src @@ -11,12 +12,14 @@ 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 + $:EOS_INLINE_USE('m_eos') + use m_helper_basic implicit none + $:EOS_INLINE_PRIVATE() + private; public :: s_initialize_acoustic_src, s_precalculate_acoustic_spatial_sources, s_acoustic_src_calculations integer, allocatable, dimension(:) :: pulse, support @@ -60,6 +63,8 @@ module m_acoustic_src contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Initialize the acoustic source module impure subroutine s_initialize_acoustic_src diff --git a/src/simulation/m_bubbles_EE.fpp b/src/simulation/m_bubbles_EE.fpp index c1f954742..2c2fe6644 100644 --- a/src/simulation/m_bubbles_EE.fpp +++ b/src/simulation/m_bubbles_EE.fpp @@ -3,6 +3,7 @@ !! @brief Contains module @ref m_bubbles_ee "m_bubbles_EE" #:include 'macros.fpp' +#:include 'inline_eos.fpp' !> @brief Computes ensemble-averaged (Euler--Euler) bubble source terms for radius, velocity, pressure, and mass transfer module m_bubbles_EE @@ -11,11 +12,13 @@ module m_bubbles_EE use m_global_parameters use m_mpi_proxy use m_variables_conversion - use m_eos + $:EOS_INLINE_USE('m_eos') use m_bubbles implicit none + $:EOS_INLINE_PRIVATE() + real(wp), allocatable, dimension(:,:,:) :: bub_adv_src real(wp), allocatable, dimension(:,:,:,:) :: bub_r_src, bub_v_src, bub_p_src, bub_m_src $:GPU_DECLARE(create='[bub_adv_src, bub_r_src, bub_v_src, bub_p_src, bub_m_src]') @@ -28,6 +31,8 @@ module m_bubbles_EE contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Initialize the Euler-Euler bubble module impure subroutine s_initialize_bubbles_EE_module diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index 0f40e6fe6..31c440644 100644 --- a/src/simulation/m_bubbles_EL.fpp +++ b/src/simulation/m_bubbles_EL.fpp @@ -3,6 +3,7 @@ !! @brief Contains module @ref m_bubbles_el "m_bubbles_EL" #:include 'macros.fpp' +#:include 'inline_eos.fpp' !> @brief Tracks Lagrangian bubbles and couples their dynamics to the Eulerian flow via volume averaging module m_bubbles_EL @@ -12,7 +13,7 @@ module m_bubbles_EL use m_bubbles_EL_kernels use m_bubbles use m_variables_conversion - use m_eos + $:EOS_INLINE_USE('m_eos') use m_compile_specific use m_boundary_common use m_helper_basic @@ -25,6 +26,8 @@ module m_bubbles_EL implicit none + $:EOS_INLINE_PRIVATE() + ! (nBub) integer, allocatable, dimension(:,:) :: lag_id !< Global and local IDs real(wp), allocatable, dimension(:) :: bub_R0 !< Initial bubble radius @@ -80,6 +83,8 @@ module m_bubbles_EL contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Initializes the lagrangian subgrid bubble solver impure subroutine s_initialize_bubbles_EL_module(q_cons_vf, bc_type) diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index 49cc12c70..8ebd4a168 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -6,6 +6,7 @@ !! boundaries #:include 'case.fpp' #:include 'macros.fpp' +#:include 'inline_eos.fpp' module m_cbc @@ -19,8 +20,12 @@ module m_cbc & get_mole_fractions use m_thermochem_state, only: get_mixture_caloric_state + $:EOS_INLINE_USE('m_eos') + implicit none + $:EOS_INLINE_PRIVATE() + private; public :: s_initialize_cbc_module, s_cbc, s_finalize_cbc_module ! The cell-average primitive variables. They are obtained by reshaping (RS) q_prim_vf in the coordinate direction normal to the @@ -85,6 +90,8 @@ module m_cbc contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Initialize the CBC module impure subroutine s_initialize_cbc_module diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 112ea52a2..d72693d66 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -3,6 +3,7 @@ !! @brief Contains module m_data_output #:include 'macros.fpp' +#:include 'inline_eos.fpp' #:include 'case.fpp' !> @brief Writes solution data, run-time stability diagnostics (ICFL, VCFL, CCFL, Rc), and probe/center-of-mass files @@ -21,8 +22,12 @@ module m_data_output use m_boundary_common use m_constants, only: model_eqns_5eq, precision_single + $:EOS_INLINE_USE('m_eos') + implicit none + $:EOS_INLINE_PRIVATE() + private public :: s_initialize_data_output_module, s_open_run_time_information_file, s_open_probe_files, & & s_write_run_time_information, s_write_data_files, s_write_serial_data_files, s_write_parallel_data_files, & @@ -54,6 +59,8 @@ module m_data_output contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Write data files. Dispatch subroutine that replaces procedure pointer. impure subroutine s_write_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta) diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp index 9da6817e8..a07f12b5b 100644 --- a/src/simulation/m_hypoelastic.fpp +++ b/src/simulation/m_hypoelastic.fpp @@ -3,6 +3,7 @@ !! @brief Contains module m_hypoelastic #:include 'macros.fpp' +#:include 'inline_eos.fpp' !> @brief Computes hypoelastic stress-rate source terms and damage-state evolution module m_hypoelastic @@ -11,10 +12,12 @@ module m_hypoelastic use m_global_parameters use m_finite_differences use m_helper - use m_eos, only: s_phase_bulk_modulus + $:EOS_INLINE_USE('m_eos', only=[]) implicit none + $:EOS_INLINE_PRIVATE() + private; public :: s_initialize_hypoelastic_module, s_finalize_hypoelastic_module, & & s_compute_hypoelastic_rhs_finite_diff_per_sweep, s_compute_hypoelastic_rhs_iface, & & s_compute_hypoelastic_rhs_axisym_geom_iface, s_compute_hypoelastic_rhs_axisym_geom_dual_pass, s_compute_damage_state, & @@ -38,6 +41,8 @@ module m_hypoelastic contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Initialize the hypoelastic module impure subroutine s_initialize_hypoelastic_module diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 614d0dc4b..3f8ff6617 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -3,6 +3,7 @@ !! @brief Contains module m_ibm #:include 'macros.fpp' +#:include 'inline_eos.fpp' !> @brief Ghost-node immersed boundary method: locates ghost/image points, computes interpolation coefficients, and corrects the !! flow state @@ -12,10 +13,10 @@ module m_ibm use m_global_parameters use m_mpi_proxy use m_variables_conversion - use m_eos + use m_constants + $:EOS_INLINE_USE('m_eos') use m_helper use m_helper_basic - use m_constants use m_compute_levelset use m_ib_patches use m_viscous @@ -29,6 +30,8 @@ module m_ibm implicit none + $:EOS_INLINE_PRIVATE() + private :: s_compute_image_points, s_compute_interpolation_coeffs, s_interpolate_image_point, s_find_ghost_points, & & s_find_num_ghost_points, s_compute_ghost_point_pressure, s_compute_ghost_point_velocity ; public :: s_initialize_ibm_module, s_ibm_setup, s_ibm_correct_state, s_finalize_ibm_module, s_report_ibm_surface @@ -70,6 +73,8 @@ module m_ibm contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Allocates memory for the variables in the IBM module impure subroutine s_initialize_ibm_module() diff --git a/src/simulation/m_igr.fpp b/src/simulation/m_igr.fpp index 846a039a9..5e27f216d 100644 --- a/src/simulation/m_igr.fpp +++ b/src/simulation/m_igr.fpp @@ -4,6 +4,7 @@ #:include 'case.fpp' #:include 'macros.fpp' +#:include 'inline_eos.fpp' !> @brief Iterative ghost rasterization (IGR) for sharp immersed boundary treatment module m_igr @@ -11,13 +12,15 @@ module m_igr use m_derived_types use m_global_parameters use m_variables_conversion - use m_eos + $:EOS_INLINE_USE('m_eos') use m_mpi_proxy use m_helper use m_boundary_common implicit none + $:EOS_INLINE_PRIVATE() + private; public :: s_initialize_igr_module, s_igr_iterative_solve, s_igr_riemann_solver, s_igr_sigma_x, s_igr_flux_add, & & s_finalize_igr_module @@ -90,6 +93,8 @@ module m_igr contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Initialize the IGR module subroutine s_initialize_igr_module() diff --git a/src/simulation/m_pressure_relaxation.fpp b/src/simulation/m_pressure_relaxation.fpp index 6c6dfe6dd..3447b79f9 100644 --- a/src/simulation/m_pressure_relaxation.fpp +++ b/src/simulation/m_pressure_relaxation.fpp @@ -4,6 +4,7 @@ #:include 'case.fpp' #:include 'macros.fpp' +#:include 'inline_eos.fpp' !> @brief Pressure relaxation for the six-equation multi-component model via Newton--Raphson equilibration and volume-fraction !! correction @@ -13,10 +14,12 @@ module m_pressure_relaxation use m_global_parameters use m_mpi_proxy, only: s_mpi_abort 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 + $:EOS_INLINE_USE('m_eos', only=['s_phase_density_on_isentrope']) implicit none + $:EOS_INLINE_PRIVATE() + private; public :: s_pressure_relaxation_procedure, s_report_pressure_relaxation !> Cell updates in which the Newton below stopped on its iteration cap rather than on the tolerance. Counted because the answer @@ -26,6 +29,8 @@ module m_pressure_relaxation contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> The main pressure relaxation procedure subroutine s_pressure_relaxation_procedure(q_cons_vf) diff --git a/src/simulation/m_qbmm.fpp b/src/simulation/m_qbmm.fpp index cde6830a7..a3914ef71 100644 --- a/src/simulation/m_qbmm.fpp +++ b/src/simulation/m_qbmm.fpp @@ -4,6 +4,7 @@ #:include 'case.fpp' #:include 'macros.fpp' +#:include 'inline_eos.fpp' !> @brief Quadrature-based moment methods (QBMM) for polydisperse bubble moment inversion and transport module m_qbmm @@ -12,13 +13,15 @@ module m_qbmm use m_global_parameters use m_mpi_proxy use m_variables_conversion - use m_eos + $:EOS_INLINE_USE('m_eos') use m_helper_basic use m_helper use m_constants, only: bubble_model_keller_miksis, bubble_model_rayleigh_plesset implicit none + $:EOS_INLINE_PRIVATE() + private; public :: s_initialize_qbmm_module, s_mom_inv, s_coeff, s_compute_qbmm_rhs real(wp), allocatable, dimension(:,:,:,:,:) :: momrhs @@ -39,6 +42,8 @@ module m_qbmm contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Initialize the QBMM module impure subroutine s_initialize_qbmm_module diff --git a/src/simulation/m_reactive_burn.fpp b/src/simulation/m_reactive_burn.fpp index 73da758cb..ee17f9967 100644 --- a/src/simulation/m_reactive_burn.fpp +++ b/src/simulation/m_reactive_burn.fpp @@ -3,6 +3,7 @@ !! @brief Contains module m_reactive_burn #:include 'macros.fpp' +#:include 'inline_eos.fpp' #:include 'case.fpp' !> @brief Condensed-phase reactive burn: a pressure-driven programmed-burn source that converts a "reactant" fluid into a "product" @@ -15,15 +16,18 @@ module m_reactive_burn use m_global_parameters - use m_variables_conversion, only: s_compute_mixture_coefficients - use m_eos, only: s_phase_temperature, f_pressure + $:EOS_INLINE_USE('m_eos', only=[]) implicit none + $:EOS_INLINE_PRIVATE() + private; public :: s_compute_reactive_burn, s_reactive_burn_substep contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Programmed-burn rate dlambda/dt for one cell state. Both the RHS source and the operator-split integrator call this, so the !! rate law is stated once. !! @param pres Mixture pressure diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index 7e74b33ed..eb02c4b6a 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -4,6 +4,7 @@ #:include 'case.fpp' #:include 'macros.fpp' +#:include 'inline_eos.fpp' !> @brief Assemble the right-hand side of the governing equations using finite-volume flux differencing, Riemann solvers, and !! physical source terms @@ -13,7 +14,7 @@ module m_rhs use m_global_parameters use m_mpi_proxy use m_variables_conversion - use m_eos + $:EOS_INLINE_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 @@ -41,6 +42,8 @@ module m_rhs implicit none + $:EOS_INLINE_PRIVATE() + private; public :: s_initialize_rhs_module, s_compute_rhs, s_finalize_rhs_module type(vector_field) :: q_cons_qp !< WENO-reconstructed cell-average conservative variables at quadrature points @@ -137,6 +140,8 @@ module m_rhs contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Initialize the RHS module impure subroutine s_initialize_rhs_module diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index 5cddef70c..752b46b91 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -5,13 +5,14 @@ !> @brief HLL approximate Riemann solver, Harten et al. SIAM Review (1983) #:include 'case.fpp' #:include 'macros.fpp' +#:include 'inline_eos.fpp' module m_riemann_solver_hll use m_derived_types use m_global_parameters use m_variables_conversion - use m_eos + $:EOS_INLINE_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 @@ -21,8 +22,12 @@ module m_riemann_solver_hll implicit none + $:EOS_INLINE_PRIVATE() + contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> HLL approximate Riemann solver, Harten et al. SIAM Review (1983) subroutine s_hll_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, & & dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, & diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index 7bb7eec1a..fe33be594 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -5,13 +5,15 @@ !> @brief HLLC Riemann solver with contact restoration, Toro et al. Shock Waves (1994) #:include 'case.fpp' #:include 'macros.fpp' +#:include 'inline_eos.fpp' +#:include 'riemann_state_inline.fpp' module m_riemann_solver_hllc use m_derived_types use m_global_parameters use m_variables_conversion - use m_eos + $:EOS_INLINE_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 @@ -20,12 +22,18 @@ module m_riemann_solver_hllc use m_chemistry 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 + $:RS_INLINE_USE('m_riemann_state') implicit none + $:EOS_INLINE_PRIVATE() + $:RS_INLINE_PRIVATE() + contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() + $:RS_INLINE_PROCS() !> HLLC Riemann solver with contact restoration, Toro et al. Shock Waves (1994) subroutine s_hllc_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, & & dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, & diff --git a/src/simulation/m_riemann_solver_hlld.fpp b/src/simulation/m_riemann_solver_hlld.fpp index bce95254d..9f201bc34 100644 --- a/src/simulation/m_riemann_solver_hlld.fpp +++ b/src/simulation/m_riemann_solver_hlld.fpp @@ -5,6 +5,7 @@ !> @brief HLLD approximate Riemann solver for MHD, Miyoshi & Kusano JCP (2005) #:include 'case.fpp' #:include 'macros.fpp' +#:include 'inline_eos.fpp' module m_riemann_solver_hlld @@ -13,10 +14,16 @@ module m_riemann_solver_hlld use m_variables_conversion use m_riemann_state + $:EOS_INLINE_USE('m_eos') + implicit none + $:EOS_INLINE_PRIVATE() + contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> HLLD Riemann solver for MHD, Miyoshi & Kusano JCP (2005) subroutine s_hlld_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, & & dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, & diff --git a/src/simulation/m_riemann_solver_hypo_hlld.fpp b/src/simulation/m_riemann_solver_hypo_hlld.fpp index 2fb798e6a..a72463015 100644 --- a/src/simulation/m_riemann_solver_hypo_hlld.fpp +++ b/src/simulation/m_riemann_solver_hypo_hlld.fpp @@ -5,23 +5,28 @@ !> @brief Dual-pass HLLD approximate Riemann solver for hypoelastic flows, with non-conservative interface-velocity coupling #:include 'case.fpp' #:include 'macros.fpp' +#:include 'inline_eos.fpp' module m_riemann_solver_hypo_hlld use m_derived_types use m_global_parameters use m_variables_conversion - use m_eos + $:EOS_INLINE_USE('m_eos') use m_riemann_state implicit none + $:EOS_INLINE_PRIVATE() + !> Four-epsilon roundoff band for G_eff: two mixture products, their sum, and the final signed sum. Relative to the operand !! scale; not a physical stiffness cutoff. real(wp), parameter :: G_eff_roundoff_factor = 4._wp contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Classify the position of the interface (xi = 0) in the five-wave HLLD fan: 0 left of S_L, 1..4 for the four inner wedges, 5 !! right of S_R. Single definition of the wave-fan geometry shared by the flux fold, the NC face-velocity export, and the !! axisymmetric face-state pick, which must stay consistent. Valid fans only: a fallback face exports the HLL/one-sided traces diff --git a/src/simulation/m_riemann_solver_lf.fpp b/src/simulation/m_riemann_solver_lf.fpp index 6893a3168..d69df05b4 100644 --- a/src/simulation/m_riemann_solver_lf.fpp +++ b/src/simulation/m_riemann_solver_lf.fpp @@ -5,6 +5,7 @@ !> @brief Lax-Friedrichs (Rusanov) approximate Riemann solver #:include 'case.fpp' #:include 'macros.fpp' +#:include 'inline_eos.fpp' module m_riemann_solver_lf @@ -16,10 +17,16 @@ module m_riemann_solver_lf use m_thermochem_state, only: get_mixture_caloric_state use m_riemann_state + $:EOS_INLINE_USE('m_eos') + implicit none + $:EOS_INLINE_PRIVATE() + contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Lax-Friedrichs (Rusanov) approximate Riemann solver subroutine s_lf_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, & & dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, flux_src_vf, & diff --git a/src/simulation/m_riemann_state.fpp b/src/simulation/m_riemann_state.fpp index aeb34c8a9..137aa9452 100644 --- a/src/simulation/m_riemann_state.fpp +++ b/src/simulation/m_riemann_state.fpp @@ -6,6 +6,7 @@ !! helpers #:include 'case.fpp' #:include 'macros.fpp' +#:include 'riemann_state_inline_procs.fpp' module m_riemann_state @@ -86,20 +87,10 @@ contains end function f_elastic_signal_speed - !> Low-Mach parameter of Thornber et al. JCP (2008): the larger of the two face Mach numbers, capped at one so the correction - !! switches itself off once the flow is no longer low speed. - function f_low_Mach_zcoef(vel_L_rms, vel_R_rms, c_L, c_R) result(zcoef) - - $:GPU_ROUTINE(function_name='f_low_Mach_zcoef', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: vel_L_rms, vel_R_rms !< Left and right squared velocity magnitudes - real(wp), intent(in) :: c_L, c_R !< Left and right sound speeds - real(wp) :: zcoef - - zcoef = min(1._wp, max(vel_L_rms**5.e-1_wp/c_L, vel_R_rms**5.e-1_wp/c_R)) - - end function f_low_Mach_zcoef - + ! Per-face Riemann helpers: the source is the Fypp macro RS_INLINE_PROCS (include/riemann_state_inline_procs.fpp), so the + ! Riemann solver modules can compile a module-local copy that inlines into their kernels. This is the public + ! definition. + $:RS_INLINE_PROCS() !> Low-Mach pressure correction added to the HLL and Lax-Friedrichs fluxes, which restores the pressure jump that the !! dissipation of those fluxes over-damps at low Mach number. Zero unless low_Mach == 1. function f_low_Mach_pcorr_hll(vel_L_rms, vel_R_rms, c_L, c_R, rho_L, rho_R, s_M, s_P) result(pcorr) @@ -119,99 +110,6 @@ contains end function f_low_Mach_pcorr_hll - !> The same correction for the HLLC flux, where the star state supplies the pressure jump directly and the correction scales - !! with the mass flux through the acoustic waves instead. Zero unless low_Mach == 1. - function f_low_Mach_pcorr_hllc(vel_L_rms, vel_R_rms, c_L, c_R, rho_L, rho_R, s_L, s_R, vel_L_norm, vel_R_norm) result(pcorr) - - $:GPU_ROUTINE(function_name='f_low_Mach_pcorr_hllc', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: vel_L_rms, vel_R_rms !< Left and right squared velocity magnitudes - real(wp), intent(in) :: c_L, c_R !< Left and right sound speeds - real(wp), intent(in) :: rho_L, rho_R !< Left and right densities - real(wp), intent(in) :: s_L, s_R !< Left and right wave speeds - real(wp), intent(in) :: vel_L_norm, vel_R_norm !< Left and right wave-normal velocities - real(wp) :: pcorr - - pcorr = 0._wp - if (low_Mach == 1) then - pcorr = rho_L*rho_R*(s_L - vel_L_norm)*(s_R - vel_R_norm)*(vel_R_norm - vel_L_norm)/(rho_R*(s_R - vel_R_norm) & - & - rho_L*(s_L - vel_L_norm))*(f_low_Mach_zcoef(vel_L_rms, vel_R_rms, c_L, c_R) - 1._wp) - end if - - end function f_low_Mach_pcorr_hllc - - !> The alternative low-Mach treatment of Thornber et al. JCP (2008) selected by low_Mach == 2: rather than correct the flux, - !! blend the wave-normal velocities towards their mean before the wave speeds are computed, which is why this mutates its - !! arguments and must be called ahead of s_L, s_R and s_S. The tangential velocities and vel_L/R_rms are deliberately left - !! untouched. - subroutine s_apply_low_Mach_velocity(vel_L_rms, vel_R_rms, c_L, c_R, vel_L_norm, vel_R_norm) - - $:GPU_ROUTINE(function_name='s_apply_low_Mach_velocity', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: vel_L_rms, vel_R_rms !< Left and right squared velocity magnitudes - real(wp), intent(in) :: c_L, c_R !< Left and right sound speeds - real(wp), intent(inout) :: vel_L_norm, vel_R_norm !< Left and right wave-normal velocities, blended in place - real(wp) :: zcoef, vel_L_tmp, vel_R_tmp - - zcoef = f_low_Mach_zcoef(vel_L_rms, vel_R_rms, c_L, c_R) - - vel_L_tmp = 5.e-1_wp*((vel_L_norm + vel_R_norm) + zcoef*(vel_L_norm - vel_R_norm)) - vel_R_tmp = 5.e-1_wp*((vel_L_norm + vel_R_norm) + zcoef*(vel_R_norm - vel_L_norm)) - - vel_L_norm = vel_L_tmp - vel_R_norm = vel_R_tmp - - end subroutine s_apply_low_Mach_velocity - - !> Interface-averaged state that the pressure-based wave-speed estimate reads. avg_state selects between the density-weighted - !! Roe average, which costs eight square roots per face, and the plain arithmetic mean; unlike the other solver switches this - !! one is not implied by the call site, so the dispatch stays here. - subroutine s_compute_average_state(rho_L, rho_R, vel_L, vel_R, H_L, H_R, gamma_L, gamma_R, qv_L, qv_R, rho_avg, vel_avg_rms, & - & H_avg, gamma_avg, qv_avg) - - $:GPU_ROUTINE(function_name='s_compute_average_state', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: rho_L, rho_R !< Left and right densities - real(wp), intent(in) :: H_L, H_R !< Left and right total enthalpies - real(wp), intent(in) :: gamma_L, gamma_R !< Left and right specific heat ratio functions - real(wp), intent(in) :: qv_L, qv_R !< Left and right reference energies - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: vel_L, vel_R - #:else - real(wp), dimension(num_vels), intent(in) :: vel_L, vel_R - #:endif - real(wp), intent(out) :: rho_avg, H_avg, gamma_avg, qv_avg - real(wp), intent(out) :: vel_avg_rms !< Squared magnitude of the averaged velocity, summed over all components - integer :: i - - vel_avg_rms = 0._wp - - if (avg_state == avg_state_roe) then - rho_avg = sqrt(rho_L*rho_R) - - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_vels - vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))**2._wp/(sqrt(rho_L) + sqrt(rho_R))**2._wp - end do - - H_avg = (sqrt(rho_L)*H_L + sqrt(rho_R)*H_R)/(sqrt(rho_L) + sqrt(rho_R)) - gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/(sqrt(rho_L) + sqrt(rho_R)) - qv_avg = (sqrt(rho_L)*qv_L + sqrt(rho_R)*qv_R)/(sqrt(rho_L) + sqrt(rho_R)) - else - rho_avg = 5.e-1_wp*(rho_L + rho_R) - - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_vels - vel_avg_rms = vel_avg_rms + (5.e-1_wp*(vel_L(i) + vel_R(i)))**2._wp - end do - - H_avg = 5.e-1_wp*(H_L + H_R) - gamma_avg = 5.e-1_wp*(gamma_L + gamma_R) - qv_avg = 5.e-1_wp*(qv_L + qv_R) - end if - - end subroutine s_compute_average_state - !> Roe-averaged reacting-mixture quantities: replaces gamma_avg with the mixture Cp/Cv and builds the c_sum_Yi_Phi term !! s_compute_speed_of_sound_avg needs. vel_avg_rms must be the full squared magnitude - its Phi_avg and vel_sum terms cancel to !! leave the Roe sound speed, and only do so for the full magnitude. @@ -1245,82 +1143,6 @@ contains end subroutine s_compute_interface_reynolds - !> Accumulate the hypoelastic stress contribution to the energies of the left and right Riemann states: mix the shear modulus - !! over the fluids, add the elastic energy of each stress component (doubled for the shear components) on each side whose - !! mixture modulus is non-negligible, then scale the returned moduli by the continuum damage state when damage is modeled - !! (energy uses the undamaged modulus; the damaged moduli feed the callers' wave speeds). The elastic shear stresses are loaded - !! from the state buffers by the caller, which reuses them for the stress fluxes and elastic wave speeds. The G > verysmall - !! per-side gate is a deliberate maintainer ruling that replaces HLL's former hard-coded G > 1000 stability floor, retiring its - !! "TODO take out if statement if stable without". - subroutine s_compute_hypoelastic_interface_energy(nf, alpha_L, alpha_R, damage_L, damage_R, tau_e_L, tau_e_R, G_L, G_R, E_L, & - & E_R) - - $:GPU_ROUTINE(function_name='s_compute_hypoelastic_interface_energy', parallelism='[seq]', cray_inline=True) - - integer, intent(in) :: nf !< Number of fluids to mix the shear modulus over - real(wp), dimension(nf), intent(in) :: alpha_L, alpha_R !< Left and right volume fractions - real(wp), intent(in) :: damage_L, damage_R !< Continuum damage states (referenced only when cont_damage) - real(wp), dimension(6), intent(in) :: tau_e_L, tau_e_R !< Left and right elastic shear stresses - real(wp), intent(out) :: G_L, G_R !< Left and right mixture shear moduli - real(wp), intent(inout) :: E_L, E_R !< Left and right state energies - integer :: i !< Loop iterator - - G_L = 0._wp; G_R = 0._wp - - $:GPU_LOOP(parallelism='[seq]') - do i = 1, nf - G_L = G_L + alpha_L(i)*Gs_rs(i) - G_R = G_R + alpha_R(i)*Gs_rs(i) - end do - - ! Elastic energy uses the undamaged modulus, so this loop precedes the damage scaling - $:GPU_LOOP(parallelism='[seq]') - do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1 - ! Elastic contribution to energy if G large enough - if (G_L > verysmall) then - E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) - ! Double for shear stresses - if (any(eqn_idx%stress%beg - 1 + i == shear_indices)) then - E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) - end if - end if - if (G_R > verysmall) then - E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) - ! Double for shear stresses - if (any(eqn_idx%stress%beg - 1 + i == shear_indices)) then - E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) - end if - end if - end do - - ! Damage scaling applies only to the returned moduli (used for wave speeds) - if (cont_damage) then - G_L = G_L*max((1._wp - damage_L), 0._wp) - G_R = G_R*max((1._wp - damage_R), 0._wp) - end if - - end subroutine s_compute_hypoelastic_interface_energy - - !> Compute the advective part of the HLLC star-state momentum flux in the wave-normal direction (pressure excluded), used to - !! assemble the geometrical source flux of the cylindrical and azimuthal sweeps. - function f_compute_hllc_star_momentum_flux(rho_L, rho_R, vel_L_norm, vel_R_norm, s_M, s_P, s_S, xi_L, xi_R, xi_M, xi_P, & - & dir_flg_norm) result(flux_mom) - - $:GPU_ROUTINE(function_name='f_compute_hllc_star_momentum_flux', parallelism='[seq]', cray_inline=True) - - real(wp), intent(in) :: rho_L, rho_R !< Left and right densities - real(wp), intent(in) :: vel_L_norm, vel_R_norm !< Left and right wave-normal velocities - real(wp), intent(in) :: s_M, s_P, s_S !< Clamped left/right and contact wave speeds - real(wp), intent(in) :: xi_L, xi_R, xi_M, xi_P !< Star-state compression factors and upwind selectors - real(wp), intent(in) :: dir_flg_norm !< Direction flag of the wave-normal direction - real(wp) :: flux_mom - - flux_mom = xi_M*(rho_L*(vel_L_norm*vel_L_norm + s_M*(xi_L*(dir_flg_norm*s_S + (1._wp - dir_flg_norm)*vel_L_norm) & - & - vel_L_norm))) + xi_P*(rho_R*(vel_R_norm*vel_R_norm + s_P*(xi_R*(dir_flg_norm*s_S + (1._wp & - & - dir_flg_norm)*vel_R_norm) - vel_R_norm))) - - end function f_compute_hllc_star_momentum_flux - !> Reshape and copy the Riemann-solver flux buffers back to the physical-space output arrays for the selected sweep direction, !! finalizing the Riemann solve. Two variants are emitted from one template so the shared unpermute logic cannot drift apart: !! the plain routine also copies the advection flux_src set and the grid_geometry==3 z-sweep geometric source flux, while the diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index a918fec7c..1d1e6f5aa 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -3,6 +3,7 @@ !! @brief Contains module m_time_steppers #:include 'macros.fpp' +#:include 'inline_eos.fpp' #:include 'case.fpp' !> @brief Total-variation-diminishing (TVD) Runge--Kutta time integrators (1st-, 2nd-, and 3rd-order SSP) @@ -32,8 +33,12 @@ module m_time_steppers use m_derived_variables use m_constants, only: model_eqns_6eq, time_stepper_rk1, time_stepper_rk2, time_stepper_rk3 + $:EOS_INLINE_USE('m_eos') + implicit none + $:EOS_INLINE_PRIVATE() + type(vector_field), allocatable, dimension(:) :: q_cons_ts !< Cell-average conservative variables at each time-stage (TS) type(scalar_field), allocatable, dimension(:) :: q_prim_vf !< Cell-average primitive variables at the current time-stage type(scalar_field), allocatable, dimension(:) :: rhs_vf !< Cell-average RHS variables at the current time-stage @@ -63,6 +68,8 @@ module m_time_steppers contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Initialize the time steppers module impure subroutine s_initialize_time_steppers_module diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp index 37afd8988..28883abd6 100644 --- a/src/simulation/m_viscous.fpp +++ b/src/simulation/m_viscous.fpp @@ -3,6 +3,7 @@ !! @brief Contains module m_viscous #:include 'case.fpp' #:include 'macros.fpp' +#:include 'inline_eos.fpp' #! Fill Re_visc(1:2) at a cylindrical-axis cell. Non-Newtonian path uses a #! Herschel-Bulkley mixture inv-Re from the local strain rate; Newtonian path is @@ -18,10 +19,12 @@ module m_viscous use m_muscl use m_helper use m_finite_differences - use m_variables_conversion, only: s_compute_mixture_coefficients + $:EOS_INLINE_USE('m_eos') use m_constants, only: model_eqns_5eq, recon_type_weno, recon_type_muscl use m_hb_function + $:EOS_INLINE_PRIVATE() + private; public s_get_viscous, s_compute_viscous_stress_cylindrical_boundary, s_initialize_viscous_module, & & s_reconstruct_cell_boundary_values_visc_deriv, s_finalize_viscous_module, s_compute_viscous_stress_tensor @@ -34,6 +37,8 @@ module m_viscous contains + ! Module-local copies of the hot EOS routines, so they inline into this module's kernels + $:EOS_INLINE_PROCS() !> Initialize the viscous module impure subroutine s_initialize_viscous_module