From 5d977c3358f008e2842022f7c85a404382d78272 Mon Sep 17 00:00:00 2001 From: Hyeoksu Lee Date: Fri, 28 Nov 2025 16:31:25 -0800 Subject: [PATCH 1/7] integrating mixture routines --- src/common/m_checker_common.fpp | 12 +- src/common/m_variables_conversion.fpp | 262 ++++---------------------- src/simulation/m_cbc.fpp | 6 +- src/simulation/m_ibm.fpp | 3 - src/simulation/m_sim_helpers.fpp | 2 - 5 files changed, 41 insertions(+), 244 deletions(-) diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp index 1ec314536a..8cbc9a7e51 100644 --- a/src/common/m_checker_common.fpp +++ b/src/common/m_checker_common.fpp @@ -312,9 +312,9 @@ contains integer :: i bub_fac = 0 - if (bubbles_euler .and. (num_fluids == 1)) bub_fac = 1 + if (bubbles_euler) bub_fac = 1 - do i = 1, num_fluids + do i = 1, num_fluids + bub_fac call s_int_to_str(i, iStr) @:PROHIBIT(.not. f_is_default(fluid_pp(i)%gamma) .and. fluid_pp(i)%gamma <= 0._wp, & "fluid_pp("//trim(iStr)//")%gamma must be positive") @@ -322,20 +322,12 @@ contains @:PROHIBIT(model_eqns == 1 .and. (.not. f_is_default(fluid_pp(i)%gamma)), & "model_eqns = 1 does not support fluid_pp("//trim(iStr)//")%gamma") - @:PROHIBIT((i <= num_fluids + bub_fac .and. fluid_pp(i)%gamma <= 0._wp) .or. & - (i > num_fluids + bub_fac .and. (.not. f_is_default(fluid_pp(i)%gamma))), & - "for fluid_pp("//trim(iStr)//")%gamma") - @:PROHIBIT(.not. f_is_default(fluid_pp(i)%pi_inf) .and. fluid_pp(i)%pi_inf < 0._wp, & "fluid_pp("//trim(iStr)//")%pi_inf must be non-negative") @:PROHIBIT(model_eqns == 1 .and. (.not. f_is_default(fluid_pp(i)%pi_inf)), & "model_eqns = 1 does not support fluid_pp("//trim(iStr)//")%pi_inf") - @:PROHIBIT((i <= num_fluids + bub_fac .and. fluid_pp(i)%pi_inf < 0._wp) .or. & - (i > num_fluids + bub_fac .and. (.not. f_is_default(fluid_pp(i)%pi_inf))), & - "for fluid_pp("//trim(iStr)//")%pi_inf") - @:PROHIBIT(fluid_pp(i)%cv < 0._wp, & "fluid_pp("//trim(iStr)//")%cv must be positive") end do diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 307406e6bd..741b08d4f9 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -33,8 +33,6 @@ module m_variables_conversion s_initialize_mv, & s_convert_to_mixture_variables, & s_convert_mixture_to_mixture_variables, & - s_convert_species_to_mixture_variables_bubbles, & - s_convert_species_to_mixture_variables_bubbles_acc, & s_convert_species_to_mixture_variables, & s_convert_species_to_mixture_variables_acc, & s_convert_conservative_to_primitive_variables, & @@ -93,11 +91,7 @@ contains call s_convert_mixture_to_mixture_variables(q_vf, i, j, k, & rho, gamma, pi_inf, qv) - else if (bubbles_euler) then - call s_convert_species_to_mixture_variables_bubbles(q_vf, i, j, k, & - rho, gamma, pi_inf, qv, Re_K) - else - ! Volume fraction model + else ! Volume fraction model call s_convert_species_to_mixture_variables(q_vf, i, j, k, & rho, gamma, pi_inf, qv, Re_K, G_K, G) end if @@ -230,129 +224,6 @@ contains end subroutine s_convert_mixture_to_mixture_variables - !> This procedure is used alongside with the gamma/pi_inf - !! model to transfer the density, the specific heat ratio - !! function and liquid stiffness function from the vector - !! of conservative or primitive variables to their scalar - !! counterparts. Specifically designed for when subgrid bubbles_euler - !! must be included. - !! @param q_vf primitive variables - !! @param j Cell index - !! @param k Cell index - !! @param l Cell index - !! @param rho density - !! @param gamma specific heat ratio - !! @param pi_inf liquid stiffness - !! @param qv fluid reference energy - subroutine s_convert_species_to_mixture_variables_bubbles(q_vf, j, k, l, & - rho, gamma, pi_inf, qv, Re_K) - - type(scalar_field), dimension(sys_size), intent(in) :: q_vf - - integer, intent(in) :: j, k, l - - real(wp), intent(out), target :: rho - real(wp), intent(out), target :: gamma - real(wp), intent(out), target :: pi_inf - real(wp), intent(out), target :: qv - - real(wp), optional, dimension(2), intent(out) :: Re_K - - integer :: i, q - real(wp), dimension(num_fluids) :: alpha_rho_K, alpha_K - - ! Constraining the partial densities and the volume fractions within - ! their physical bounds to make sure that any mixture variables that - ! are derived from them result within the limits that are set by the - ! fluids physical parameters that make up the mixture - do i = 1, num_fluids - alpha_rho_K(i) = q_vf(i)%sf(j, k, l) - alpha_K(i) = q_vf(advxb + i - 1)%sf(j, k, l) - end do - - if (mpp_lim) then - - 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) - end do - - alpha_K = alpha_K/max(sum(alpha_K), 1.e-16_wp) - - end if - - ! Performing the transfer of the density, the specific heat ratio - ! function as well as the liquid stiffness function, respectively - - if (model_eqns == 4) then - rho = q_vf(1)%sf(j, k, l) - gamma = fluid_pp(1)%gamma !qK_vf(gamma_idx)%sf(i,j,k) - pi_inf = fluid_pp(1)%pi_inf !qK_vf(pi_inf_idx)%sf(i,j,k) - qv = fluid_pp(1)%qv - else if ((model_eqns == 2) .and. bubbles_euler) then - rho = 0._wp; gamma = 0._wp; pi_inf = 0._wp; qv = 0._wp - - if (mpp_lim .and. (num_fluids > 2)) then - do i = 1, num_fluids - rho = rho + q_vf(i)%sf(j, k, l) - gamma = gamma + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%gamma - pi_inf = pi_inf + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%pi_inf - qv = qv + q_vf(i)%sf(j, k, l)*fluid_pp(i)%qv - end do - else if (num_fluids == 2) then - rho = q_vf(1)%sf(j, k, l) - gamma = fluid_pp(1)%gamma - pi_inf = fluid_pp(1)%pi_inf - qv = fluid_pp(1)%qv - else if (num_fluids > 2) then - !TODO: This may need fixing for hypo + bubbles_euler - do i = 1, num_fluids - 1 !leave out bubble part of mixture - rho = rho + q_vf(i)%sf(j, k, l) - gamma = gamma + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%gamma - pi_inf = pi_inf + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%pi_inf - qv = qv + q_vf(i)%sf(j, k, l)*fluid_pp(i)%qv - end do - ! rho = qK_vf(1)%sf(j,k,l) - ! gamma_K = fluid_pp(1)%gamma - ! pi_inf_K = fluid_pp(1)%pi_inf - else - rho = q_vf(1)%sf(j, k, l) - gamma = fluid_pp(1)%gamma - pi_inf = fluid_pp(1)%pi_inf - qv = fluid_pp(1)%qv - end if - end if - -#ifdef MFC_SIMULATION - ! Computing the shear and bulk Reynolds numbers from species analogs - if (viscous) then - if (num_fluids == 1) then ! need to consider case with num_fluids >= 2 - do i = 1, 2 - - Re_K(i) = dflt_real; if (Re_size(i) > 0) Re_K(i) = 0._wp - - do q = 1, Re_size(i) - Re_K(i) = (1 - alpha_K(Re_idx(i, q)))/fluid_pp(Re_idx(i, q))%Re(i) & - + Re_K(i) - end do - - Re_K(i) = 1._wp/max(Re_K(i), sgm_eps) - - end do - end if - end if -#endif - - ! Post process requires rho_sf/gamma_sf/pi_inf_sf/qv_sf to also be updated -#ifdef MFC_POST_PROCESS - rho_sf(j, k, l) = rho - gamma_sf(j, k, l) = gamma - pi_inf_sf(j, k, l) = pi_inf - qv_sf(j, k, l) = qv -#endif - - end subroutine s_convert_species_to_mixture_variables_bubbles - !> This subroutine is designed for the volume fraction model !! and provided a set of either conservative or primitive !! variables, computes the density, the specific heat ratio @@ -417,7 +288,6 @@ contains end do alpha_K = alpha_K/max(sum(alpha_K), 1.e-16_wp) - end if ! Calculating the density, the specific heat ratio function, the @@ -425,26 +295,36 @@ contains ! respectively, from the species analogs rho = 0._wp; gamma = 0._wp; pi_inf = 0._wp; qv = 0._wp - do i = 1, num_fluids - rho = rho + alpha_rho_K(i) - gamma = gamma + alpha_K(i)*gammas(i) - pi_inf = pi_inf + alpha_K(i)*pi_infs(i) - qv = qv + alpha_rho_K(i)*qvs(i) - end do + if (bubbles_euler) then + rho = alpha_rho_K(1) + gamma = gammas(1) + pi_inf = pi_infs(1) + qv = qvs(1) + else + do i = 1, num_fluids + rho = rho + alpha_rho_K(i) + gamma = gamma + alpha_K(i)*gammas(i) + pi_inf = pi_inf + alpha_K(i)*pi_infs(i) + qv = qv + alpha_rho_K(i)*qvs(i) + end do + end if + #ifdef MFC_SIMULATION ! Computing the shear and bulk Reynolds numbers from species analogs - do i = 1, 2 + if (viscous) then + do i = 1, 2 - Re_K(i) = dflt_real; if (Re_size(i) > 0) Re_K(i) = 0._wp + Re_K(i) = dflt_real; if (Re_size(i) > 0) Re_K(i) = 0._wp - do j = 1, Re_size(i) - Re_K(i) = alpha_K(Re_idx(i, j))/fluid_pp(Re_idx(i, j))%Re(i) & - + Re_K(i) - end do + do j = 1, Re_size(i) + Re_K(i) = alpha_K(Re_idx(i, j))/fluid_pp(Re_idx(i, j))%Re(i) & + + Re_K(i) + end do - Re_K(i) = 1._wp/max(Re_K(i), sgm_eps) + Re_K(i) = 1._wp/max(Re_K(i), sgm_eps) - end do + end do + end if #endif if (present(G_K)) then @@ -504,15 +384,21 @@ contains end do alpha_K = alpha_K/max(alpha_K_sum, sgm_eps) - end if - 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 + if (bubbles_euler) then + rho_K = alpha_rho_K(1) + gamma_K = gammas(1) + pi_inf_K = pi_infs(1) + qv_K = qvs(1) + else + 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 if (present(G_K)) then G_K = 0._wp @@ -525,7 +411,6 @@ contains end if if (viscous) then - do i = 1, 2 Re_K(i) = dflt_real @@ -537,77 +422,12 @@ contains end do Re_K(i) = 1._wp/max(Re_K(i), sgm_eps) - end do end if #endif end subroutine s_convert_species_to_mixture_variables_acc - subroutine s_convert_species_to_mixture_variables_bubbles_acc(rho_K, & - gamma_K, pi_inf_K, qv_K, & - alpha_K, alpha_rho_K, Re_K) - $:GPU_ROUTINE(function_name='s_convert_species_to_mixture_variables_bubbles_acc', & - & parallelism='[seq]', cray_inline=True) - - real(wp), intent(inout) :: rho_K, gamma_K, pi_inf_K, qv_K - - real(wp), dimension(num_fluids), intent(in) :: alpha_K, alpha_rho_K !< - !! Partial densities and volume fractions - - real(wp), dimension(2), intent(out) :: Re_K - - integer :: i, j !< Generic loop iterators - -#ifdef MFC_SIMULATION - rho_K = 0._wp - gamma_K = 0._wp - pi_inf_K = 0._wp - qv_K = 0._wp - - if (mpp_lim .and. (model_eqns == 2) .and. (num_fluids > 2)) then - 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 - else if ((model_eqns == 2) .and. (num_fluids > 2)) then - do i = 1, num_fluids - 1 - 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 - else - rho_K = alpha_rho_K(1) - gamma_K = gammas(1) - pi_inf_K = pi_infs(1) - qv_K = qvs(1) - end if - - if (viscous) then - if (num_fluids == 1) then ! need to consider case with num_fluids >= 2 - - do i = 1, 2 - Re_K(i) = dflt_real - - if (Re_size(i) > 0) Re_K(i) = 0._wp - - do j = 1, Re_size(i) - Re_K(i) = (1._wp - alpha_K(Re_idx(i, j)))/Res_vc(i, j) & - + Re_K(i) - end do - - Re_K(i) = 1._wp/max(Re_K(i), sgm_eps) - - end do - end if - end if -#endif - - end subroutine s_convert_species_to_mixture_variables_bubbles_acc - !> The computation of parameters, the allocation of memory, !! the association of pointers and/or the execution of any !! other procedures that are necessary to setup the module. @@ -883,9 +703,6 @@ contains if (elasticity) then call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, & alpha_rho_K, Re_K, G_K, Gs_vc) - else if (bubbles_euler) then - call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, pi_inf_K, qv_K, & - alpha_K, alpha_rho_K, Re_K) else call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, & alpha_K, alpha_rho_K, Re_K) @@ -1495,9 +1312,6 @@ contains call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, & alpha_K, alpha_rho_K, Re_K, & G_K, Gs_vc) - else if (bubbles_euler) then - call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, & - pi_inf_K, qv_K, alpha_K, alpha_rho_K, Re_K) else call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, & alpha_K, alpha_rho_K, Re_K) diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index 93bf4ae64a..b5027b9fd1 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -813,11 +813,7 @@ contains adv_local(i) = q_prim_rs${XYZ}$_vf(0, k, r, E_idx + i) end do - if (bubbles_euler) then - call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, adv_local, alpha_rho, Re_cbc) - else - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, adv_local, alpha_rho, Re_cbc) - end if + call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, adv_local, alpha_rho, Re_cbc) $:GPU_LOOP(parallelism='[seq]') do i = 1, contxe diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 1f3c9a6237..aff2a1882a 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -249,9 +249,6 @@ contains if (elasticity) then call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, & alpha_rho_IP, Re_K, G_K, Gs) - else if (bubbles_euler) then - call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv_K, alpha_IP, & - alpha_rho_IP, Re_K) else call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, & alpha_rho_IP, Re_K) diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index 4e42768b21..054c1f1222 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -134,8 +134,6 @@ contains if (elasticity) then call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, & alpha_rho, Re, G_local, Gs) - elseif (bubbles_euler) then - call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re) else call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re) end if From ae827a53f61f4492a753b782c63e7f58c5c97e64 Mon Sep 17 00:00:00 2001 From: Hyeoksu Lee Date: Fri, 28 Nov 2025 21:57:27 -0800 Subject: [PATCH 2/7] update mixture rules --- src/common/m_variables_conversion.fpp | 176 +++++++++----------------- src/simulation/m_bubbles_EE.fpp | 26 +--- src/simulation/m_sim_helpers.fpp | 24 +--- 3 files changed, 67 insertions(+), 159 deletions(-) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 741b08d4f9..582dc89ce4 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -39,6 +39,7 @@ module m_variables_conversion s_convert_primitive_to_conservative_variables, & s_convert_primitive_to_flux_variables, & s_compute_pressure, & + s_compute_species_fraction, & #ifndef MFC_PRE_PROCESS s_compute_speed_of_sound, & s_compute_fast_magnetosonic_speed, & @@ -250,70 +251,31 @@ contains real(wp), intent(out), target :: qv real(wp), optional, dimension(2), intent(out) :: Re_K - !! Partial densities and volume fractions real(wp), optional, intent(out) :: G_K real(wp), optional, dimension(num_fluids), intent(in) :: G - real(wp), dimension(num_fluids) :: alpha_rho_K, alpha_K !< integer :: i, j !< Generic loop iterator ! Computing the density, the specific heat ratio function and the ! liquid stiffness function, respectively - - if (igr) then - if (num_fluids == 1) then - alpha_rho_K(1) = q_vf(contxb)%sf(k, l, r) - alpha_K(1) = 1._wp - else - do i = 1, num_fluids - 1 - alpha_rho_K(i) = q_vf(i)%sf(k, l, r) - alpha_K(i) = q_vf(advxb + i - 1)%sf(k, l, r) - end do - - alpha_rho_K(num_fluids) = q_vf(num_fluids)%sf(k, l, r) - alpha_K(num_fluids) = 1._wp - sum(alpha_K(1:num_fluids - 1)) - end if - else - do i = 1, num_fluids - alpha_rho_K(i) = q_vf(i)%sf(k, l, r) - alpha_K(i) = q_vf(advxb + i - 1)%sf(k, l, r) - end do - end if - - if (mpp_lim) then - 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) - end do - - alpha_K = alpha_K/max(sum(alpha_K), 1.e-16_wp) - end if + call s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K) ! Calculating the density, the specific heat ratio function, the ! liquid stiffness function, and the energy reference function, ! respectively, from the species analogs rho = 0._wp; gamma = 0._wp; pi_inf = 0._wp; qv = 0._wp - - if (bubbles_euler) then - rho = alpha_rho_K(1) - gamma = gammas(1) - pi_inf = pi_infs(1) - qv = qvs(1) - else - do i = 1, num_fluids - rho = rho + alpha_rho_K(i) - gamma = gamma + alpha_K(i)*gammas(i) - pi_inf = pi_inf + alpha_K(i)*pi_infs(i) - qv = qv + alpha_rho_K(i)*qvs(i) - end do - end if + do i = 1, num_fluids + rho = rho + alpha_rho_K(i) + gamma = gamma + alpha_K(i)*gammas(i) + pi_inf = pi_inf + alpha_K(i)*pi_infs(i) + qv = qv + alpha_rho_K(i)*qvs(i) + end do #ifdef MFC_SIMULATION ! Computing the shear and bulk Reynolds numbers from species analogs if (viscous) then do i = 1, 2 - Re_K(i) = dflt_real; if (Re_size(i) > 0) Re_K(i) = 0._wp do j = 1, Re_size(i) @@ -322,7 +284,6 @@ contains end do Re_K(i) = 1._wp/max(Re_K(i), sgm_eps) - end do end if #endif @@ -353,52 +314,25 @@ contains & parallelism='[seq]', cray_inline=True) real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K - real(wp), dimension(num_fluids), intent(inout) :: alpha_rho_K, alpha_K !< real(wp), dimension(2), intent(out) :: Re_K - !! Partial densities and volume fractions - real(wp), optional, intent(out) :: G_K real(wp), optional, dimension(num_fluids), intent(in) :: G integer :: i, j !< Generic loop iterators - real(wp) :: alpha_K_sum #ifdef MFC_SIMULATION ! Constraining the partial densities and the volume fractions within ! their physical bounds to make sure that any mixture variables that ! are derived from them result within the limits that are set by the ! fluids physical parameters that make up the mixture - rho_K = 0._wp - gamma_K = 0._wp - pi_inf_K = 0._wp - qv_K = 0._wp - - alpha_K_sum = 0._wp - - if (mpp_lim) then - 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 - - if (bubbles_euler) then - rho_K = alpha_rho_K(1) - gamma_K = gammas(1) - pi_inf_K = pi_infs(1) - qv_K = qvs(1) - else - 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 + rho_K = 0._wp; gamma_K = 0._wp; pi_inf_K = 0._wp; qv_K = 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) + end do if (present(G_K)) then G_K = 0._wp @@ -437,7 +371,6 @@ contains $:GPU_ENTER_DATA(copyin='[is1b,is1e,is2b,is2e,is3b,is3e]') -#ifdef MFC_SIMULATION @:ALLOCATE(gammas (1:num_fluids)) @:ALLOCATE(gs_min (1:num_fluids)) @:ALLOCATE(pi_infs(1:num_fluids)) @@ -446,16 +379,6 @@ contains @:ALLOCATE(qvs (1:num_fluids)) @:ALLOCATE(qvps (1:num_fluids)) @:ALLOCATE(Gs_vc (1:num_fluids)) -#else - @:ALLOCATE(gammas (1:num_fluids)) - @:ALLOCATE(gs_min (1:num_fluids)) - @:ALLOCATE(pi_infs(1:num_fluids)) - @:ALLOCATE(ps_inf(1:num_fluids)) - @:ALLOCATE(cvs (1:num_fluids)) - @:ALLOCATE(qvs (1:num_fluids)) - @:ALLOCATE(qvps (1:num_fluids)) - @:ALLOCATE(Gs_vc (1:num_fluids)) -#endif do i = 1, num_fluids gammas(i) = fluid_pp(i)%gamma @@ -484,12 +407,7 @@ contains #endif if (bubbles_euler) then -#ifdef MFC_SIMULATION @:ALLOCATE(bubrs_vc(1:nb)) -#else - @:ALLOCATE(bubrs_vc(1:nb)) -#endif - do i = 1, nb bubrs_vc(i) = bub_idx%rs(i) end do @@ -675,27 +593,7 @@ contains do j = ibounds(1)%beg, ibounds(1)%end dyn_pres_K = 0._wp - if (igr) then - if (num_fluids == 1) then - alpha_rho_K(1) = qK_cons_vf(contxb)%sf(j, k, l) - alpha_K(1) = 1._wp - else - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - 1 - alpha_rho_K(i) = qK_cons_vf(i)%sf(j, k, l) - alpha_K(i) = qK_cons_vf(advxb + i - 1)%sf(j, k, l) - end do - - alpha_rho_K(num_fluids) = qK_cons_vf(num_fluids)%sf(j, k, l) - alpha_K(num_fluids) = 1._wp - sum(alpha_K(1:num_fluids - 1)) - end if - else - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - alpha_rho_K(i) = qK_cons_vf(i)%sf(j, k, l) - alpha_K(i) = qK_cons_vf(advxb + i - 1)%sf(j, k, l) - end do - end if + call s_compute_species_fraction(qK_cons_vf, j, k, l, alpha_rho_K, alpha_K) if (model_eqns /= 4) then #ifdef MFC_SIMULATION @@ -1296,6 +1194,7 @@ contains do i = advxb, advxe alpha_K(i - E_idx) = qK_prim_vf(j, k, l, i) end do + $:GPU_LOOP(parallelism='[seq]') do i = 1, num_vels vel_K(i) = qK_prim_vf(j, k, l, contxe + i) @@ -1389,6 +1288,47 @@ contains #endif end subroutine s_convert_primitive_to_flux_variables + subroutine s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K) + $:GPU_ROUTINE(function_name='s_compute_species_fraction', & + & parallelism='[seq]', cray_inline=True) + type(scalar_field), dimension(sys_size), intent(in) :: q_vf + integer, intent(in) :: k, l, r + real(wp), dimension(num_fluids), intent(out) :: alpha_rho_K, alpha_K + integer :: i + + if (num_fluids == 1) then + alpha_rho_K(1) = q_vf(contxb)%sf(k, l, r) + if (igr .or. bubbles_euler) then + alpha_K(1) = 1._wp + else + alpha_K(1) = q_vf(advxb)%sf(k, l, r) + end if + else + if (igr) then + do i = 1, num_fluids - 1 + alpha_rho_K(i) = q_vf(i)%sf(k, l, r) + alpha_K(i) = q_vf(advxb + i - 1)%sf(k, l, r) + end do + alpha_rho_K(num_fluids) = q_vf(num_fluids)%sf(k, l, r) + alpha_K(num_fluids) = 1._wp - sum(alpha_K(1:num_fluids - 1)) + else + do i = 1, num_fluids + alpha_rho_K(i) = q_vf(i)%sf(k, l, r) + alpha_K(i) = q_vf(advxb + i - 1)%sf(k, l, r) + end do + end if + end if + + if (mpp_lim) then + 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) + end do + alpha_K = alpha_K/max(sum(alpha_K), 1.e-16_wp) + end if + + end subroutine s_compute_species_fraction + impure subroutine s_finalize_variables_conversion_module() ! Deallocating the density, the specific heat ratio function and the diff --git a/src/simulation/m_bubbles_EE.fpp b/src/simulation/m_bubbles_EE.fpp index d8bed3ad16..d7a5b3f376 100644 --- a/src/simulation/m_bubbles_EE.fpp +++ b/src/simulation/m_bubbles_EE.fpp @@ -246,30 +246,16 @@ contains n_tait = 0._wp B_tait = 0._wp - if (mpp_lim .and. (num_fluids > 2)) then - $:GPU_LOOP(parallelism='[seq]') - do ii = 1, num_fluids - myRho = myRho + myalpha_rho(ii) - n_tait = n_tait + myalpha(ii)*gammas(ii) - B_tait = B_tait + myalpha(ii)*pi_infs(ii) - end do - else if (num_fluids > 2) then - $:GPU_LOOP(parallelism='[seq]') - do ii = 1, num_fluids - 1 - myRho = myRho + myalpha_rho(ii) - n_tait = n_tait + myalpha(ii)*gammas(ii) - B_tait = B_tait + myalpha(ii)*pi_infs(ii) - end do - else - myRho = myalpha_rho(1) - n_tait = gammas(1) - B_tait = pi_infs(1)/pi_fac - end if + $:GPU_LOOP(parallelism='[seq]') + do ii = 1, num_fluids + myRho = myRho + myalpha_rho(ii) + n_tait = n_tait + myalpha(ii)*gammas(ii) + B_tait = B_tait + myalpha(ii)*pi_infs(ii)/pi_fac + end do n_tait = 1._wp/n_tait + 1._wp !make this the usual little 'gamma' B_tait = B_tait*(n_tait - 1)/n_tait ! make this the usual pi_inf - myRho = q_prim_vf(1)%sf(j, k, l) myP = q_prim_vf(E_idx)%sf(j, k, l) alf = q_prim_vf(alf_idx)%sf(j, k, l) myR = q_prim_vf(rs(q))%sf(j, k, l) diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index 054c1f1222..d4c7e5a87c 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -109,27 +109,7 @@ contains integer :: i - if (igr) then - if (num_fluids == 1) then - alpha_rho(1) = q_prim_vf(contxb)%sf(j, k, l) - alpha(1) = 1._wp - else - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - 1 - alpha_rho(i) = q_prim_vf(i)%sf(j, k, l) - alpha(i) = q_prim_vf(advxb + i - 1)%sf(j, k, l) - end do - - alpha_rho(num_fluids) = q_prim_vf(num_fluids)%sf(j, k, l) - alpha(num_fluids) = 1._wp - sum(alpha(1:num_fluids - 1)) - end if - else - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - alpha_rho(i) = q_prim_vf(i)%sf(j, k, l) - alpha(i) = q_prim_vf(advxb + i - 1)%sf(j, k, l) - end do - end if + call s_compute_species_fraction(q_prim_vf, j, k, l, alpha_rho, alpha) if (elasticity) then call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, & @@ -138,6 +118,8 @@ contains call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re) end if + if (num_fluids == 1 .and. bubbles_euler) alpha = q_prim_vf(alf_idx)%sf(j, k, l) + if (igr) then $:GPU_LOOP(parallelism='[seq]') do i = 1, num_vels From e1e953cf35e81acbb294a67b3c663b8b2bfc8735 Mon Sep 17 00:00:00 2001 From: Hyeoksu Lee Date: Fri, 28 Nov 2025 22:09:19 -0800 Subject: [PATCH 3/7] fix mixture rule --- src/simulation/m_bubbles_EE.fpp | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/src/simulation/m_bubbles_EE.fpp b/src/simulation/m_bubbles_EE.fpp index d7a5b3f376..9aa5645907 100644 --- a/src/simulation/m_bubbles_EE.fpp +++ b/src/simulation/m_bubbles_EE.fpp @@ -242,16 +242,22 @@ contains myalpha(ii) = q_cons_vf(advxb + ii - 1)%sf(j, k, l) end do - myRho = 0._wp - n_tait = 0._wp - B_tait = 0._wp - - $:GPU_LOOP(parallelism='[seq]') - do ii = 1, num_fluids - myRho = myRho + myalpha_rho(ii) - n_tait = n_tait + myalpha(ii)*gammas(ii) - B_tait = B_tait + myalpha(ii)*pi_infs(ii)/pi_fac - end do + if (num_fluids == 1) then + myRho = myalpha_rho(1) + n_tait = gammas(1) + B_tait = pi_infs(1)/pi_fac + else + myRho = 0._wp + n_tait = 0._wp + B_tait = 0._wp + + $:GPU_LOOP(parallelism='[seq]') + do ii = 1, num_fluids + myRho = myRho + myalpha_rho(ii) + n_tait = n_tait + myalpha(ii)*gammas(ii) + B_tait = B_tait + myalpha(ii)*pi_infs(ii)/pi_fac + end do + end if n_tait = 1._wp/n_tait + 1._wp !make this the usual little 'gamma' B_tait = B_tait*(n_tait - 1)/n_tait ! make this the usual pi_inf From 50f78b8ad187feba85170de8cc4cdc45e9242414 Mon Sep 17 00:00:00 2001 From: Hyeoksu Lee Date: Sat, 29 Nov 2025 11:41:46 -0800 Subject: [PATCH 4/7] correction for vf of EE bubble --- src/common/m_variables_conversion.fpp | 47 +++++++++++++++++++-------- src/simulation/m_sim_helpers.fpp | 2 -- 2 files changed, 33 insertions(+), 16 deletions(-) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 582dc89ce4..398a9fa023 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -264,13 +264,20 @@ contains ! Calculating the density, the specific heat ratio function, the ! liquid stiffness function, and the energy reference function, ! respectively, from the species analogs - rho = 0._wp; gamma = 0._wp; pi_inf = 0._wp; qv = 0._wp - do i = 1, num_fluids - rho = rho + alpha_rho_K(i) - gamma = gamma + alpha_K(i)*gammas(i) - pi_inf = pi_inf + alpha_K(i)*pi_infs(i) - qv = qv + alpha_rho_K(i)*qvs(i) - end do + if (num_fluids == 1 .and. bubbles_euler) then + rho = alpha_rho_K(1) + gamma = gammas(1) + pi_inf = pi_infs(1) + qv = qvs(1) + else + rho = 0._wp; gamma = 0._wp; pi_inf = 0._wp; qv = 0._wp + do i = 1, num_fluids + rho = rho + alpha_rho_K(i) + gamma = gamma + alpha_K(i)*gammas(i) + pi_inf = pi_inf + alpha_K(i)*pi_infs(i) + qv = qv + alpha_rho_K(i)*qvs(i) + end do + end if #ifdef MFC_SIMULATION ! Computing the shear and bulk Reynolds numbers from species analogs @@ -326,13 +333,20 @@ contains ! their physical bounds to make sure that any mixture variables that ! are derived from them result within the limits that are set by the ! fluids physical parameters that make up the mixture - rho_K = 0._wp; gamma_K = 0._wp; pi_inf_K = 0._wp; qv_K = 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) - end do + 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 + 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 if (present(G_K)) then G_K = 0._wp @@ -595,6 +609,8 @@ contains call s_compute_species_fraction(qK_cons_vf, j, k, l, alpha_rho_K, alpha_K) + ! print *, j, k, l, alpha_rho_K, qK_cons_vf(1)%sf(j, k, l), alpha_K, qK_cons_vf(advxb)%sf(j, k, l) + if (model_eqns /= 4) then #ifdef MFC_SIMULATION ! If in simulation, use acc mixture subroutines @@ -1324,9 +1340,12 @@ contains alpha_rho_K(i) = max(0._wp, alpha_rho_K(i)) alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp) end do + alpha_K = alpha_K/max(sum(alpha_K), 1.e-16_wp) end if + if (num_fluids == 1 .and. bubbles_euler) alpha_K(1) = q_vf(advxb)%sf(k, l, r) + end subroutine s_compute_species_fraction impure subroutine s_finalize_variables_conversion_module() diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index d4c7e5a87c..0aa632d3ef 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -118,8 +118,6 @@ contains call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re) end if - if (num_fluids == 1 .and. bubbles_euler) alpha = q_prim_vf(alf_idx)%sf(j, k, l) - if (igr) then $:GPU_LOOP(parallelism='[seq]') do i = 1, num_vels From 1568228bcf805d030cf6a99a544e3fe66924e1d0 Mon Sep 17 00:00:00 2001 From: Hyeoksu Lee Date: Sat, 29 Nov 2025 15:54:00 -0800 Subject: [PATCH 5/7] add mpp_lim to mixture rule acc --- src/common/m_variables_conversion.fpp | 10 +++++++--- src/simulation/m_bubbles_EL.fpp | 6 +----- src/simulation/m_hyperelastic.fpp | 8 +++----- 3 files changed, 11 insertions(+), 13 deletions(-) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 398a9fa023..3f4e487c0b 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -339,6 +339,13 @@ contains pi_inf_K = pi_infs(1) qv_K = qvs(1) else + if (mpp_lim) then + 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) + end do + alpha_K = alpha_K/max(sum(alpha_K), sgm_eps) + end if rho_K = 0._wp; gamma_K = 0._wp; pi_inf_K = 0._wp; qv_K = 0._wp do i = 1, num_fluids rho_K = rho_K + alpha_rho_K(i) @@ -609,8 +616,6 @@ contains call s_compute_species_fraction(qK_cons_vf, j, k, l, alpha_rho_K, alpha_K) - ! print *, j, k, l, alpha_rho_K, qK_cons_vf(1)%sf(j, k, l), alpha_K, qK_cons_vf(advxb)%sf(j, k, l) - if (model_eqns /= 4) then #ifdef MFC_SIMULATION ! If in simulation, use acc mixture subroutines @@ -1340,7 +1345,6 @@ contains alpha_rho_K(i) = max(0._wp, alpha_rho_K(i)) alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp) end do - alpha_K = alpha_K/max(sum(alpha_K), 1.e-16_wp) end if diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index 9b42dc207b..47e0709788 100644 --- a/src/simulation/m_bubbles_EL.fpp +++ b/src/simulation/m_bubbles_EL.fpp @@ -642,11 +642,7 @@ contains call s_get_pinf(k, q_prim_vf, 1, myPinf, cell, aux1, aux2) ! Obtain liquid density and computing speed of sound from pinf - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - myalpha_rho(i) = q_prim_vf(i)%sf(cell(1), cell(2), cell(3)) - myalpha(i) = q_prim_vf(E_idx + i)%sf(cell(1), cell(2), cell(3)) - end do + call s_compute_species_fraction(q_prim_vf, cell(1), cell(2), cell(3), myalpha_rho, myalpha) call s_convert_species_to_mixture_variables_acc(myRho, gamma, pi_inf, qv, myalpha, & myalpha_rho, Re) call s_compute_cson_from_pinf(q_prim_vf, myPinf, cell, myRho, gamma, pi_inf, myCson) diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index 9e21ef8820..21625069b2 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -110,11 +110,9 @@ contains do l = 0, p do k = 0, n do j = 0, m - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l) - alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l) - end do + + call s_compute_species_fraction(q_cons_vf, j, k, l, alpha_rho_k, alpha_k) + ! If in simulation, use acc mixture subroutines call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha_k, & alpha_rho_k, Re, G_local, Gs_hyper) From 4b7ab1c884a27f86fe1f945b9e06040459d0cb94 Mon Sep 17 00:00:00 2001 From: Hyeoksu Lee Date: Sat, 29 Nov 2025 16:09:10 -0800 Subject: [PATCH 6/7] add comments --- src/common/m_variables_conversion.fpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 3f4e487c0b..0f5654f3e5 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -1309,6 +1309,7 @@ contains #endif end subroutine s_convert_primitive_to_flux_variables + !> This subroutine computes partial densities and volume fractions subroutine s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K) $:GPU_ROUTINE(function_name='s_compute_species_fraction', & & parallelism='[seq]', cray_inline=True) From f882347e90c14dccd1fa4199ad46f02e8aeb2634 Mon Sep 17 00:00:00 2001 From: Hyeoksu Lee Date: Sat, 29 Nov 2025 22:23:48 -0800 Subject: [PATCH 7/7] update bubble case validation --- toolchain/mfc/case_validator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 4eea0e2940..085754abdc 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -369,8 +369,8 @@ def check_stiffened_eos(self): if num_fluids is None: return - # Allow one extra fluid property slot when using bubbles_euler with single-component flows - bub_fac = 1 if (bubbles_euler and num_fluids == 1) else 0 + # Allow one extra fluid property slot when using bubbles_euler + bub_fac = 1 if (bubbles_euler) else 0 for i in range(1, num_fluids + 1 + bub_fac): gamma = self.get(f'fluid_pp({i})%gamma')