From b20291f46ebb01ae92756f30779e47e2b9d443ab Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:19:50 -0500 Subject: [PATCH 1/4] Fix MUSCL THINC right-state using already-overwritten left-state values The right reconstruction divides by left-state values that were already overwritten during left reconstruction. Save original density ratios (contxb/advxb and contxe/(1-advxb)) before left reconstruction and reuse them for both left and right states. Co-Authored-By: Claude Opus 4.6 --- src/simulation/m_muscl.fpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/simulation/m_muscl.fpp b/src/simulation/m_muscl.fpp index de1af10f24..f0b1fa9e14 100644 --- a/src/simulation/m_muscl.fpp +++ b/src/simulation/m_muscl.fpp @@ -248,11 +248,12 @@ contains integer :: j, k, l real(wp) :: aCL, aCR, aC, aTHINC, qmin, qmax, A, B, C, sign, moncon + real(wp) :: rho_b, rho_e #:for MUSCL_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')] if (muscl_dir == ${MUSCL_DIR}$) then - $:GPU_PARALLEL_LOOP(collapse=3,private='[j,k,l,aCL,aC,aCR,aTHINC,moncon,sign,qmin,qmax]') + $:GPU_PARALLEL_LOOP(collapse=3,private='[j,k,l,aCL,aC,aCR,aTHINC,moncon,sign,qmin,qmax,rho_b,rho_e]') do l = is3_muscl%beg, is3_muscl%end do k = is2_muscl%beg, is2_muscl%end do j = is1_muscl%beg, is1_muscl%end @@ -278,14 +279,16 @@ contains B = exp(sign*ic_beta*(2._wp*C - 1._wp)) A = (B/cosh(ic_beta) - 1._wp)/tanh(ic_beta) + ! Save original density ratios before THINC overwrites them + rho_b = vL_rs_vf_${XYZ}$ (j, k, l, contxb)/vL_rs_vf_${XYZ}$ (j, k, l, advxb) + rho_e = vL_rs_vf_${XYZ}$ (j, k, l, contxe)/(1._wp - vL_rs_vf_${XYZ}$ (j, k, l, advxb)) + ! Left reconstruction aTHINC = qmin + 5e-1_wp*qmax*(1._wp + sign*A) if (aTHINC < ic_eps) aTHINC = ic_eps if (aTHINC > 1 - ic_eps) aTHINC = 1 - ic_eps - vL_rs_vf_${XYZ}$ (j, k, l, contxb) = vL_rs_vf_${XYZ}$ (j, k, l, contxb)/ & - vL_rs_vf_${XYZ}$ (j, k, l, advxb)*aTHINC - vL_rs_vf_${XYZ}$ (j, k, l, contxe) = vL_rs_vf_${XYZ}$ (j, k, l, contxe)/ & - (1._wp - vL_rs_vf_${XYZ}$ (j, k, l, advxb))*(1._wp - aTHINC) + vL_rs_vf_${XYZ}$ (j, k, l, contxb) = rho_b*aTHINC + vL_rs_vf_${XYZ}$ (j, k, l, contxe) = rho_e*(1._wp - aTHINC) vL_rs_vf_${XYZ}$ (j, k, l, advxb) = aTHINC vL_rs_vf_${XYZ}$ (j, k, l, advxe) = 1 - aTHINC @@ -293,10 +296,8 @@ contains aTHINC = qmin + 5e-1_wp*qmax*(1._wp + sign*(tanh(ic_beta) + A)/(1._wp + A*tanh(ic_beta))) if (aTHINC < ic_eps) aTHINC = ic_eps if (aTHINC > 1 - ic_eps) aTHINC = 1 - ic_eps - vR_rs_vf_${XYZ}$ (j, k, l, contxb) = vL_rs_vf_${XYZ}$ (j, k, l, contxb)/ & - vL_rs_vf_${XYZ}$ (j, k, l, advxb)*aTHINC - vR_rs_vf_${XYZ}$ (j, k, l, contxe) = vL_rs_vf_${XYZ}$ (j, k, l, contxe)/ & - (1._wp - vL_rs_vf_${XYZ}$ (j, k, l, advxb))*(1._wp - aTHINC) + vR_rs_vf_${XYZ}$ (j, k, l, contxb) = rho_b*aTHINC + vR_rs_vf_${XYZ}$ (j, k, l, contxe) = rho_e*(1._wp - aTHINC) vR_rs_vf_${XYZ}$ (j, k, l, advxb) = aTHINC vR_rs_vf_${XYZ}$ (j, k, l, advxe) = 1 - aTHINC From e7c7ea8977174b8abcb966fd66de28d119f3ee38 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Mon, 23 Feb 2026 09:38:33 -0500 Subject: [PATCH 2/4] Add missing A, B, C to THINC GPU private clause These scalars are written per-thread inside the parallel loop (lines 278-280) but were missing from the private list, causing a GPU data race on all backends. Co-Authored-By: Claude Opus 4.6 --- src/simulation/m_muscl.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_muscl.fpp b/src/simulation/m_muscl.fpp index f0b1fa9e14..0db1063246 100644 --- a/src/simulation/m_muscl.fpp +++ b/src/simulation/m_muscl.fpp @@ -253,7 +253,7 @@ contains #:for MUSCL_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')] if (muscl_dir == ${MUSCL_DIR}$) then - $:GPU_PARALLEL_LOOP(collapse=3,private='[j,k,l,aCL,aC,aCR,aTHINC,moncon,sign,qmin,qmax,rho_b,rho_e]') + $:GPU_PARALLEL_LOOP(collapse=3,private='[j,k,l,aCL,aC,aCR,aTHINC,moncon,sign,qmin,qmax,A,B,C,rho_b,rho_e]') do l = is3_muscl%beg, is3_muscl%end do k = is2_muscl%beg, is2_muscl%end do j = is1_muscl%beg, is1_muscl%end From c689ec7166792a910fe99c346a785301107af7f0 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 24 Feb 2026 18:22:10 -0500 Subject: [PATCH 3/4] Guard THINC density ratio computation against divide-by-zero When advxb approaches 0 or 1 at an interface, rho_b and rho_e could divide by zero. Guard with sgm_eps (1e-16) using the same pattern already used for C = (aC - qmin + sgm_eps)/(qmax + sgm_eps). Co-Authored-By: Claude Sonnet 4.6 --- src/simulation/m_muscl.fpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/simulation/m_muscl.fpp b/src/simulation/m_muscl.fpp index 0db1063246..8518b9a394 100644 --- a/src/simulation/m_muscl.fpp +++ b/src/simulation/m_muscl.fpp @@ -280,8 +280,10 @@ contains A = (B/cosh(ic_beta) - 1._wp)/tanh(ic_beta) ! Save original density ratios before THINC overwrites them - rho_b = vL_rs_vf_${XYZ}$ (j, k, l, contxb)/vL_rs_vf_${XYZ}$ (j, k, l, advxb) - rho_e = vL_rs_vf_${XYZ}$ (j, k, l, contxe)/(1._wp - vL_rs_vf_${XYZ}$ (j, k, l, advxb)) + rho_b = vL_rs_vf_${XYZ}$ (j, k, l, contxb)/ & + max(vL_rs_vf_${XYZ}$ (j, k, l, advxb), sgm_eps) + rho_e = vL_rs_vf_${XYZ}$ (j, k, l, contxe)/ & + max(1._wp - vL_rs_vf_${XYZ}$ (j, k, l, advxb), sgm_eps) ! Left reconstruction aTHINC = qmin + 5e-1_wp*qmax*(1._wp + sign*A) From 5766cec2b87d5d93054f417fe0dd1e23b23e8b16 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 24 Feb 2026 19:44:08 -0500 Subject: [PATCH 4/4] Revert "Guard THINC density ratio computation against divide-by-zero" This reverts commit c689ec71. Co-Authored-By: Claude Sonnet 4.6 --- src/simulation/m_muscl.fpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/simulation/m_muscl.fpp b/src/simulation/m_muscl.fpp index 8518b9a394..0db1063246 100644 --- a/src/simulation/m_muscl.fpp +++ b/src/simulation/m_muscl.fpp @@ -280,10 +280,8 @@ contains A = (B/cosh(ic_beta) - 1._wp)/tanh(ic_beta) ! Save original density ratios before THINC overwrites them - rho_b = vL_rs_vf_${XYZ}$ (j, k, l, contxb)/ & - max(vL_rs_vf_${XYZ}$ (j, k, l, advxb), sgm_eps) - rho_e = vL_rs_vf_${XYZ}$ (j, k, l, contxe)/ & - max(1._wp - vL_rs_vf_${XYZ}$ (j, k, l, advxb), sgm_eps) + rho_b = vL_rs_vf_${XYZ}$ (j, k, l, contxb)/vL_rs_vf_${XYZ}$ (j, k, l, advxb) + rho_e = vL_rs_vf_${XYZ}$ (j, k, l, contxe)/(1._wp - vL_rs_vf_${XYZ}$ (j, k, l, advxb)) ! Left reconstruction aTHINC = qmin + 5e-1_wp*qmax*(1._wp + sign*A)