diff --git a/docs/documentation/gpuParallelization.md b/docs/documentation/gpuParallelization.md index 3023ac788..d59b84fc5 100644 --- a/docs/documentation/gpuParallelization.md +++ b/docs/documentation/gpuParallelization.md @@ -614,15 +614,14 @@ Does not do anything for OpenMP currently pure, sequential per-thread helpers is: ```fortran -subroutine s_accumulate_mixture_properties(nf, alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K) +subroutine s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K) - $:GPU_ROUTINE(function_name='s_accumulate_mixture_properties', parallelism='[seq]', cray_inline=True) + $:GPU_ROUTINE(function_name='s_compute_mixture_coefficients', parallelism='[seq]', cray_inline=True) - integer, intent(in) :: nf - real(wp), dimension(nf), intent(in) :: alpha_rho_K, alpha_K - real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K + real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K + real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K ... -end subroutine s_accumulate_mixture_properties +end subroutine s_compute_mixture_coefficients ``` **When to use it.** Extract a block into a `GPU_ROUTINE` helper when: diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 051a5a6d3..5f90e9988 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -25,8 +25,9 @@ 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_finalize_variables_conversion_module, gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps + & s_compute_species_fraction, s_compute_mixture_coefficients, s_compute_energy, s_compute_speed_of_sound, & + & s_compute_speed_of_sound_avg, s_compute_fast_magnetosonic_speed, s_finalize_variables_conversion_module, gammas, & + & gs_min, pi_infs, ps_inf, cvs, qvs, qvps real(wp), allocatable, dimension(:) :: Gs_vc integer, allocatable, dimension(:) :: bubrs_vc @@ -44,7 +45,6 @@ module m_variables_conversion real(wp), allocatable, dimension(:,:,:), public :: rho_sf !< Scalar density function real(wp), allocatable, dimension(:,:,:), public :: gamma_sf !< Scalar sp. heat ratio function real(wp), allocatable, dimension(:,:,:), public :: pi_inf_sf !< Scalar liquid stiffness function - real(wp), allocatable, dimension(:,:,:), public :: qv_sf !< Scalar liquid energy reference function contains @@ -155,7 +155,6 @@ contains rho_sf(i, j, k) = rho gamma_sf(i, j, k) = gamma pi_inf_sf(i, j, k) = pi_inf - qv_sf(i, j, k) = qv end if end subroutine s_convert_mixture_to_mixture_variables @@ -189,7 +188,6 @@ contains rho_sf(k, l, r) = rho gamma_sf(k, l, r) = gamma pi_inf_sf(k, l, r) = pi_inf - qv_sf(k, l, r) = qv end if end subroutine s_convert_species_to_mixture_variables @@ -220,29 +218,16 @@ contains if (present(G_K)) G_K = 0._wp ! Constrain partial densities and volume fractions within physical bounds - if (num_fluids == 1 .and. bubbles_euler) then - rho_K = alpha_rho_K(1) - gamma_K = gammas(1) - pi_inf_K = pi_infs(1) - qv_K = qvs(1) - else - if (mpp_lim) then - alpha_K_sum = 0._wp - do i = 1, num_fluids - alpha_rho_K(i) = max(0._wp, alpha_rho_K(i)) - alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp) - alpha_K_sum = alpha_K_sum + alpha_K(i) - end do - alpha_K = alpha_K/max(alpha_K_sum, sgm_eps) - end if - rho_K = 0._wp; gamma_K = 0._wp; pi_inf_K = 0._wp; qv_K = 0._wp + if (mpp_lim) then + alpha_K_sum = 0._wp do i = 1, num_fluids - rho_K = rho_K + alpha_rho_K(i) - gamma_K = gamma_K + alpha_K(i)*gammas(i) - pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i) - qv_K = qv_K + alpha_rho_K(i)*qvs(i) + alpha_rho_K(i) = max(0._wp, alpha_rho_K(i)) + alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp) + alpha_K_sum = alpha_K_sum + alpha_K(i) end do + alpha_K = alpha_K/max(alpha_K_sum, sgm_eps) end if + call s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K) if (present(G_K)) then G_K = 0._wp @@ -339,18 +324,15 @@ contains allocate (rho_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) allocate (gamma_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) allocate (pi_inf_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) - allocate (qv_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) else allocate (rho_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) allocate (gamma_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) allocate (pi_inf_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) - allocate (qv_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) end if else allocate (rho_sf(-buff_size:m + buff_size,0:0,0:0)) allocate (gamma_sf(-buff_size:m + buff_size,0:0,0:0)) allocate (pi_inf_sf(-buff_size:m + buff_size,0:0,0:0)) - allocate (qv_sf(-buff_size:m + buff_size,0:0,0:0)) end if end if @@ -1192,7 +1174,7 @@ contains !> Deallocate fluid property arrays and post-processing fields allocated during module initialization. impure subroutine s_finalize_variables_conversion_module() - if (allocated(rho_sf)) deallocate (rho_sf, gamma_sf, pi_inf_sf, qv_sf) + if (allocated(rho_sf)) deallocate (rho_sf, gamma_sf, pi_inf_sf) @:DEALLOCATE(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps, Gs_vc) if (allocated(bubrs_vc)) then @@ -1204,33 +1186,88 @@ contains end subroutine s_finalize_variables_conversion_module - !> 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) + !> 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 + integer :: i !< Loop iterator over fluids + + if (num_fluids == 1 .and. bubbles_euler) then + rho_K = alpha_rho_K(1) + gamma_K = gammas(1) + pi_inf_K = pi_infs(1) + qv_K = 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) + gamma_K = gamma_K + alpha_K(i)*gammas(i) + pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i) + qv_K = qv_K + alpha_rho_K(i)*qvs(i) + end do + end if + + end subroutine s_compute_mixture_coefficients + + !> 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 = gamma*pres + pi_inf + 5.e-1_wp*rho*vel_sum + qv + + end subroutine s_compute_energy + + !> Compute the speed of sound of a thermodynamic state. + !! + !! Enthalpy is not an argument: substituting H = ((Gamma + 1)p + Pi + qv)/rho + |u|^2/2 into + !! c^2 = (H - |u|^2/2 - qv/rho)/Gamma leaves c^2 = ((Gamma + 1)p + Pi)/(Gamma rho), so H, |u|^2 + !! and qv all cancel. Averaged states, whose enthalpy is a free input, use + !! s_compute_speed_of_sound_avg. + subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, 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 + 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(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) + 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/(1._wp + (gamma + 1._wp)*pres/rho)) else if (alt_soundspeed) then ! Wood's mixture sound speed via bulk moduli blkmod1 = ((gammas(1) + 1._wp)*pres + pi_infs(1))/gammas(1) @@ -1251,8 +1288,8 @@ contains 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 + else ! Stiffened-gas mixture, with H, |u|^2 and qv cancelled out + c = ((gamma + 1._wp)*pres + pi_inf)/(gamma*rho) end if if (mixture_err .and. c < 0._wp) then @@ -1264,6 +1301,46 @@ contains end subroutine s_compute_speed_of_sound + !> Compute the 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 it, along with |u|^2 and qv. Only the enthalpy-reading branches + !! differ; the rest defer to 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) + + $: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 (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) + 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)) then + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c) + 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 + !> 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_data_output.fpp b/src/post_process/m_data_output.fpp index 1d5ddfedb..9f171bf9b 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -1238,11 +1238,11 @@ contains impure subroutine s_write_energy_data_file(q_prim_vf, q_cons_vf) type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf, q_cons_vf - real(wp) :: Elk, Egk, Elp, Egint, Vb, Vl, pres_av, Et - real(wp) :: rho, pres, dV, tmp, gamma, pi_inf, MaxMa, MaxMa_glb, maxvel, c, Ma, H, qv - real(wp), dimension(num_vels) :: vel - real(wp), dimension(num_fluids) :: adv - integer :: i, j, k, l, s !< looping indices + real(wp) :: Elk, Egk, Elp, Egint, Vb, Vl, pres_av, Et + real(wp) :: rho, pres, dV, tmp, gamma, pi_inf, MaxMa, MaxMa_glb, maxvel, c, Ma + real(wp), dimension(num_vels) :: vel + real(wp), dimension(num_fluids) :: adv + integer :: i, j, k, l, s !< looping indices Egk = 0._wp Elp = 0._wp @@ -1267,7 +1267,6 @@ contains rho = 0._wp gamma = 0._wp pi_inf = 0._wp - qv = 0._wp pres = q_prim_vf(eqn_idx%E)%sf(i, j, k) Egint = Egint + q_prim_vf(eqn_idx%E + 2)%sf(i, j, k)*(gammas(2)*pres)*dV do s = 1, num_vels @@ -1283,12 +1282,9 @@ contains gamma = gamma + adv(l)*gammas(l) pi_inf = pi_inf + adv(l)*pi_infs(l) rho = rho + adv(l)*q_prim_vf(l)%sf(i, j, k) - qv = qv + adv(l)*q_prim_vf(l)%sf(i, j, k)*qvs(l) end do - 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_compute_speed_of_sound(pres, rho, gamma, pi_inf, 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..9e6b1ee6c 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -177,11 +177,11 @@ contains end subroutine s_perform_time_step !> Derive requested flow quantities from primitive variables and write them to the formatted database files. - impure subroutine s_save_data(t_step, varname, pres, c, H) + impure subroutine s_save_data(t_step, varname, pres, c) integer, intent(inout) :: t_step character(LEN=name_len), intent(inout) :: varname - real(wp), intent(inout) :: pres, c, H + real(wp), intent(inout) :: pres, c real(wp), dimension(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end, & & -offset_z%beg:p + offset_z%end) :: liutex_mag @@ -531,10 +531,7 @@ contains pres = q_prim_vf(eqn_idx%E)%sf(i, j, k) - 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_compute_speed_of_sound(pres, rho_sf(i, j, k), gamma_sf(i, j, k), pi_inf_sf(i, j, k), adv, c) out%q_sf(i, j, k) = c end do diff --git a/src/post_process/p_main.fpp b/src/post_process/p_main.fpp index 20883a42a..897a8ebdc 100644 --- a/src/post_process/p_main.fpp +++ b/src/post_process/p_main.fpp @@ -15,7 +15,6 @@ program p_main character(LEN=name_len) :: varname real(wp) :: pres real(wp) :: c - real(wp) :: H real(wp) :: start, finish call s_initialize_mpi_domain() @@ -41,7 +40,7 @@ program p_main call s_perform_time_step(t_step) - call s_save_data(t_step, varname, pres, c, H) + call s_save_data(t_step, varname, pres, c) call cpu_time(finish) diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index bc72c87bd..5b305b4dc 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -501,7 +501,6 @@ contains real(wp) :: rho !< Cell averaged density real(wp) :: pres !< Cell averaged pressure real(wp) :: E !< Cell averaged energy - real(wp) :: H !< Cell averaged enthalpy real(wp) :: gamma !< Cell averaged specific heat ratio real(wp) :: pi_inf !< Cell averaged liquid stiffness real(wp) :: qv !< Cell averaged fluid reference energy @@ -595,7 +594,7 @@ contains ! FD2 or FD4 of RHS at j = 0 $:GPU_PARALLEL_LOOP(collapse=2, private='[r, k, alpha_rho, vel, adv_local, mf, dvel_ds, dadv_ds, Re_cbc, & & 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, & + & h_k, Cp_i, Gamma_i, Xs, drho_dt, dpres_dt, dpi_inf_dt, dqv_dt, dgamma_dt, rho, pres, E, & & 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]') do r = is3%beg, is3%end @@ -655,13 +654,10 @@ contains gamma = 1.0_wp/(Cp/Cv - 1.0_wp) end if else - E = gamma*pres + pi_inf + 5.e-1_wp*rho*vel_K_sum end if - 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_compute_speed_of_sound(pres, rho, gamma, pi_inf, 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 22bde22e1..1f511d24f 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -198,7 +198,7 @@ contains 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_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c) if (any_non_newtonian) then Re(1) = 0._wp @@ -1291,8 +1291,7 @@ 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_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c) if (hypoelasticity) c = sqrt(c*c + (4._wp/3._wp)*G_local/rho) accel = accel_mag(j - 2, k, l) @@ -1375,8 +1374,7 @@ 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_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c) if (hypoelasticity) c = sqrt(c*c + (4._wp/3._wp)*G_local/rho) end if end if @@ -1445,8 +1443,7 @@ 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_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c) if (hypoelasticity) c = sqrt(c*c + (4._wp/3._wp)*G_local/rho) accel = accel_mag(j - 2, k - 2, l - 2) diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index 28c0db107..c5c65368e 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -199,8 +199,8 @@ contains alpha_R = alpha_R/max(alpha_R_sum, sgm_eps) end if - call s_accumulate_mixture_properties(num_fluids, alpha_rho_L, alpha_L, rho_L, gamma_L, pi_inf_L, qv_L) - call s_accumulate_mixture_properties(num_fluids, alpha_rho_R, alpha_R, rho_R, gamma_R, pi_inf_R, qv_R) + call s_compute_mixture_coefficients(alpha_rho_L, alpha_L, rho_L, gamma_L, pi_inf_L, qv_L) + call s_compute_mixture_coefficients(alpha_rho_R, alpha_R, rho_R, gamma_R, pi_inf_R, qv_R) if (viscous) then call s_compute_interface_reynolds(alpha_L, Re_L, Re_size_loc1, Re_size_loc2) @@ -285,15 +285,17 @@ contains pres_mag%L = 0.5_wp*(B%L(1)**2._wp + B%L(2)**2._wp + B%L(3)**2._wp) pres_mag%R = 0.5_wp*(B%R(1)**2._wp + B%R(2)**2._wp + B%R(3)**2._wp) #:endif - E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L + pres_mag%L + call s_compute_energy(pres_L, alpha_rho_L, alpha_L, vel_L_rms, E_L) + E_L = E_L + pres_mag%L ! includes magnetic energy - E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R + pres_mag%R + call s_compute_energy(pres_R, alpha_rho_R, alpha_R, vel_R_rms, E_R) + E_R = E_R + pres_mag%R H_L = (E_L + pres_L - pres_mag%L)/rho_L ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound) H_R = (E_R + pres_R - pres_mag%R)/rho_R else - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R + call s_compute_energy(pres_L, alpha_rho_L, alpha_L, vel_L_rms, E_L) + call s_compute_energy(pres_R, alpha_rho_R, alpha_R, vel_R_rms, E_R) H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R end if @@ -320,17 +322,15 @@ 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) + call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_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) + call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_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) + call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, & + & c_sum_Yi_Phi, 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..e285f0727 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -256,34 +256,30 @@ contains alpha_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%adv%beg + i - 1) end do - call s_accumulate_mixture_properties(num_fluids, alpha_rho_L, alpha_L, rho_L, gamma_L, pi_inf_L, & - & qv_L) - call s_accumulate_mixture_properties(num_fluids, alpha_rho_R, alpha_R, rho_R, gamma_R, pi_inf_R, & - & qv_R) + call s_compute_mixture_coefficients(alpha_rho_L, alpha_L, rho_L, gamma_L, pi_inf_L, qv_L) + call s_compute_mixture_coefficients(alpha_rho_R, alpha_R, rho_R, gamma_R, pi_inf_R, qv_R) if (viscous) then call s_compute_interface_reynolds(alpha_L, Re_L, Re_size_loc1, Re_size_loc2) call s_compute_interface_reynolds(alpha_R, Re_R, Re_size_loc1, Re_size_loc2) end if - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R + call s_compute_energy(pres_L, alpha_rho_L, alpha_L, vel_L_rms, E_L) + call s_compute_energy(pres_R, alpha_rho_R, alpha_R, vel_R_rms, E_R) H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @: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) + call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_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) + call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_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) + call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, & + & H_avg, 0._wp, alpha_R, c_avg) if (viscous) then $:GPU_LOOP(parallelism='[seq]') @@ -502,27 +498,8 @@ contains vel_R_rms = vel_R_rms + vel_R(i)**2._wp end do - ! Retain this in the refactor - if (mpp_lim .and. (num_fluids > 2)) then - call s_accumulate_mixture_properties(num_fluids, alpha_rho_L, alpha_L, rho_L, gamma_L, & - & pi_inf_L, qv_L) - call s_accumulate_mixture_properties(num_fluids, alpha_rho_R, alpha_R, rho_R, gamma_R, & - & pi_inf_R, qv_R) - else if (num_fluids > 2) then - call s_accumulate_mixture_properties(num_fluids - 1, alpha_rho_L, alpha_L, rho_L, gamma_L, & - & pi_inf_L, qv_L) - call s_accumulate_mixture_properties(num_fluids - 1, alpha_rho_R, alpha_R, rho_R, gamma_R, & - & pi_inf_R, qv_R) - else - rho_L = qL_prim_rsx_vf(${SF('')}$, 1) - gamma_L = gammas(1) - pi_inf_L = pi_infs(1) - qv_L = qvs(1) - rho_R = qR_prim_rsx_vf(${SF(' + 1')}$, 1) - gamma_R = gammas(1) - pi_inf_R = pi_infs(1) - qv_R = qvs(1) - end if + call s_compute_mixture_coefficients(alpha_rho_L, alpha_L, rho_L, gamma_L, pi_inf_L, qv_L) + call s_compute_mixture_coefficients(alpha_rho_R, alpha_R, rho_R, gamma_R, pi_inf_R, qv_R) if (viscous) then if (num_fluids == 1) then ! Need to consider case with num_fluids >= 2 @@ -551,8 +528,8 @@ contains pres_L = qL_prim_rsx_vf(${SF('')}$, eqn_idx%E) pres_R = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%E) - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms - E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms + call s_compute_energy(pres_L, alpha_rho_L, alpha_L, vel_L_rms, E_L) + call s_compute_energy(pres_R, alpha_rho_R, alpha_R, vel_R_rms, E_R) H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -646,16 +623,14 @@ 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) + call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_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) + call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_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) + call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, & + & H_avg, 0._wp, alpha_R, c_avg) if (viscous) then $:GPU_LOOP(parallelism='[seq]') @@ -967,10 +942,8 @@ contains alpha_lim_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%E + i) end do - call s_accumulate_mixture_properties(num_fluids, alpha_rho_L, alpha_lim_L, rho_L, gamma_L, & - & pi_inf_L, qv_L) - call s_accumulate_mixture_properties(num_fluids, alpha_rho_R, alpha_lim_R, rho_R, gamma_R, & - & pi_inf_R, qv_R) + call s_compute_mixture_coefficients(alpha_rho_L, alpha_lim_L, rho_L, gamma_L, pi_inf_L, qv_L) + call s_compute_mixture_coefficients(alpha_rho_R, alpha_lim_R, rho_R, gamma_R, pi_inf_R, qv_R) if (viscous) then call s_compute_interface_reynolds(alpha_L, Re_L, Re_size_loc1, Re_size_loc2) @@ -1026,8 +999,8 @@ contains H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R else - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R + call s_compute_energy(pres_L, alpha_rho_L, alpha_lim_L, vel_L_rms, E_L) + call s_compute_energy(pres_R, alpha_rho_R, alpha_lim_R, vel_R_rms, E_R) H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -1059,16 +1032,14 @@ 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) + call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_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) + call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_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) + call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, & + & H_avg, c_sum_Yi_Phi, 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..bce95254d 100644 --- a/src/simulation/m_riemann_solver_hlld.fpp +++ b/src/simulation/m_riemann_solver_hlld.fpp @@ -122,35 +122,22 @@ contains end if end if - ! Sum properties of all fluid components - rho%L = 0._wp; gamma%L = 0._wp; pi_inf%L = 0._wp; qv%L = 0._wp - rho%R = 0._wp; gamma%R = 0._wp; pi_inf%R = 0._wp; qv%R = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - rho%L = rho%L + alpha_rho_L(i) - gamma%L = gamma%L + alpha_L(i)*gammas(i) - pi_inf%L = pi_inf%L + alpha_L(i)*pi_infs(i) - qv%L = qv%L + alpha_rho_L(i)*qvs(i) - - rho%R = rho%R + alpha_rho_R(i) - gamma%R = gamma%R + alpha_R(i)*gammas(i) - pi_inf%R = pi_inf%R + alpha_R(i)*pi_infs(i) - qv%R = qv%R + alpha_rho_R(i)*qvs(i) - end do + call s_compute_mixture_coefficients(alpha_rho_L, alpha_L, rho%L, gamma%L, pi_inf%L, qv%L) + call s_compute_mixture_coefficients(alpha_rho_R, alpha_R, rho%R, gamma%R, pi_inf%R, qv%R) pres_mag%L = 0.5_wp*sum(B%L**2._wp) pres_mag%R = 0.5_wp*sum(B%R**2._wp) - E%L = gamma%L*pres%L + pi_inf%L + 0.5_wp*rho%L*vel_rms%L + qv%L + pres_mag%L - E%R = gamma%R*pres%R + pi_inf%R + 0.5_wp*rho%R*vel_rms%R + qv%R + pres_mag%R ! includes magnetic energy + call s_compute_energy(pres%L, alpha_rho_L, alpha_L, vel_rms%L, E%L) + E%L = E%L + pres_mag%L + call s_compute_energy(pres%R, alpha_rho_R, alpha_R, vel_rms%R, E%R) + E%R = E%R + pres_mag%R ! includes magnetic energy H_no_mag%L = (E%L + pres%L - pres_mag%L)/rho%L ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound) 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) + call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, alpha_L, c%L) + call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%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 52316daa2..ac56d0a55 100644 --- a/src/simulation/m_riemann_solver_hypo_hlld.fpp +++ b/src/simulation/m_riemann_solver_hypo_hlld.fpp @@ -317,21 +317,8 @@ contains pTot_R = pTot_L end if - ! Sum properties of all fluid components - rho%L = 0._wp; gamma%L = 0._wp; pi_inf%L = 0._wp; qv%L = 0._wp - rho%R = 0._wp; gamma%R = 0._wp; pi_inf%R = 0._wp; qv%R = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - rho%L = rho%L + alpha_rho_L(i) - gamma%L = gamma%L + alpha_L(i)*gammas(i) - pi_inf%L = pi_inf%L + alpha_L(i)*pi_infs(i) - qv%L = qv%L + alpha_rho_L(i)*qvs(i) - - rho%R = rho%R + alpha_rho_R(i) - gamma%R = gamma%R + alpha_R(i)*gammas(i) - pi_inf%R = pi_inf%R + alpha_R(i)*pi_infs(i) - qv%R = qv%R + alpha_rho_R(i)*qvs(i) - end do + call s_compute_mixture_coefficients(alpha_rho_L, alpha_L, rho%L, gamma%L, pi_inf%L, qv%L) + call s_compute_mixture_coefficients(alpha_rho_R, alpha_R, rho%R, gamma%R, pi_inf%R, qv%R) G_L = 0._wp; G_R = 0._wp $:GPU_LOOP(parallelism='[seq]') @@ -340,8 +327,8 @@ contains G_R = G_R + alpha_R(i)*Gs_rs(i) end do - E%L = gamma%L*pres%L + pi_inf%L + 5e-1_wp*rho%L*vel_rms%L + qv%L - E%R = gamma%R*pres%R + pi_inf%R + 5e-1_wp*rho%R*vel_rms%R + qv%R + call s_compute_energy(pres%L, alpha_rho_L, alpha_L, vel_rms%L, E%L) + call s_compute_energy(pres%R, alpha_rho_R, alpha_R, vel_rms%R, E%R) ! Freeze the thermal/kinetic enthalpy used by the EOS sound-speed call before ! adding hypoelastic strain energy to the conservative total energy. @@ -357,10 +344,8 @@ 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) + call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, alpha_L, c%L) + call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%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 1ea7ffa12..694fc2e15 100644 --- a/src/simulation/m_riemann_solver_lf.fpp +++ b/src/simulation/m_riemann_solver_lf.fpp @@ -166,8 +166,8 @@ contains alpha_R = alpha_R/max(alpha_R_sum, sgm_eps) end if - call s_accumulate_mixture_properties(num_fluids, alpha_rho_L, alpha_L, rho_L, gamma_L, pi_inf_L, qv_L) - call s_accumulate_mixture_properties(num_fluids, alpha_rho_R, alpha_R, rho_R, gamma_R, pi_inf_R, qv_R) + call s_compute_mixture_coefficients(alpha_rho_L, alpha_L, rho_L, gamma_L, pi_inf_L, qv_L) + call s_compute_mixture_coefficients(alpha_rho_R, alpha_R, rho_R, gamma_R, pi_inf_R, qv_R) if (viscous) then call s_compute_interface_reynolds(alpha_L, Re_L, Re_size_loc1, Re_size_loc2) @@ -223,17 +223,15 @@ contains H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R else - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R + call s_compute_energy(pres_L, alpha_rho_L, alpha_L, vel_L_rms, E_L) + call s_compute_energy(pres_R, alpha_rho_R, alpha_R, vel_R_rms, E_R) H_L = (E_L + pres_L)/rho_L 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) + call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_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) + call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R) s_L = 0._wp; s_R = 0._wp diff --git a/src/simulation/m_riemann_state.fpp b/src/simulation/m_riemann_state.fpp index abb1733b7..bd978308f 100644 --- a/src/simulation/m_riemann_state.fpp +++ b/src/simulation/m_riemann_state.fpp @@ -1030,33 +1030,6 @@ contains end subroutine s_calculate_bulk_stress_tensor - !> Accumulate the mixture density, specific heat ratio function, liquid stiffness function, and internal energy reference of one - !! Riemann state from its partial densities and volume fractions. The number of fluids is an explicit argument because the - !! 5-equation bubble model accumulates over num_fluids - 1 fluids. - subroutine s_accumulate_mixture_properties(nf, alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K) - - $:GPU_ROUTINE(function_name='s_accumulate_mixture_properties', parallelism='[seq]', cray_inline=True) - - integer, intent(in) :: nf !< Number of fluids to accumulate over - real(wp), dimension(nf), intent(in) :: alpha_rho_K, alpha_K - real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K - integer :: i !< Loop iterator over fluids - - rho_K = 0._wp - gamma_K = 0._wp - pi_inf_K = 0._wp - qv_K = 0._wp - - $:GPU_LOOP(parallelism='[seq]') - do i = 1, nf - rho_K = rho_K + alpha_rho_K(i) - gamma_K = gamma_K + alpha_K(i)*gammas(i) - pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i) - qv_K = qv_K + alpha_rho_K(i)*qvs(i) - end do - - end subroutine s_accumulate_mixture_properties - !> Compute the shear and volume Reynolds numbers of one Riemann state by inverse-weighting the fluid Reynolds numbers with the !! volume fractions. subroutine s_compute_interface_reynolds(alpha_K, Re_K, Re_size_loc1, Re_size_loc2) diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index daf1c2a73..775acf33a 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -677,7 +677,7 @@ 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_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c) if (any_non_newtonian) then Re(1) = 0._wp