Skip to content

fix(sim): Euler–Euler bubbles GPU data race on shared bubble-wall scratch #1648

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
sbryngelson merged 2 commits into from
Jul 17, 2026
+18 −19
Merged

fix(sim): Euler–Euler bubbles GPU data race on shared bubble-wall scratch #1648

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
288a3d7
bubbles: fix Euler-Euler GPU data race on shared bubble-wall scratch
sbryngelson Jul 11, 2026
3d2551c
style: format m_bubbles_EE.fpp (ffmt line wrap)
sbryngelson Jul 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 7 additions & 10 deletions src/simulation/m_bubbles.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Comment on lines 260 to +288

@sbryngelson sbryngelson Jul 16, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The thermal == 3 .and. .not. bubbles_lagrange branch that reads fchi_vw/frho_mw is only reachable from the Euler–Euler source loop (m_bubbles_EE), which always computes these via s_bwproperty into per-thread scratch and passes them (fchi_vw=chi_vw_l, frho_mw=rho_mw_l). The other callers — the Lagrangian path (m_bubbles_EL) and the adaptive substep — run with bubbles_lagrange = .true., so they take the bubbles_lagrange branch above and never touch these optionals. So the branch is never entered with the optionals absent.

These are $:GPU_ROUTINE(parallelism='[seq]') device routines, so an @:ASSERT(present(...)) guard isn't available on-device, and a local s_bwproperty fallback would duplicate the bubble-wall computation into the device path for an unreachable case; the coupling is enforced by the caller structure instead. (PR is now rebased onto master post-#1290, which reworked these routines to non-elemental — the fix re-applies cleanly.)

end if
else
! polytropic
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Comment on lines 307 to 325

@sbryngelson sbryngelson Jul 16, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same reasoning as the s_vflux thread: the thermal == 3 non-Lagrangian path that uses fk_mw is only reached from the Euler–Euler loop (m_bubbles_EE), which passes fk_mw=k_mw_l; the Lagrangian/adaptive callers take the bubbles_lagrange branch (which returns before fk_mw is used). Device $:GPU_ROUTINE, so no on-device assert; coupling enforced by the callers.

f_bpres_dot = -3._wp*gam_m*fV/fR*(fpb - pv)
end if
Expand Down
20 changes: 11 additions & 9 deletions src/simulation/m_bubbles_EE.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading