diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 872c04b39..bb01c8a89 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -121,6 +121,29 @@ module m_derived_types type(int_bounds_info) :: x, y, z end type bc_xyz_info + !> Thermodynamic state handed to the equation-of-state routines. + !> + !> Carries scalars only. The volume-fraction array `adv` stays a separate argument: a derived-type + !> component cannot have a runtime extent, and `num_fluids` is a parameter only under case + !> optimization, so `dimension(num_fluids)` will not compile in a general build. Padding to + !> num_fluids_max would compile but put ten reals into a per-cell private struct on device. + !> + !> `H` is the specific total enthalpy and must include `qv`, because the sound-speed relation + !> subtracts `qv/rho`. Build states with s_eos_state so that invariant holds by construction + !> rather than by convention; s_eos_state_roe exists for the Roe-averaged paths, which supply an + !> `H` that is deliberately not the exact state enthalpy. + !> Contains no allocatable members - safe to use inside device routines. + type eos_state + real(wp) :: rho !< Mixture density + real(wp) :: pres !< Pressure + real(wp) :: gamma !< Stiffened-gas gamma (1/(Gamma-1)) + real(wp) :: pi_inf !< Stiffened-gas stiffness + real(wp) :: qv !< Heat of formation (volumetric) + real(wp) :: vel_sum !< |u|^2 + real(wp) :: H !< Specific total enthalpy, including qv + real(wp) :: c_c !< Roe-averaged chemistry sound-speed term (0 when unused) + end type eos_state + !> QBMM moment index mappings - separate from bub beg/end so eqn_idx contains no allocatables. type qbmm_idx_info integer, dimension(:), allocatable :: rs !< R moment indices per bubble bin diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index e56eb78df..0cf680b4c 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -25,7 +25,7 @@ module m_variables_conversion & s_convert_mixture_to_mixture_variables, s_convert_species_to_mixture_variables, & & s_convert_species_to_mixture_variables_kernel, s_convert_conservative_to_primitive_variables, & & s_convert_primitive_to_conservative_variables, s_convert_primitive_to_flux_variables, s_compute_pressure, & - & s_compute_species_fraction, s_compute_speed_of_sound, s_compute_fast_magnetosonic_speed, & + & s_compute_species_fraction, s_compute_speed_of_sound, s_compute_fast_magnetosonic_speed, s_eos_state, s_eos_state_roe, & & s_finalize_variables_conversion_module, gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps real(wp), allocatable, dimension(:) :: Gs_vc @@ -1202,55 +1202,92 @@ contains end subroutine s_finalize_variables_conversion_module + !> Build an exact thermodynamic state. The specific total enthalpy is derived here rather than supplied, so it cannot disagree + !! with qv: H = (E + p)/rho with E = gamma*p + pi_inf + qv + rho*|u|^2/2. Every caller that previously open-coded the closed + !! form should use this, which makes the defect in #1707 unrepresentable. + subroutine s_eos_state(s, pres, rho, gamma, pi_inf, qv, vel_sum) + + $:GPU_ROUTINE(function_name='s_eos_state', parallelism='[seq]', cray_inline=True) + + type(eos_state), intent(out) :: s + real(wp), intent(in) :: pres, rho, gamma, pi_inf, qv, vel_sum + + s%pres = pres + s%rho = rho + s%gamma = gamma + s%pi_inf = pi_inf + s%qv = qv + s%vel_sum = vel_sum + s%c_c = 0._wp + s%H = ((gamma + 1._wp)*pres + pi_inf + qv)/rho + 5.e-1_wp*vel_sum + + end subroutine s_eos_state + + !> Build a state whose enthalpy is supplied by the caller. The Roe-averaged Riemann paths, the chemistry Roe branch and the + !! relativistic branch all pass an H that is deliberately not the exact enthalpy of the state, so they cannot use s_eos_state. + subroutine s_eos_state_roe(s, pres, rho, gamma, pi_inf, qv, vel_sum, H) + + $:GPU_ROUTINE(function_name='s_eos_state_roe', parallelism='[seq]', cray_inline=True) + + type(eos_state), intent(out) :: s + real(wp), intent(in) :: pres, rho, gamma, pi_inf, qv, vel_sum, H + + s%pres = pres + s%rho = rho + s%gamma = gamma + s%pi_inf = pi_inf + s%qv = qv + s%vel_sum = vel_sum + s%H = H + s%c_c = 0._wp + + end subroutine s_eos_state_roe + !> 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) + subroutine s_compute_speed_of_sound(s, adv, c) $:GPU_ROUTINE(parallelism='[seq]') - real(wp), intent(in) :: pres - real(wp), intent(in) :: rho, gamma, pi_inf, qv - real(wp), intent(in) :: H + type(eos_state), intent(in) :: s #: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)) + if (avg_state == avg_state_roe .and. abs(s%c_c) > verysmall) then + c = sqrt(s%c_c - (s%gamma - 1.0_wp)*(s%vel_sum - s%H)) else - c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho) + c = sqrt((1.0_wp + 1.0_wp/s%gamma)*s%pres/s%rho) end if else if (relativity) then ! Relativistic sound speed - c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H) + c = sqrt((1._wp + 1._wp/s%gamma)*s%pres/s%rho/s%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))) + blkmod1 = ((gammas(1) + 1._wp)*s%pres + pi_infs(1))/gammas(1) + blkmod2 = ((gammas(2) + 1._wp)*s%pres + pi_infs(2))/gammas(2) + c = (1._wp/(s%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)) + c = c + adv(q)*gs_min(q)*(s%pres + pi_infs(q)/(gammas(q) + 1._wp)) end do - c = c/rho + c = c/s%rho else if (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 + c = (1._wp/s%gamma + 1._wp)*(s%pres + s%pi_inf/(s%gamma + 1._wp))/s%rho else - c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/(rho*(1._wp - adv(num_fluids))) + c = (1._wp/s%gamma + 1._wp)*(s%pres + s%pi_inf/(s%gamma + 1._wp))/(s%rho*(1._wp - adv(num_fluids))) end if else - c = (H - 5.e-1*vel_sum - qv/rho)/gamma + c = (s%H - 5.e-1*s%vel_sum - s%qv/s%rho)/s%gamma end if if (mixture_err .and. c < 0._wp) then diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index 1d5ddfedb..0a36631ad 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -1243,6 +1243,7 @@ contains real(wp), dimension(num_vels) :: vel real(wp), dimension(num_fluids) :: adv integer :: i, j, k, l, s !< looping indices + type(eos_state) :: eos_s Egk = 0._wp Elp = 0._wp @@ -1288,7 +1289,8 @@ contains H = ((gamma + 1._wp)*pres + pi_inf + qv)/rho - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, 0._wp, 0._wp, c, qv) + call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp) + call s_compute_speed_of_sound(eos_s, adv, c) Ma = maxvel/c if (Ma > MaxMa .and. (adv(1) > (1.0_wp - 1.0e-10_wp))) then diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 9c697c0ea..0c7fa0859 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -187,10 +187,11 @@ contains & -offset_z%beg:p + offset_z%end) :: liutex_mag real(wp), dimension(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end,-offset_z%beg:p + offset_z%end, & & 3) :: liutex_axis - integer :: i, j, k, l, kx, ky, kz, kf, j_glb, k_glb, l_glb - character(50) :: filename - logical :: file_exists - integer :: x_beg, x_end, y_beg, y_end, z_beg, z_end + integer :: i, j, k, l, kx, ky, kz, kf, j_glb, k_glb, l_glb + character(50) :: filename + logical :: file_exists + integer :: x_beg, x_end, y_beg, y_end, z_beg, z_end + type(eos_state) :: eos_s if (output_partial_domain) then call s_define_output_region @@ -533,8 +534,9 @@ contains H = ((gamma_sf(i, j, k) + 1._wp)*pres + pi_inf_sf(i, j, k) + qv_sf(i, j, k))/rho_sf(i, j, k) - call s_compute_speed_of_sound(pres, rho_sf(i, j, k), gamma_sf(i, j, k), pi_inf_sf(i, j, k), H, adv, & - & 0._wp, 0._wp, c, qv_sf(i, j, k)) + call s_eos_state_roe(eos_s, pres, rho_sf(i, j, k), gamma_sf(i, j, k), pi_inf_sf(i, j, k), qv_sf(i, j, k), & + & 0._wp, H) + call s_compute_speed_of_sound(eos_s, adv, c) out%q_sf(i, j, k) = c end do diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index bc72c87bd..87b927ebb 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -511,6 +511,7 @@ contains real(wp) :: Cv, Cp, e_mix, Mw, R_gas real(wp) :: vel_K_sum, vel_dv_dt_sum integer :: i, j, k, r !< Generic loop iterators + type(eos_state) :: eos_s ! Reshaping of inputted data and association of the FD and PI coefficients, or CBC coefficients, respectively, hinging on ! selected CBC coordinate direction @@ -597,7 +598,7 @@ contains & dalpha_rho_ds, dpres_ds, dvel_dt, dadv_dt, dalpha_rho_dt, L, lambda, Ys, dYs_dt, dYs_ds, & & h_k, Cp_i, Gamma_i, Xs, drho_dt, dpres_dt, dpi_inf_dt, dqv_dt, dgamma_dt, rho, pres, E, H, & & gamma, pi_inf, qv, c, Ma, T, sum_Enthalpies, Cv, Cp, e_mix, Mw, R_gas, vel_K_sum, & - & vel_dv_dt_sum, i, j]', copyin='[dir_idx]') + & vel_dv_dt_sum, i, j, eos_s]', copyin='[dir_idx]') do r = is3%beg, is3%end do k = is2%beg, is2%end ! Transferring the Primitive Variables @@ -661,7 +662,8 @@ contains H = (E + pres)/rho ! Compute mixture sound speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c, qv) + call s_eos_state_roe(eos_s, pres, rho, gamma, pi_inf, qv, vel_K_sum, H) + call s_compute_speed_of_sound(eos_s, adv_local, c) ! First-Order Spatial Derivatives of Primitive Variables diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 362c43e33..946e8dbcc 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -193,6 +193,7 @@ contains real(wp) :: Rc_min_loc, Rc_min_glb !< Rc stability extrema on local and global grids real(wp) :: icfl, vcfl, ccfl, Rc integer :: fl !< Fluid loop iterator + type(eos_state) :: eos_s icfl_max_loc = 0._wp vcfl_max_loc = 0._wp @@ -200,14 +201,15 @@ contains Rc_min_loc = huge(1.0_wp) ! Computing Stability Criteria at Current Time-step $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, icfl, & - & vcfl, Rc, ccfl, fl]', reduction='[[icfl_max_loc, vcfl_max_loc, ccfl_max_loc], [Rc_min_loc]]', & - & reductionOp='[max, min]') + & vcfl, Rc, ccfl, fl, eos_s]', reduction='[[icfl_max_loc, vcfl_max_loc, & + & ccfl_max_loc], [Rc_min_loc]]', reductionOp='[max, min]') do l = 0, p do k = 0, n do j = 0, m call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l) - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv) + call s_eos_state_roe(eos_s, pres, rho, gamma, pi_inf, qv, vel_sum, H) + call s_compute_speed_of_sound(eos_s, alpha, c) if (any_non_newtonian) then Re(1) = 0._wp @@ -1162,6 +1164,7 @@ contains real(wp) :: rad, thickness !< For integral quantities logical :: trigger !< For integral quantities real(wp) :: rhoYks(1:num_species) + type(eos_state) :: eos_s T = dflt_T_guess @@ -1299,8 +1302,8 @@ contains end if ! Compute mixture sound Speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, & - & 0._wp, c, qv) + call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp) + call s_compute_speed_of_sound(eos_s, alpha, c) accel = accel_mag(j - 2, k, l) end if @@ -1382,8 +1385,8 @@ contains Rdot(:) = nRdot(:)/nbub end if ! Compute mixture sound speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, & - & 0._wp, 0._wp, c, qv) + call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp) + call s_compute_speed_of_sound(eos_s, alpha, c) end if end if else @@ -1445,8 +1448,8 @@ contains end if ! Compute mixture sound speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, & - & 0._wp, 0._wp, c, qv) + call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp) + call s_compute_speed_of_sound(eos_s, alpha, c) accel = accel_mag(j - 2, k - 2, l - 2) end if diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index 28c0db107..9dfe05d30 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -97,6 +97,7 @@ contains type(riemann_states_vec3) :: cm !< Conservative momentum variables integer :: i, j, k, l !< Generic loop iterators integer :: Re_size_loc1, Re_size_loc2 !< host copies of Re_size; amdflang reads the declare-target original stale cross-TU + type(eos_state) :: eos_s_L, eos_s_R, eos_s_avg ! Populating the buffers of the left and right Riemann problem states variables, based on the choice of boundary conditions call s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, & @@ -120,8 +121,8 @@ contains & Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, & & gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, damage_L, damage_R, & & rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, & - & Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R]', & - & copyin='[norm_dir]', firstprivate='[Re_size_loc1, Re_size_loc2]') + & Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R, eos_s_L, & + & eos_s_R, eos_s_avg]', copyin='[norm_dir]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = ${Z_BND}$%beg, ${Z_BND}$%end do k = ${Y_BND}$%beg, ${Y_BND}$%end do j = ${X_BND}$%beg, ${X_BND}$%end @@ -320,17 +321,17 @@ contains @:compute_average_state() - call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, 0._wp, c_L, & - & qv_L) + eos_s_L = eos_state(rho_L, pres_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L, 0._wp) + call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L) - call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, 0._wp, c_R, & - & qv_R) + eos_s_R = eos_state(rho_R, pres_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R, 0._wp) + call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R) !> The computation of c_avg does not require all the variables, and therefore the non '_avg' ! variables are placeholders to call the subroutine. - call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, vel_avg_rms, & - & c_sum_Yi_Phi, c_avg, qv_avg) + eos_s_avg = eos_state(rho_avg, pres_R, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, c_sum_Yi_Phi) + call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg) if (mhd) then call s_compute_fast_magnetosonic_speed(rho_L, c_L, B%L, norm_dir, c_fast%L, H_L) diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index e760f5874..dfbf5abb9 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -117,6 +117,7 @@ contains real(wp) :: zcoef, pcorr !< low Mach number correction integer :: i, j, k, l, q !< Generic loop iterators integer :: Re_size_loc1, Re_size_loc2 !< host copy of Re_size; amdflang reads the declare-target original stale cross-TU + type(eos_state) :: eos_s_L, eos_s_R, eos_s_avg ! HLLC star-state helpers #:if not MFC_CASE_OPTIMIZATION and USING_AMD @@ -187,7 +188,7 @@ contains & ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, & & Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, rho_Star, E_Star, p_Star, p_K_Star, & & vel_K_star, s_L, s_R, s_M, s_P, s_S, xi_M, xi_P, xi_L, xi_R, xi_L_m1, xi_R_m1, xi_MP, & - & xi_PP]', firstprivate='[Re_size_loc1, Re_size_loc2]') + & xi_PP, eos_s_L, eos_s_R, eos_s_avg]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = ${Z_BND}$%beg, ${Z_BND}$%end do k = ${Y_BND}$%beg, ${Y_BND}$%end do j = ${X_BND}$%beg, ${X_BND}$%end @@ -274,16 +275,16 @@ contains @:compute_average_state() - call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, 0._wp, & - & c_L, qv_L) + eos_s_L = eos_state(rho_L, pres_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L, 0._wp) + call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L) - call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, 0._wp, & - & c_R, qv_R) + eos_s_R = eos_state(rho_R, pres_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R, 0._wp) + call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R) !> The computation of c_avg does not require all the variables, and therefore the non '_avg' ! variables are placeholders to call the subroutine. - call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, vel_avg_rms, & - & 0._wp, c_avg, qv_avg) + eos_s_avg = eos_state(rho_avg, pres_R, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, 0._wp) + call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg) if (viscous) then $:GPU_LOOP(parallelism='[seq]') @@ -474,7 +475,7 @@ contains & alpha_R_sum, s_L, s_R, s_M, s_P, s_S, xi_M, xi_P, xi_L, xi_R, xi_L_m1, xi_R_m1, xi_MP, & & xi_PP, nbub_L, nbub_R, PbwR3Lbar, PbwR3Rbar, R3Lbar, R3Rbar, R3V2Lbar, R3V2Rbar, Ys_L, & & Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Yi_avg, Phi_avg, h_iL, h_iR, & - & h_avg_2]', firstprivate='[Re_size_loc1, Re_size_loc2]') + & h_avg_2, eos_s_L, eos_s_R, eos_s_avg]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = ${Z_BND}$%beg, ${Z_BND}$%end do k = ${Y_BND}$%beg, ${Y_BND}$%end do j = ${X_BND}$%beg, ${X_BND}$%end @@ -646,16 +647,16 @@ contains end do end if - call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, 0._wp, & - & c_L, qv_L) + eos_s_L = eos_state(rho_L, pres_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L, 0._wp) + call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L) - call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, 0._wp, & - & c_R, qv_R) + eos_s_R = eos_state(rho_R, pres_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R, 0._wp) + call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R) !> The computation of c_avg does not require all the variables, and therefore the non '_avg' ! variables are placeholders to call the subroutine. - call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, vel_avg_rms, & - & 0._wp, c_avg, qv_avg) + eos_s_avg = eos_state(rho_avg, pres_R, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, 0._wp) + call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg) if (viscous) then $:GPU_LOOP(parallelism='[seq]') @@ -857,11 +858,11 @@ contains ! Private list split across _hllc_p1/p2/p3 for Fypp line-length limits #:set _hllc_p1 = '[i, j, k, l, q, T_L, T_R, vel_L_rms, vel_R_rms, pres_L, pres_R, rho_L, gamma_L, pi_inf_L, qv_L, rho_R, gamma_R, pi_inf_R, qv_R, alpha_L_sum, alpha_R_sum, E_L, E_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, Gamm_L, Gamm_R, Y_L, Y_R, H_L, H_R, qv_avg, rho_avg, gamma_avg, H_avg, c_L, c_R, c_avg, s_P, s_M, xi_P, xi_M, xi_L, xi_R, xi_L_m1, xi_R_m1, Ms_L, Ms_R, pres_SL, pres_SR, vel_L, vel_R, Re_L, Re_R, alpha_L, alpha_R, alpha_rho_L, alpha_rho_R, alpha_lim_L, alpha_lim_R, s_L, s_R, s_S, vel_avg_rms, pcorr, zcoef, ptilde_L, ptilde_R, vel_L_tmp, vel_R_tmp, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, tau_e_L, tau_e_R, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, G_L, G_R, damage_L, damage_R,' #:set _hllc_p2 = 'U_L, U_R, F_L, F_R, F_star_L, F_star_R, F_HLLC, u_n_HLLC, u_t_HLLC, u_t2_HLLC, pres_tot_L, pres_tot_R, u_n_L, u_n_R, u_t_L, u_t_R, u_t2_L, u_t2_R, tau_nn_L, tau_nn_R, tau_nt_L, tau_nt_R, tau_tt_L, tau_tt_R, tau_nt2_L, tau_nt2_R, tau_t2t2_L, tau_t2t2_R, tau_t1t2_L, tau_t1t2_R, tau_qq_L, tau_qq_R, p_face, tau_qq_face, A_L, A_R, denom_A, u_t_star, tau_nt_star, u_t2_star, tau_nt2_star, pres_tot_star,' - #:set _hllc_p3 = 'F_HLL, u_n_HLL_trace, u_t_HLL_trace, u_t2_HLL_trace, p_face_HLL, tau_qq_face_HLL, tau_nn_HLL, phi, Sigma_L, Sigma_R, dSigma, Sigma_ref, a_L_ref, a_R_ref, a_ref, du_t, dtau_nt, du_t2, dtau_nt2, sensor_ptot, sensor_vt, sensor_tnt, sensor_combined, idx_phys]' + #:set _hllc_p3 = 'F_HLL, u_n_HLL_trace, u_t_HLL_trace, u_t2_HLL_trace, p_face_HLL, tau_qq_face_HLL, tau_nn_HLL, phi, Sigma_L, Sigma_R, dSigma, Sigma_ref, a_L_ref, a_R_ref, a_ref, du_t, dtau_nt, du_t2, dtau_nt2, sensor_ptot, sensor_vt, sensor_tnt, sensor_combined, idx_phys, eos_s_L, eos_s_R, eos_s_avg]' #:set _hllc_priv = _hllc_p1 + _hllc_p2 + _hllc_p3 #:else ! Master's pure-fluid private list, unchanged - #:set _hllc_priv = '[i, T_L, T_R, vel_L_rms, vel_R_rms, pres_L, pres_R, rho_L, gamma_L, pi_inf_L, qv_L, rho_R, gamma_R, pi_inf_R, qv_R, alpha_L_sum, alpha_R_sum, E_L, E_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, Y_L, Y_R, H_L, H_R, qv_avg, rho_avg, gamma_avg, H_avg, c_L, c_R, c_avg, s_P, s_M, xi_P, xi_M, xi_L, xi_R, xi_L_m1, xi_R_m1, Ms_L, Ms_R, pres_SL, pres_SR, vel_L, vel_R, Re_L, Re_R, alpha_L, alpha_R, alpha_rho_L, alpha_rho_R, alpha_lim_L, alpha_lim_R, s_L, s_R, s_S, vel_avg_rms, pcorr, zcoef, vel_L_tmp, vel_R_tmp, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2]' + #:set _hllc_priv = '[i, T_L, T_R, vel_L_rms, vel_R_rms, pres_L, pres_R, rho_L, gamma_L, pi_inf_L, qv_L, rho_R, gamma_R, pi_inf_R, qv_R, alpha_L_sum, alpha_R_sum, E_L, E_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, Y_L, Y_R, H_L, H_R, qv_avg, rho_avg, gamma_avg, H_avg, c_L, c_R, c_avg, s_P, s_M, xi_P, xi_M, xi_L, xi_R, xi_L_m1, xi_R_m1, Ms_L, Ms_R, pres_SL, pres_SR, vel_L, vel_R, Re_L, Re_R, alpha_L, alpha_R, alpha_rho_L, alpha_rho_R, alpha_lim_L, alpha_lim_R, s_L, s_R, s_S, vel_avg_rms, pcorr, zcoef, vel_L_tmp, vel_R_tmp, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, eos_s_L, eos_s_R, eos_s_avg]' #:endif ! The two calls below are identical on purpose. An offload kernel is named ! after the .fpp line of its GPU_PARALLEL_LOOP, so one shared call would give @@ -1059,16 +1060,17 @@ contains @:compute_average_state() - call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, & - & 0._wp, c_L, qv_L) + eos_s_L = eos_state(rho_L, pres_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L, 0._wp) + call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L) - call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, & - & 0._wp, c_R, qv_R) + eos_s_R = eos_state(rho_R, pres_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R, 0._wp) + call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R) !> The computation of c_avg does not require all the variables, and therefore the non '_avg' ! variables are placeholders to call the subroutine. - call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & - & vel_avg_rms, c_sum_Yi_Phi, c_avg, qv_avg) + eos_s_avg = eos_state(rho_avg, pres_R, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, & + & c_sum_Yi_Phi) + call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg) if (viscous) then if (chemistry) then diff --git a/src/simulation/m_riemann_solver_hlld.fpp b/src/simulation/m_riemann_solver_hlld.fpp index c49992bed..37d1b893e 100644 --- a/src/simulation/m_riemann_solver_hlld.fpp +++ b/src/simulation/m_riemann_solver_hlld.fpp @@ -56,10 +56,11 @@ contains ! normal velocity, and x is the normal direction Note: Bx is omitted as the magnetic flux is always zero in the normal ! direction - real(wp) :: sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx - real(wp) :: vL_star, vR_star, wL_star, wR_star - real(wp) :: v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double - integer :: i, j, k, l + real(wp) :: sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx + real(wp) :: vL_star, vR_star, wL_star, wR_star + real(wp) :: v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double + integer :: i, j, k, l + type(eos_state) :: eos_s_L, eos_s_R call s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, & & qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, norm_dir, ix, iy, iz) @@ -78,7 +79,7 @@ contains & U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, F_hlld, s_L, s_R, s_M, s_starL, s_starR, & & pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, & & sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, v_double, w_double, & - & By_double, Bz_double, E_doubleL, E_doubleR, E_double]', copyin='[norm_dir]') + & By_double, Bz_double, E_doubleL, E_doubleR, E_double, eos_s_L, eos_s_R]', copyin='[norm_dir]') do l = ${Z_BND}$%beg, ${Z_BND}$%end do k = ${Y_BND}$%beg, ${Y_BND}$%end do j = ${X_BND}$%beg, ${X_BND}$%end @@ -147,10 +148,10 @@ contains H_no_mag%R = (E%R + pres%R - pres_mag%R)/rho%R ! (2) Compute fast wave speeds - call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, H_no_mag%L, alpha_L, vel_rms%L, & - & 0._wp, c%L, qv%L) - call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, H_no_mag%R, alpha_R, vel_rms%R, & - & 0._wp, c%R, qv%R) + eos_s_L = eos_state(rho%L, pres%L, gamma%L, pi_inf%L, qv%L, vel_rms%L, H_no_mag%L, 0._wp) + call s_compute_speed_of_sound(eos_s_L, alpha_L, c%L) + eos_s_R = eos_state(rho%R, pres%R, gamma%R, pi_inf%R, qv%R, vel_rms%R, H_no_mag%R, 0._wp) + call s_compute_speed_of_sound(eos_s_R, alpha_R, c%R) call s_compute_fast_magnetosonic_speed(rho%L, c%L, B%L, norm_dir, c_fast%L, H_no_mag%L) call s_compute_fast_magnetosonic_speed(rho%R, c%R, B%R, norm_dir, c_fast%R, H_no_mag%R) diff --git a/src/simulation/m_riemann_solver_hypo_hlld.fpp b/src/simulation/m_riemann_solver_hypo_hlld.fpp index 444be1709..3065a57d1 100644 --- a/src/simulation/m_riemann_solver_hypo_hlld.fpp +++ b/src/simulation/m_riemann_solver_hypo_hlld.fpp @@ -103,14 +103,15 @@ contains ! HLLD Hypo variables - real(wp) :: G_eff, G_eff_tol, C_NC, sqrtC_NC - real(wp) :: A_L, A_R, denomA, fac_L, fac_R - real(wp) :: u_n_L, u_t_L, u_n_R, u_t_R - real(wp) :: u_t2_L, u_t2_R - real(wp) :: tau_nn_L, tau_nt_L, tau_tt_L, tau_nn_R, tau_nt_R, tau_tt_R - real(wp) :: tau_nt2_L, tau_nt2_R, tau_t2t2_L, tau_t2t2_R, tau_t1t2_L, tau_t1t2_R - real(wp) :: tau_qq_L, tau_qq_R - real(wp) :: G_L, G_R + real(wp) :: G_eff, G_eff_tol, C_NC, sqrtC_NC + real(wp) :: A_L, A_R, denomA, fac_L, fac_R + real(wp) :: u_n_L, u_t_L, u_n_R, u_t_R + real(wp) :: u_t2_L, u_t2_R + real(wp) :: tau_nn_L, tau_nt_L, tau_tt_L, tau_nn_R, tau_nt_R, tau_tt_R + real(wp) :: tau_nt2_L, tau_nt2_R, tau_t2t2_L, tau_t2t2_R, tau_t1t2_L, tau_t1t2_R + real(wp) :: tau_qq_L, tau_qq_R + real(wp) :: G_L, G_R + type(eos_state) :: eos_s_L, eos_s_R #:if not MFC_CASE_OPTIMIZATION and USING_AMD real(wp), dimension(6) :: tau_e_L, tau_e_R #:else @@ -183,7 +184,7 @@ contains #:set _hlld_p1 = '[i,j,k,l,ipass,degenerate,shear_degenerate,fan_fallback,alpha_rho_L,alpha_rho_R,vel,alpha_L,alpha_R,rho,pres,E,H,gamma,pi_inf,qv,vel_rms,c,S_L,S_R,s_M,S_Lstar,S_Rstar,pTot_L,pTot_R,rhoL_star,rhoR_star,U_L,U_R,F_L,F_R,F_hlld,us_c,uss_c,zone,F_HLL_c,U_HLL_c,rho_HLL,u_n_HLL_cons,tau_nn_HLL,u_n_HLL_trace,u_t_HLL_trace,p_face_HLL,tau_qq_face_HLL,ncomp,G_eff,G_eff_tol,C_NC,sqrtC_NC,A_L,A_R,denomA,fac_L,fac_R,' #:set _hlld_p2 = 'u_n_L,u_t_L,u_n_R,u_t_R,u_t2_L,u_t2_R,tau_nn_L,tau_nt_L,tau_tt_L,tau_nn_R,tau_nt_R,tau_tt_R,tau_nt2_L,tau_nt2_R,tau_t2t2_L,tau_t2t2_R,tau_t1t2_L,tau_t1t2_R,tau_qq_L,tau_qq_R,G_L,G_R,tau_e_L,tau_e_R,alpha1_L_star,alpha1_R_star,alpha2_L_star,alpha2_R_star,u_t_star,tau_nt_star,u_t2_star,tau_nt2_star,tau_nn_L_star,tau_nn_R_star,tau_tt_L_star,tau_tt_R_star,tau_tt_L_starstar,tau_tt_R_starstar,' #:set _hlld_p3 = 'tau_t2t2_L_star,tau_t2t2_R_star,tau_t2t2_L_starstar,tau_t2t2_R_starstar,tau_t1t2_L_star,tau_t1t2_R_star,tau_t1t2_L_starstar,tau_t1t2_R_starstar,tau_qq_L_star,tau_qq_R_star,pTot_star,E_L_star,E_R_star,E_L_starstar,E_R_starstar,p_face,tau_qq_face,u_n_face,u_t_face,G_hat,rho_hat,tau_nn_hat,tau_nt_hat,tau_tt_hat,tau_qq_hat,tau_nt2_hat,tau_t2t2_hat,tau_t1t2_hat,' - #:set _hlld_p4 = 'alpha_hat,alpha_rho_hat,tau_e_hat,pres_hat,blkmod1_hat,blkmod2_hat,K_hat,C_hat_1,C_hat_2,Sigma_L,Sigma_R,dSigma,Sigma_ref,a_L_ref,a_R_ref,a_ref,du_t,dtau_nt,du_t2,dtau_nt2,sensor_ptot,sensor_vt,sensor_tnt,sensor_combined,phi,alpha_L_sum,alpha_R_sum]' + #:set _hlld_p4 = 'alpha_hat,alpha_rho_hat,tau_e_hat,pres_hat,blkmod1_hat,blkmod2_hat,K_hat,C_hat_1,C_hat_2,Sigma_L,Sigma_R,dSigma,Sigma_ref,a_L_ref,a_R_ref,a_ref,du_t,dtau_nt,du_t2,dtau_nt2,sensor_ptot,sensor_vt,sensor_tnt,sensor_combined,phi,alpha_L_sum,alpha_R_sum,eos_s_L,eos_s_R]' ! Wave-fan side table for the per-component F_hlld fold below: side name, the side's two zones, ! its starstar zone, and the outer/inner wave speeds. The L and R sides are mirror images. #:set HLLD_FAN_SIDES = [('L', 1, 2, 2, 'S_L', 'S_Lstar'), ('R', 3, 4, 3, 'S_R', 'S_Rstar')] @@ -357,10 +358,11 @@ contains ! Compute Riemann states - call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, H%L, alpha_L, vel_rms%L, 0._wp, c%L, & - & qv%L) - call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, H%R, alpha_R, vel_rms%R, 0._wp, c%R, & - & qv%R) + eos_s_L = eos_state(rho%L, pres%L, gamma%L, pi_inf%L, qv%L, vel_rms%L, H%L, 0._wp) + call s_compute_speed_of_sound(eos_s_L, alpha_L, c%L) + + eos_s_R = eos_state(rho%R, pres%R, gamma%R, pi_inf%R, qv%R, vel_rms%R, H%R, 0._wp) + call s_compute_speed_of_sound(eos_s_R, alpha_R, c%R) S_L = min(u_n_L - sqrt(max(verysmall, c%L*c%L + ((4._wp/3._wp)*G_L + tau_nn_L)/rho%L)), & & u_n_R - sqrt(max(verysmall, c%R*c%R + ((4._wp/3._wp)*G_R + tau_nn_R)/rho%R))) diff --git a/src/simulation/m_riemann_solver_lf.fpp b/src/simulation/m_riemann_solver_lf.fpp index 512a576e7..9501a8100 100644 --- a/src/simulation/m_riemann_solver_lf.fpp +++ b/src/simulation/m_riemann_solver_lf.fpp @@ -93,6 +93,7 @@ contains integer :: i, j, k, l !< Generic loop iterators integer :: Re_size_loc1, Re_size_loc2 !< host copies of Re_size; amdflang reads the declare-target original stale cross-TU integer, dimension(3) :: idx_right_phys !< Physical (j,k,l) indices for right state. + type(eos_state) :: eos_s_L, eos_s_R ! Populating the buffers of the left and right Riemann problem states variables, based on the choice of boundary conditions call s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, & @@ -116,7 +117,7 @@ contains & c_avg, pres_L, pres_R, rho_L, rho_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, c_L, & & c_R, E_L, E_R, H_L, H_R, ptilde_L, ptilde_R, s_M, s_P, xi_M, xi_P, Cp_avg, Cv_avg, T_avg, & & eps, c_sum_Yi_Phi, Cp_L, Cp_R, Cv_L, Cv_R, R_gas_L, R_gas_R, MW_L, MW_R, T_L, T_R, Y_L, & - & Y_R]', firstprivate='[Re_size_loc1, Re_size_loc2]') + & Y_R, eos_s_L, eos_s_R]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = ${Z_BND}$%beg, ${Z_BND}$%end do k = ${Y_BND}$%beg, ${Y_BND}$%end do j = ${X_BND}$%beg, ${X_BND}$%end @@ -291,11 +292,11 @@ contains H_R = (E_R + pres_R)/rho_R end if - call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, 0._wp, c_L, & - & qv_L) + eos_s_L = eos_state(rho_L, pres_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L, 0._wp) + call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L) - call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, 0._wp, c_R, & - & qv_R) + eos_s_R = eos_state(rho_R, pres_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R, 0._wp) + call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R) if (mhd) then call s_compute_fast_magnetosonic_speed(rho_L, c_L, B%L, norm_dir, c_fast%L, H_L) diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index daf1c2a73..aa685807d 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -659,14 +659,15 @@ contains real(wp) :: dt_local integer :: j, k, l !< Generic loop iterators integer :: fl !< Fluid loop iterator + type(eos_state) :: eos_s if (.not. igr) then call s_convert_conservative_to_primitive_variables(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, idwint) end if dt_local = huge(1.0_wp) - $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, fl, max_dt]', & - & reduction='[[dt_local]]', reductionOp='[min]') + $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, fl, max_dt, & + & eos_s]', reduction='[[dt_local]]', reductionOp='[min]') do l = 0, p do k = 0, n do j = 0, m @@ -677,7 +678,8 @@ contains end if ! Compute mixture sound speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv) + call s_eos_state_roe(eos_s, pres, rho, gamma, pi_inf, qv, vel_sum, H) + call s_compute_speed_of_sound(eos_s, alpha, c) if (any_non_newtonian) then Re(1) = 0._wp