diff --git a/docs/module_categories.json b/docs/module_categories.json index 1391908cad..bc8550454a 100644 --- a/docs/module_categories.json +++ b/docs/module_categories.json @@ -13,6 +13,7 @@ "m_riemann_solver_hllc", "m_muscl", "m_variables_conversion", + "m_thermodynamics", "m_thinc" ] }, diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index 13e84c8b79..5aaae06f21 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/src/common/m_global_parameters_common.fpp @@ -123,6 +123,15 @@ module m_global_parameters_common integer :: mpi_info_int !> @} + !> @name Per-fluid property tables. In the simulation target these live in m_global_parameters; for pre and post process they + !! live here so the thermodynamics interface can read them from global parameters without a module cycle. + !> @{ +#ifndef MFC_SIMULATION + real(wp), allocatable, public, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps + $:GPU_DECLARE(create='[gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps]') +#endif + !> @} + contains !> Initialize equation-index state (eqn_idx, sys_size, b_size, tensor_size) from the namelist parameters. This is the shared diff --git a/src/common/m_thermodynamics.fpp b/src/common/m_thermodynamics.fpp new file mode 100644 index 0000000000..e032b750c5 --- /dev/null +++ b/src/common/m_thermodynamics.fpp @@ -0,0 +1,248 @@ +!> +!! @file m_thermodynamics.fpp +!! @brief Contains module m_thermodynamics + +#:include 'macros.fpp' +#:include 'case.fpp' + +!> @brief Central equation-of-state interface. Every genuine thermodynamic query (pressure, internal energy, temperature, sound +!! speed, and the rho c^2 derivative) is answered through one module so no solver path carries its own EOS algebra. Compile-time +!! chemistry selects the adapter: the stiffened_gas adapter (analytic mixture EOS) or the pyrometheus adapter +!! (temperature-then-query against the thermochemistry surface). This is a leaf module (it uses global parameters and the +!! thermochemistry surface, never m_variables_conversion) so no cycle forms. +module m_thermodynamics + + use m_derived_types + use m_global_parameters + use m_constants, only: model_eqns_4eq, model_eqns_5eq, model_eqns_6eq, avg_state_roe, verysmall, sgm_eps + use m_thermochem, only: num_species, get_temperature, get_pressure, gas_constant, get_mixture_molecular_weight, & + & get_mixture_energy_mass + + implicit none + + private + public :: s_compute_pressure + public :: s_compute_internal_energy + public :: f_compute_temperature + public :: f_validate_state +#ifndef MFC_PRE_PROCESS + public :: s_compute_speed_of_sound + public :: f_bulk_modulus +#endif + +contains + + !> Pressure from the conserved energy and mixture properties. The #:if chemistry split selects the stiffened-gas adapter + !! (analytic inverse of the mixture EOS) or the Pyrometheus adapter (solve temperature, then query pressure). Arithmetic is + !! relocated unchanged from the conversion module. + subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, T, stress, mom, G, pres_mag) + + $:GPU_ROUTINE(function_name='s_compute_pressure',parallelism='[seq]', cray_noinline=True) + + real(stp), intent(in) :: energy, alf + real(wp), intent(in) :: dyn_p + real(wp), intent(in) :: pi_inf, gamma, rho, qv + real(wp), intent(out) :: pres + real(wp), intent(inout) :: T + real(stp), intent(in), optional :: stress, mom + real(wp), intent(in), optional :: G, pres_mag + + ! Chemistry + real(wp), dimension(1:num_species), intent(in) :: rhoYks + real(wp), dimension(1:num_species) :: Y_rs + real(wp) :: E_e + real(wp) :: e_Per_Kg, Pdyn_Per_Kg + real(wp) :: T_guess + integer :: s !< Generic loop iterator + #:if not chemistry + ! Depending on model_eqns and bubbles_euler, the appropriate procedure for computing pressure is targeted by the + ! procedure pointer + + if (mhd) then + ! MHD pressure: subtract magnetic pressure from total energy + pres = (energy - dyn_p - pi_inf - qv - pres_mag)/gamma + else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then + ! Gamma/pi_inf model or five-equation model (Allaire et al. JCP 2002): p from mixture EOS + pres = (energy - dyn_p - pi_inf - qv)/gamma + else if ((model_eqns /= model_eqns_4eq) .and. bubbles_euler) then + ! Bubble-augmented pressure with void fraction correction + pres = ((energy - dyn_p)/(1._wp - alf) - pi_inf - qv)/gamma + else + ! Four-equation model (Kapila et al. PoF 2001): Tait EOS inversion + pres = (pref + pi_inf)*(energy/(rhoref*(1 - alf)))**(1/gamma + 1) - pi_inf + end if + + if (hypoelasticity .and. present(G)) then + ! Subtract elastic strain energy before computing pressure (hypoelastic model) + E_e = 0._wp + do s = eqn_idx%stress%beg, eqn_idx%stress%end + if (G > 0) then + E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G) + ! Double for shear stresses + if (any(s == shear_indices)) then + E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G) + end if + end if + end do + + pres = (energy - 0.5_wp*(mom**2._wp)/rho - pi_inf - qv - E_e)/gamma + end if + #:else + ! Reacting mixture pressure from temperature and species + Y_rs(:) = rhoYks(:)/rho + e_Per_Kg = energy/rho + Pdyn_Per_Kg = dyn_p/rho + + T_guess = T + + call get_temperature(e_Per_Kg - Pdyn_Per_Kg, T_guess, Y_rs, .true., T) + call get_pressure(rho, T, Y_rs, pres) + #:endif + + end subroutine s_compute_pressure + + !> Conserved energy from the primitive pressure and mixture properties: the inverse of s_compute_pressure over the same model + !! branches. The chemistry adapter builds the mixture internal energy from temperature and species. Arithmetic is relocated + !! unchanged from the conversion module. + subroutine s_compute_internal_energy(pres, alf, dyn_p, pi_inf, gamma, rho, qv, Ys, energy, pres_mag) + + $:GPU_ROUTINE(function_name='s_compute_internal_energy',parallelism='[seq]', cray_noinline=True) + + real(stp), intent(in) :: pres, alf + real(wp), intent(in) :: dyn_p, pi_inf, gamma, rho, qv + real(stp), intent(out) :: energy + real(wp), intent(in), optional :: pres_mag + real(wp), dimension(1:num_species), intent(in) :: Ys + real(wp) :: e_mix, mix_mol_weight, T + #:if not chemistry + ! Computing the energy from the pressure + if (mhd) then + ! MHD energy includes magnetic pressure contribution + energy = gamma*pres + dyn_p + pres_mag + pi_inf + qv + else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then + ! Five-equation model (Allaire et al. JCP 2002): E = Gamma*p + 0.5*rho*|u|^2 + pi_inf + qv + energy = gamma*pres + dyn_p + pi_inf + qv + else if ((model_eqns /= model_eqns_4eq) .and. bubbles_euler) then + ! Bubble-augmented energy with void fraction correction + energy = dyn_p + (1._wp - alf)*(gamma*pres + pi_inf) + else + ! Four-equation model (Kapila et al. PoF 2001): Tait EOS, no conserved energy variable + energy = 0._wp + end if + #:else + ! Reacting mixture: compute conserved energy from species mass fractions and temperature + call get_mixture_molecular_weight(Ys, mix_mol_weight) + T = pres*mix_mol_weight/(gas_constant*rho) + call get_mixture_energy_mass(T, Ys, e_mix) + energy = dyn_p + rho*e_mix + #:endif + + end subroutine s_compute_internal_energy + + !> Mixture temperature from pressure via the ideal-gas relation T = p/(rho R_gas), where R_gas is the specific gas constant of + !! the local composition. The one temperature query shared by the conversion routines, so the chemistry cons->prim path no + !! longer carries its own copy. Arithmetic is relocated unchanged from the conversion module. + pure function f_compute_temperature(pres, rho, R_gas) result(T) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: pres, rho, R_gas + real(wp) :: T + + T = pres/rho/R_gas + + end function f_compute_temperature + + !> Thermodynamic-state admissibility: a state is physical only when its squared sound speed is non-negative, i.e. the wave speed + !! is real. The mixture_err sound-speed guards in the interface and in post-processing both consult this one predicate instead + !! of testing the sign inline, so the validity criterion lives in one place. Pure and device-safe: it returns a flag and never + !! aborts, because GPU device code cannot. Later adapters (Mie-Gruneisen, JWL) tighten the admissible domain here. + pure function f_validate_state(sound_speed_sq) result(is_admissible) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: sound_speed_sq + logical :: is_admissible + + is_admissible = (sound_speed_sq >= 0._wp) + + end function f_validate_state + +#ifndef MFC_PRE_PROCESS + !> Stiffened-gas bulk modulus, the thermodynamic derivative $\rho c^2 = ((\gamma+1)p+\pi_\infty)/\gamma$. Shared by the + !! alt_soundspeed Wood mixture sound speed here, the five-equation K-div-u source in m_rhs, and the post-process derived sound + !! speed, so the one derivative lives in one place. + pure function f_bulk_modulus(gamma, pres, pi_inf) result(blkmod) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: gamma, pres, pi_inf + real(wp) :: blkmod + + blkmod = ((gamma + 1._wp)*pres + pi_inf)/gamma + + end function f_bulk_modulus + + !> Speed of sound from thermodynamic state, supporting the several equation-of-state models. Arithmetic is relocated unchanged + !! from the conversion module: the non-chemistry branches hold c squared until a single final sqrt, the chemistry branch returns + !! c already square-rooted, and the mixture_err floor is preserved. + subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c, qv) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: pres + real(wp), intent(in) :: rho, gamma, pi_inf, qv + real(wp), intent(in) :: H + #: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(in) :: vel_sum + real(wp), intent(in) :: c_c + real(wp), intent(out) :: c + real(wp) :: blkmod1, blkmod2 + integer :: q + + 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 + c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho) + end if + else if (relativity) then ! Relativistic sound speed + c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H) + else + if (alt_soundspeed) then ! Wood's mixture sound speed via bulk moduli + blkmod1 = f_bulk_modulus(gammas(1), pres, pi_infs(1)) + blkmod2 = f_bulk_modulus(gammas(2), pres, pi_infs(2)) + c = (1._wp/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2))) + else if (model_eqns == model_eqns_6eq) then ! Six-equation model sound speed + c = 0._wp + $:GPU_LOOP(parallelism='[seq]') + do q = 1, num_fluids + c = c + adv(q)*gs_min(q)*(pres + pi_infs(q)/(gammas(q) + 1._wp)) + end do + c = c/rho + else if (((model_eqns == model_eqns_4eq) .or. (model_eqns == model_eqns_5eq .and. bubbles_euler))) then + ! Sound speed for bubble mixture to order O(\alpha) + + if (mpp_lim .and. (num_fluids > 1)) then + c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/rho + else + c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/(rho*(1._wp - adv(num_fluids))) + end if + else + c = (H - 5.e-1*vel_sum - qv/rho)/gamma + end if + + if (mixture_err .and. .not. f_validate_state(c)) then + c = 100._wp*sgm_eps + else + c = sqrt(c) + end if + end if + + end subroutine s_compute_speed_of_sound +#endif +end module m_thermodynamics diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index fe17b89c27..4f15e2b067 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -14,10 +14,15 @@ module m_variables_conversion use m_helper_basic use m_helper use m_constants, only: riemann_solver_hll, riemann_solver_hlld, model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, & - & model_eqns_4eq, avg_state_roe + & model_eqns_4eq use m_thermochem, only: num_species, get_temperature, get_pressure, gas_constant, get_mixture_molecular_weight, & & get_mixture_energy_mass + use m_thermodynamics, only: s_compute_pressure, s_compute_internal_energy, f_compute_temperature +#ifndef MFC_PRE_PROCESS + use m_thermodynamics, only: s_compute_speed_of_sound +#endif + implicit none private @@ -39,12 +44,6 @@ module m_variables_conversion #endif s_finalize_variables_conversion_module - ! In simulation, gammas, pi_infs, and qvs are already declared in m_global_variables -#ifndef MFC_SIMULATION - real(wp), allocatable, public, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps - $:GPU_DECLARE(create='[gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps]') -#endif - real(wp), allocatable, dimension(:) :: Gs_vc integer, allocatable, dimension(:) :: bubrs_vc real(wp), allocatable, dimension(:,:) :: Res_vc @@ -79,73 +78,6 @@ contains end subroutine s_convert_to_mixture_variables - !> Compute the pressure from the appropriate equation of state - subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, T, stress, mom, G, pres_mag) - - $:GPU_ROUTINE(function_name='s_compute_pressure',parallelism='[seq]', cray_noinline=True) - - real(stp), intent(in) :: energy, alf - real(wp), intent(in) :: dyn_p - real(wp), intent(in) :: pi_inf, gamma, rho, qv - real(wp), intent(out) :: pres - real(wp), intent(inout) :: T - real(stp), intent(in), optional :: stress, mom - real(wp), intent(in), optional :: G, pres_mag - - ! Chemistry - real(wp), dimension(1:num_species), intent(in) :: rhoYks - real(wp), dimension(1:num_species) :: Y_rs - real(wp) :: E_e - real(wp) :: e_Per_Kg, Pdyn_Per_Kg - real(wp) :: T_guess - integer :: s !< Generic loop iterator - #:if not chemistry - ! Depending on model_eqns and bubbles_euler, the appropriate procedure for computing pressure is targeted by the - ! procedure pointer - - if (mhd) then - ! MHD pressure: subtract magnetic pressure from total energy - pres = (energy - dyn_p - pi_inf - qv - pres_mag)/gamma - else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then - ! Gamma/pi_inf model or five-equation model (Allaire et al. JCP 2002): p from mixture EOS - pres = (energy - dyn_p - pi_inf - qv)/gamma - else if ((model_eqns /= model_eqns_4eq) .and. bubbles_euler) then - ! Bubble-augmented pressure with void fraction correction - pres = ((energy - dyn_p)/(1._wp - alf) - pi_inf - qv)/gamma - else - ! Four-equation model (Kapila et al. PoF 2001): Tait EOS inversion - pres = (pref + pi_inf)*(energy/(rhoref*(1 - alf)))**(1/gamma + 1) - pi_inf - end if - - if (hypoelasticity .and. present(G)) then - ! Subtract elastic strain energy before computing pressure (hypoelastic model) - E_e = 0._wp - do s = eqn_idx%stress%beg, eqn_idx%stress%end - if (G > 0) then - E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G) - ! Double for shear stresses - if (any(s == shear_indices)) then - E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G) - end if - end if - end do - - pres = (energy - 0.5_wp*(mom**2._wp)/rho - pi_inf - qv - E_e)/gamma - end if - #:else - ! Reacting mixture pressure from temperature and species - Y_rs(:) = rhoYks(:)/rho - e_Per_Kg = energy/rho - Pdyn_Per_Kg = dyn_p/rho - - T_guess = T - - call get_temperature(e_Per_Kg - Pdyn_Per_Kg, T_guess, Y_rs, .true., T) - call get_pressure(rho, T, Y_rs, pres) - #:endif - - end subroutine s_compute_pressure - !> Convert mixture variables to density, gamma, pi_inf, and qv for the gamma/pi_inf model. Given conservative or primitive !! variables, transfers the density, specific heat ratio function and the liquid stiffness function from q_vf to rho, gamma and !! pi_inf. @@ -799,7 +731,6 @@ contains real(wp), dimension(2) :: Re_K integer :: i, j, k, l !< Generic loop iterators real(wp), dimension(num_species) :: Ys - real(wp) :: e_mix, mix_mol_weight, T real(wp) :: pres_mag real(wp) :: Ga !< Lorentz factor (gamma in relativity) real(wp) :: h !< relativistic enthalpy @@ -897,43 +828,29 @@ contains end do if (chemistry) then - ! Reacting mixture: compute conserved energy from species mass fractions and temperature do i = eqn_idx%species%beg, eqn_idx%species%end Ys(i - eqn_idx%species%beg + 1) = q_prim_vf(i)%sf(j, k, l) q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l) end do + end if - call get_mixture_molecular_weight(Ys, mix_mol_weight) - T = q_prim_vf(eqn_idx%E)%sf(j, k, l)*mix_mol_weight/(gas_constant*rho) - call get_mixture_energy_mass(T, Ys, e_mix) - - q_cons_vf(eqn_idx%E)%sf(j, k, l) = dyn_pres + rho*e_mix - else - ! Computing the energy from the pressure - if (mhd) then - if (n == 0) then - pres_mag = 0.5_wp*(Bx0**2 + q_prim_vf(eqn_idx%B%beg)%sf(j, k, & - & l)**2 + q_prim_vf(eqn_idx%B%beg + 1)%sf(j, k, l)**2) - else - pres_mag = 0.5_wp*(q_prim_vf(eqn_idx%B%beg)%sf(j, k, l)**2 + q_prim_vf(eqn_idx%B%beg + 1)%sf(j, & - & k, l)**2 + q_prim_vf(eqn_idx%B%beg + 2)%sf(j, k, l)**2) - end if - ! MHD energy includes magnetic pressure contribution - q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, & - & l) + dyn_pres + pres_mag + pi_inf + qv - else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then - ! Five-equation model (Allaire et al. JCP 2002): E = Gamma*p + 0.5*rho*|u|^2 + pi_inf + qv - q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + dyn_pres + pi_inf + qv - else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler)) then - ! Bubble-augmented energy with void fraction correction - q_cons_vf(eqn_idx%E)%sf(j, k, l) = dyn_pres + (1._wp - q_prim_vf(eqn_idx%alf)%sf(j, k, & - & l))*(gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + pi_inf) + if (mhd) then + if (n == 0) then + pres_mag = 0.5_wp*(Bx0**2 + q_prim_vf(eqn_idx%B%beg)%sf(j, k, & + & l)**2 + q_prim_vf(eqn_idx%B%beg + 1)%sf(j, k, l)**2) else - ! Four-equation model (Kapila et al. PoF 2001): Tait EOS, no conserved energy variable - q_cons_vf(eqn_idx%E)%sf(j, k, l) = 0._wp + pres_mag = 0.5_wp*(q_prim_vf(eqn_idx%B%beg)%sf(j, k, l)**2 + q_prim_vf(eqn_idx%B%beg + 1)%sf(j, k, & + & l)**2 + q_prim_vf(eqn_idx%B%beg + 2)%sf(j, k, l)**2) end if + else + pres_mag = 0._wp end if + ! eqn_idx%alf = 0 for non-bubble configs; clamp the index (value read only in the bubble branch) + call s_compute_internal_energy(q_prim_vf(eqn_idx%E)%sf(j, k, l), q_prim_vf(max(1, eqn_idx%alf))%sf(j, k, l), & + & dyn_pres, pi_inf, gamma, rho, qv, Ys, q_cons_vf(eqn_idx%E)%sf(j, k, l), & + & pres_mag=pres_mag) + ! Six-equation model (Saurel et al. JCP 2009): compute per-phase internal energies if (model_eqns == model_eqns_6eq) then do i = 1, num_fluids @@ -1115,7 +1032,7 @@ contains ! Computing the energy from the internal energy of the mixture call get_mixture_molecular_weight(Y_k, mix_mol_weight) R_gas = gas_constant/mix_mol_weight - T_K = pres_K/rho_K/R_gas + T_K = f_compute_temperature(pres_K, rho_K, R_gas) call get_mixture_energy_mass(T_K, Y_K, E_K) E_K = rho_K*E_K + 5.e-1_wp*rho_K*vel_K_sum else @@ -1249,68 +1166,6 @@ contains end subroutine s_finalize_variables_conversion_module -#ifndef MFC_PRE_PROCESS - !> Compute the speed of sound from thermodynamic state variables, supporting multiple equation-of-state models. - subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c, qv) - - $:GPU_ROUTINE(parallelism='[seq]') - - real(wp), intent(in) :: pres - real(wp), intent(in) :: rho, gamma, pi_inf, qv - real(wp), intent(in) :: H - #: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(in) :: vel_sum - real(wp), intent(in) :: c_c - real(wp), intent(out) :: c - real(wp) :: blkmod1, blkmod2 - integer :: q - - 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 - c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho) - end if - else if (relativity) then ! Relativistic sound speed - c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H) - else - if (alt_soundspeed) then ! Wood's mixture sound speed via bulk moduli - blkmod1 = ((gammas(1) + 1._wp)*pres + pi_infs(1))/gammas(1) - blkmod2 = ((gammas(2) + 1._wp)*pres + pi_infs(2))/gammas(2) - c = (1._wp/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2))) - else if (model_eqns == model_eqns_6eq) then ! Six-equation model sound speed - c = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do q = 1, num_fluids - c = c + adv(q)*gs_min(q)*(pres + pi_infs(q)/(gammas(q) + 1._wp)) - end do - c = c/rho - else if (((model_eqns == model_eqns_4eq) .or. (model_eqns == model_eqns_5eq .and. bubbles_euler))) then - ! Sound speed for bubble mixture to order O(\alpha) - - if (mpp_lim .and. (num_fluids > 1)) then - c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/rho - else - c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/(rho*(1._wp - adv(num_fluids))) - end if - else - c = (H - 5.e-1*vel_sum - qv/rho)/gamma - end if - - if (mixture_err .and. c < 0._wp) then - c = 100._wp*sgm_eps - else - c = sqrt(c) - end if - end if - - end subroutine s_compute_speed_of_sound -#endif - #ifndef MFC_PRE_PROCESS !> Compute the fast magnetosonic wave speed from the sound speed, density, and magnetic field components. subroutine s_compute_fast_magnetosonic_speed(rho, c, B, norm, c_fast, h) diff --git a/src/post_process/m_derived_variables.fpp b/src/post_process/m_derived_variables.fpp index 05f1fc1805..82e91b99c9 100644 --- a/src/post_process/m_derived_variables.fpp +++ b/src/post_process/m_derived_variables.fpp @@ -11,6 +11,7 @@ module m_derived_variables use m_mpi_proxy use m_helper_basic use m_variables_conversion + use m_thermodynamics, only: f_bulk_modulus, f_validate_state use m_constants, only: model_eqns_gamma_law implicit none @@ -102,16 +103,16 @@ contains do j = -offset_y%beg, n + offset_y%end do i = -offset_x%beg, m + offset_x%end if (alt_soundspeed .neqv. .true.) then - q_sf(i, j, k) = (((gamma_sf(i, j, k) + 1._wp)*q_prim_vf(eqn_idx%E)%sf(i, j, k) + pi_inf_sf(i, j, & - & k))/(gamma_sf(i, j, k)*rho_sf(i, j, k))) + q_sf(i, j, k) = f_bulk_modulus(gamma_sf(i, j, k), q_prim_vf(eqn_idx%E)%sf(i, j, k), pi_inf_sf(i, j, & + & k))/rho_sf(i, j, k) else - blkmod1 = ((gammas(1) + 1._wp)*q_prim_vf(eqn_idx%E)%sf(i, j, k) + pi_infs(1))/gammas(1) - blkmod2 = ((gammas(2) + 1._wp)*q_prim_vf(eqn_idx%E)%sf(i, j, k) + pi_infs(2))/gammas(2) + blkmod1 = f_bulk_modulus(gammas(1), q_prim_vf(eqn_idx%E)%sf(i, j, k), pi_infs(1)) + blkmod2 = f_bulk_modulus(gammas(2), q_prim_vf(eqn_idx%E)%sf(i, j, k), pi_infs(2)) q_sf(i, j, k) = (1._wp/(rho_sf(i, j, k)*(q_prim_vf(eqn_idx%adv%beg)%sf(i, j, & & k)/blkmod1 + (1._wp - q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k))/blkmod2))) end if - if (mixture_err .and. q_sf(i, j, k) < 0._wp) then + if (mixture_err .and. .not. f_validate_state(q_sf(i, j, k))) then q_sf(i, j, k) = 1.e-16_wp else q_sf(i, j, k) = sqrt(q_sf(i, j, k)) diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index 8176f7b838..0f316c7315 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -13,6 +13,7 @@ module m_rhs use m_global_parameters use m_mpi_proxy use m_variables_conversion + use m_thermodynamics, only: f_bulk_modulus 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 @@ -870,10 +871,10 @@ contains do q_loop = 0, p do l_loop = 0, n do k_loop = 0, m - blkmod1(k_loop, l_loop, q_loop) = ((gammas(1) + 1._wp)*q_prim_vf%vf(eqn_idx%E)%sf(k_loop, l_loop, & - & q_loop) + pi_infs(1))/gammas(1) - blkmod2(k_loop, l_loop, q_loop) = ((gammas(2) + 1._wp)*q_prim_vf%vf(eqn_idx%E)%sf(k_loop, l_loop, & - & q_loop) + pi_infs(2))/gammas(2) + blkmod1(k_loop, l_loop, q_loop) = f_bulk_modulus(gammas(1), q_prim_vf%vf(eqn_idx%E)%sf(k_loop, l_loop, & + & q_loop), pi_infs(1)) + blkmod2(k_loop, l_loop, q_loop) = f_bulk_modulus(gammas(2), q_prim_vf%vf(eqn_idx%E)%sf(k_loop, l_loop, & + & q_loop), pi_infs(2)) alpha1(k_loop, l_loop, q_loop) = q_cons_vf%vf(eqn_idx%adv%beg)%sf(k_loop, l_loop, q_loop) if (bubbles_euler) then