diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 307406e6b..0f5654f3e 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -33,14 +33,13 @@ 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, & 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, & @@ -93,11 +92,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 +225,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 @@ -379,72 +251,48 @@ 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 + call s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K) - 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 + ! Calculating the density, the specific heat ratio function, the + ! liquid stiffness function, and the energy reference function, + ! respectively, from the species analogs + 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 - 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) + 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 - - alpha_K = alpha_K/max(sum(alpha_K), 1.e-16_wp) - end if - ! 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 #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) + Re_K(i) = 1._wp/max(Re_K(i), sgm_eps) end do - - Re_K(i) = 1._wp/max(Re_K(i), sgm_eps) - - end do + end if #endif if (present(G_K)) then @@ -473,47 +321,40 @@ 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 + 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 + 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 - 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) + 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 - - 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 (present(G_K)) then G_K = 0._wp do i = 1, num_fluids @@ -525,7 +366,6 @@ contains end if if (viscous) then - do i = 1, 2 Re_K(i) = dflt_real @@ -537,77 +377,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. @@ -617,7 +392,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)) @@ -626,16 +400,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 @@ -664,12 +428,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 @@ -855,27 +614,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 @@ -883,9 +622,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) @@ -1479,6 +1215,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) @@ -1495,9 +1232,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) @@ -1575,6 +1309,50 @@ 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) + 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 + + 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() ! 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 d8bed3ad1..9aa564590 100644 --- a/src/simulation/m_bubbles_EE.fpp +++ b/src/simulation/m_bubbles_EE.fpp @@ -242,34 +242,26 @@ 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 + 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 - 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) + B_tait = B_tait + myalpha(ii)*pi_infs(ii)/pi_fac end do - else - myRho = myalpha_rho(1) - n_tait = gammas(1) - B_tait = pi_infs(1)/pi_fac 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 - 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_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index 9b42dc207..47e070978 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_cbc.fpp b/src/simulation/m_cbc.fpp index 93bf4ae64..b5027b9fd 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_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index 9e21ef882..21625069b 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) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 1f3c9a623..aff2a1882 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 4e42768b2..0aa632d3e 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -109,33 +109,11 @@ 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, & 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 diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 4eea0e294..085754abd 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')