diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp index 8a0eb0d20d..639c2aa27b 100644 --- a/src/simulation/m_bubbles.fpp +++ b/src/simulation/m_bubbles.fpp @@ -17,11 +17,6 @@ module m_bubbles implicit none - real(wp) :: chi_vw !< Bubble wall properties (Ando 2010) - real(wp) :: k_mw !< Bubble wall properties (Ando 2010) - real(wp) :: rho_mw !< Bubble wall properties (Ando 2010) - $:GPU_DECLARE(create='[chi_vw, k_mw, rho_mw]') - contains !> Compute the bubble radial acceleration based on the selected bubble model @@ -251,7 +246,7 @@ contains end subroutine s_bwproperty !> Compute the vapour flux - subroutine s_vflux(fR, fV, fpb, fmass_v, iR0, vflux, fmass_g, fbeta_c, fR_m, fgamma_m) + subroutine s_vflux(fR, fV, fpb, fmass_v, iR0, vflux, fmass_g, fbeta_c, fR_m, fgamma_m, fchi_vw, frho_mw) $:GPU_ROUTINE(parallelism='[seq]') real(wp), intent(in) :: fR @@ -262,6 +257,7 @@ contains real(wp), intent(out) :: vflux real(wp), intent(in), optional :: fmass_g, fbeta_c real(wp), intent(out), optional :: fR_m, fgamma_m + real(wp), intent(in), optional :: fchi_vw, frho_mw !< Per-thread bubble-wall scratch (Euler-Euler path) real(wp) :: chi_bar real(wp) :: rho_mw_lag real(wp) :: grad_chi @@ -288,8 +284,8 @@ contains end if else chi_bar = fmass_v/(fmass_v + mass_g0(iR0)) - grad_chi = -Re_trans_c(iR0)*(chi_bar - chi_vw) - vflux = rho_mw*grad_chi/Pe_c/(1._wp - chi_vw)/fR + grad_chi = -Re_trans_c(iR0)*(chi_bar - fchi_vw) + vflux = frho_mw*grad_chi/Pe_c/(1._wp - fchi_vw)/fR end if else ! polytropic @@ -299,7 +295,7 @@ contains end subroutine s_vflux !> Compute the time derivative of the internal bubble pressure - function f_bpres_dot(fvflux, fR, fV, fpb, fmass_v, iR0, fbeta_t, fR_m, fgamma_m) + function f_bpres_dot(fvflux, fR, fV, fpb, fmass_v, iR0, fbeta_t, fR_m, fgamma_m, fk_mw) $:GPU_ROUTINE(parallelism='[seq]') real(wp), intent(in) :: fvflux @@ -309,6 +305,7 @@ contains real(wp), intent(in) :: fmass_v integer, intent(in) :: iR0 real(wp), intent(in), optional :: fbeta_t, fR_m, fgamma_m + real(wp), intent(in), optional :: fk_mw !< Per-thread bubble-wall conductivity scratch (Euler-Euler path) real(wp) :: T_bar real(wp) :: grad_T real(wp) :: f_bpres_dot @@ -324,7 +321,7 @@ contains end if grad_T = -Re_trans_T(iR0)*((fpb/pb0(iR0))*(fR/R0(iR0))**3*(mass_g0(iR0) + mass_v0(iR0))/(mass_g0(iR0) + fmass_v) & & - 1._wp) - f_bpres_dot = 3._wp*gam_m*(-fV*fpb + fvflux*R_v*Tw + pb0(iR0)*k_mw*grad_T/Pe_T(iR0)/fR)/fR + f_bpres_dot = 3._wp*gam_m*(-fV*fpb + fvflux*R_v*Tw + pb0(iR0)*fk_mw*grad_T/Pe_T(iR0)/fR)/fR else f_bpres_dot = -3._wp*gam_m*fV/fR*(fpb - pv) end if diff --git a/src/simulation/m_bubbles_EE.fpp b/src/simulation/m_bubbles_EE.fpp index b07bf35f07..9e3f2f29ec 100644 --- a/src/simulation/m_bubbles_EE.fpp +++ b/src/simulation/m_bubbles_EE.fpp @@ -141,12 +141,13 @@ contains impure subroutine s_compute_bubble_EE_source(q_cons_vf, q_prim_vf, rhs_vf, divu_in) type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf - type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf + type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf - type(scalar_field), intent(in) :: divu_in !< matrix for div(u) - real(wp) :: rddot - real(wp) :: pb_local, mv_local, vflux, pbdot - real(wp) :: n_tait, B_tait + type(scalar_field), intent(in) :: divu_in !< matrix for div(u) + real(wp) :: rddot + real(wp) :: pb_local, mv_local, vflux, pbdot + real(wp) :: n_tait, B_tait + real(wp) :: chi_vw_l, k_mw_l, rho_mw_l !< Per-thread bubble-wall scratch (avoid module-scalar race) #:if not MFC_CASE_OPTIMIZATION and USING_AMD real(wp), dimension(3) :: Rtmp, Vtmp @@ -182,7 +183,8 @@ contains adap_dt_stop_sum = 0 $:GPU_PARALLEL_LOOP(private='[j, k, l, Rtmp, Vtmp, myalpha_rho, myalpha, myR, myV, alf, myP, myRho, R2Vav, R3, nbub, & - & pb_local, mv_local, vflux, pbdot, rddot, n_tait, B_tait, adap_dt_stop]', collapse=3, copy='[adap_dt_stop_sum]') + & pb_local, mv_local, vflux, pbdot, rddot, n_tait, B_tait, adap_dt_stop, chi_vw_l, k_mw_l, & + & rho_mw_l]', collapse=3, copy='[adap_dt_stop_sum]') do l = 0, p do k = 0, n do j = 0, m @@ -261,9 +263,9 @@ contains if (.not. polytropic) then pb_local = q_prim_vf(ps(q))%sf(j, k, l) mv_local = q_prim_vf(ms(q))%sf(j, k, l) - call s_bwproperty(pb_local, q, chi_vw, k_mw, rho_mw) - call s_vflux(myR, myV, pb_local, mv_local, q, vflux) - pbdot = f_bpres_dot(vflux, myR, myV, pb_local, mv_local, q) + call s_bwproperty(pb_local, q, chi_vw_l, k_mw_l, rho_mw_l) + call s_vflux(myR, myV, pb_local, mv_local, q, vflux, fchi_vw=chi_vw_l, frho_mw=rho_mw_l) + pbdot = f_bpres_dot(vflux, myR, myV, pb_local, mv_local, q, fk_mw=k_mw_l) bub_p_src(j, k, l, q) = nbub*pbdot bub_m_src(j, k, l, q) = nbub*vflux*4._wp*pi*(myR**2._wp) else