Pass scalars into the EOS device routines: Cray OpenACC misaddresses a device-array element passed by reference

Author: sbryngelson

.claude/rules/common-pitfalls.md

@@ -60,6 +60,15 @@ covered in `docs/documentation/contributing.md`.
   QBMM/viscous and MHD HLLD, while both Lagrange bubble cases *complete* with out-of-tolerance
   answers. Measured 2026-08-29 on MI210. A compile-only check returns green, so any future attempt to
   drop these must run the tests, not just build.
+- Pass **scalars** into device routines, never an element of a device-resident array: a `declare create`
+  module array (`blkmod1(k,l,q)`, `gammas(i)`), an attached field (`q_prim_vf%vf(i)%sf(j,k,l)`), or a
+  dummy sized by a device global (`alpha_K(i)` with `dimension(num_fluids)`). Copy the element to a
+  local first and receive outputs into a local. CCE OpenACC compiles it clean and, unless it inlines
+  the callee, reads and writes the element at the wrong address: on PR #1811 a kernel wrote
+  `blkmod1(:,:,:)` through such a call and every cell stayed 0 while the host got 1.4 from the same
+  fields; every test through the path ended in `NaN(s) in timestep output` while CCE OpenMP, amdflang
+  and nvfortran were bit-identical. A leaf callee that gets inlined hides the bug, which is why the
+  parent PR passed. Expressions (`a/max(b, eps)`) are temporaries and are fine.
 - The same "call it from the loop body" rule covers `m_thermochem`: calling `get_species_*` from
   inside a `GPU_ROUTINE` rather than from the kernel gave CCE OpenMP a runtime
   `Memory access fault by GPU node-N ... Reason: Unknown` on the first step (exit 134), while every

src/common/m_variables_conversion.fpp

@@ -752,6 +752,7 @@ contains
         real(wp)                         :: rho
         real(wp)                         :: gamma
         real(wp)                         :: pi_inf
+        real(wp)                         :: pres_i, alpha_i, alpha_rho_i, e_i
         real(wp)                         :: qv
         real(wp)                         :: dyn_pres
         real(wp)                         :: nbub, R3tmp
@@ -894,10 +895,11 @@ contains
                     ! Six-equation model (Saurel et al. JCP 2009): compute per-phase internal energies
                     if (model_eqns == model_eqns_6eq) then
                         do i = 1, num_fluids
-                            call s_phase_internal_energy(q_prim_vf(eqn_idx%E)%sf(j, k, l), &
-                                                         & q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l), &
-                                                         & q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l), i, &
-                                                         & q_cons_vf(i + eqn_idx%int_en%beg - 1)%sf(j, k, l))
+                            pres_i = q_prim_vf(eqn_idx%E)%sf(j, k, l)
+                            alpha_i = q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l)
+                            alpha_rho_i = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)
+                            call s_phase_internal_energy(pres_i, alpha_i, alpha_rho_i, i, e_i)
+                            q_cons_vf(i + eqn_idx%int_en%beg - 1)%sf(j, k, l) = e_i
                         end do
                     end if
 
@@ -1241,7 +1243,7 @@ contains
         #:endif
         real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K
         real(wp)              :: gamma_i, pi_inf_i, dpi_i, dgamma_i
-        real(wp)              :: rho_i
+        real(wp)              :: rho_i, alpha_i, alpha_rho_i
         integer               :: i  !< Loop iterator over fluids
 
         ! The bubbly closure is written for one carrier liquid, which keeps its own coefficients
@@ -1263,7 +1265,9 @@ contains
             $:GPU_LOOP(parallelism='[seq]')
             do i = 1, num_fluids
                 rho_K = rho_K + alpha_rho_K(i)
-                call s_phase_coefficients(alpha_rho_K(i), alpha_K(i), i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i)
+                alpha_rho_i = alpha_rho_K(i)
+                alpha_i = alpha_K(i)
+                call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i)
                 gamma_K = gamma_K + alpha_K(i)*gamma_i
                 pi_inf_K = pi_inf_K + alpha_K(i)*pi_inf_i
                 qv_K = qv_K + alpha_rho_K(i)*qvs(i)
@@ -1283,7 +1287,7 @@ contains
             real(wp), dimension(num_fluids), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv
         #:endif
         real(wp), intent(out) :: drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt
-        real(wp)              :: rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i
+        real(wp)              :: rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i, alpha_i, alpha_rho_i
         integer               :: i  !< Loop iterator over fluids
 
         dgamma_dt = 0._wp
@@ -1299,7 +1303,9 @@ contains
             $:GPU_LOOP(parallelism='[seq]')
             do i = 1, num_fluids
                 drho_dt = drho_dt + dalpha_rho_dt(i)
-                call s_phase_coefficients(alpha_rho(i), adv(i), i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i)
+                alpha_rho_i = alpha_rho(i)
+                alpha_i = adv(i)
+                call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i)
                 ! d(alpha X(rho_i))/dt with rho_i = alpha_rho/alpha; the alpha in dX/dt cancels
                 dgamma_dt = dgamma_dt + dadv_dt(i)*gamma_i + dgamma_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i))
                 dpi_inf_dt = dpi_inf_dt + dadv_dt(i)*pi_inf_i + dpi_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i))
@@ -1599,11 +1605,12 @@ contains
         real(wp), intent(in)  :: rho, pres
         integer, intent(in)   :: i
         real(wp), intent(out) :: T
-        real(wp)              :: p_ref, e_ref, dp_drho, de_drho, G0, dG0, T_ref
+        real(wp)              :: p_ref, e_ref, dp_drho, de_drho, G0, dG0, T0, T_ref
 
         if (f_is_state_dependent(i)) then
             call s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, G0, dG0)
-            call s_rk4(ode_reference_temperature, i, 1._wp/rho0s(i), t0s(i), 1._wp/rho, T_ref)
+            T0 = t0s(i)
+            call s_rk4(ode_reference_temperature, i, 1._wp/rho0s(i), T0, 1._wp/rho, T_ref)
             T = T_ref + (pres - p_ref)/(rho*G0*cvs(i))
         else
             T = (pres + isentrope_B(i))/((isentrope_n(i) - 1._wp)*cvs(i)*rho)
@@ -1758,7 +1765,7 @@ contains
             real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
         #:endif
         real(wp) :: alf  !< Subgrid void fraction; dilute by construction
-        real(wp) :: blkmod_q
+        real(wp) :: blkmod_q, alpha_q, alpha_rho_q, gamma_q, pi_inf_q
         integer  :: q
 
         if (chemistry) then  ! Reacting mixture sound speed
@@ -1772,7 +1779,9 @@ contains
                 c = 0._wp
                 $:GPU_LOOP(parallelism='[seq]')
                 do q = 1, num_fluids
-                    call s_phase_bulk_modulus(pres, adv(q), alpha_rho(q), q, blkmod_q)
+                    alpha_q = adv(q)
+                    alpha_rho_q = alpha_rho(q)
+                    call s_phase_bulk_modulus(pres, alpha_q, alpha_rho_q, q, blkmod_q)
                     if (alt_soundspeed) then
                         c = c + adv(q)/blkmod_q
                     else
@@ -1788,14 +1797,18 @@ contains
                 c = 0._wp
                 $:GPU_LOOP(parallelism='[seq]')
                 do q = 1, num_fluids
-                    c = c + adv(q)/f_bulk_modulus(pres, gammas(q), pi_infs(q))
+                    gamma_q = gammas(q)
+                    pi_inf_q = pi_infs(q)
+                    c = c + adv(q)/f_bulk_modulus(pres, gamma_q, pi_inf_q)
                 end do
                 c = 1._wp/(rho*c)
             else if (model_eqns == model_eqns_6eq) then  ! volume-weighted arithmetic mean
                 c = 0._wp
                 $:GPU_LOOP(parallelism='[seq]')
                 do q = 1, num_fluids
-                    c = c + adv(q)*f_bulk_modulus(pres, gammas(q), pi_infs(q))
+                    gamma_q = gammas(q)
+                    pi_inf_q = pi_infs(q)
+                    c = c + adv(q)*f_bulk_modulus(pres, gamma_q, pi_inf_q)
                 end do
                 c = c/rho
             else  ! the mixture coefficients already carry the mixing

src/simulation/m_acoustic_src.fpp

@@ -136,7 +136,7 @@ contains
         #:endif
         real(wp)                            :: myRho, pi_inf_mix, qv_dummy
         real(wp)                            :: sim_time, c, gamma_mix
-        real(wp)                            :: blkmod_q
+        real(wp)                            :: blkmod_q, pres_q, alpha_q, alpha_rho_q
         real(wp)                            :: frequency_local, gauss_sigma_time_local
         real(wp)                            :: mass_src_diff, mom_src_diff
         real(wp)                            :: source_temporal
@@ -206,9 +206,10 @@ contains
 
                 deallocate (phi_rn)
 
-                $:GPU_PARALLEL_LOOP(private='[myalpha, myalpha_rho, myRho, pi_inf_mix, qv_dummy, c, blkmod_q, gamma_mix, &
-                                    & frequency_local, gauss_sigma_time_local, mass_src_diff, mom_src_diff, source_temporal, j, &
-                                    & k, l, q]', copyin = '[sum_BB, freq_conv_flag, gauss_conv_flag, sim_time]')
+                $:GPU_PARALLEL_LOOP(private='[myalpha, myalpha_rho, myRho, pi_inf_mix, qv_dummy, c, blkmod_q, pres_q, alpha_q, &
+                                    & alpha_rho_q, gamma_mix, frequency_local, gauss_sigma_time_local, mass_src_diff, &
+                                    & mom_src_diff, source_temporal, j, k, l, q]', copyin = '[sum_BB, freq_conv_flag, &
+                                    & gauss_conv_flag, sim_time]')
                 do i = 1, num_points
                     j = source_spatials(ai)%coord(1, i)
                     k = source_spatials(ai)%coord(2, i)
@@ -229,7 +230,10 @@ contains
                         c = 0._wp
                         $:GPU_LOOP(parallelism='[seq]')
                         do q = 1, num_fluids
-                            call s_phase_bulk_modulus(q_prim_vf(eqn_idx%E)%sf(j, k, l), myalpha(q), myalpha_rho(q), q, blkmod_q)
+                            pres_q = q_prim_vf(eqn_idx%E)%sf(j, k, l)
+                            alpha_q = myalpha(q)
+                            alpha_rho_q = myalpha_rho(q)
+                            call s_phase_bulk_modulus(pres_q, alpha_q, alpha_rho_q, q, blkmod_q)
                             c = c + myalpha(q)*blkmod_q
                         end do
                         c = c/myRho

src/simulation/m_hypoelastic.fpp

@@ -599,14 +599,14 @@ contains
     !! @param nc_iface_vel_y_hatR_vf hat_R-pass radial-direction interface velocities
     subroutine s_compute_hypoelastic_rhs_axisym_geom_dual_pass(q_prim_vf, rhs_vf, nc_iface_vel_y_vf, nc_iface_vel_y_hatR_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), dimension(:), intent(in)           :: nc_iface_vel_y_vf
-        type(scalar_field), dimension(:), intent(in)           :: nc_iface_vel_y_hatR_vf
-        real(wp)                                               :: rho_K, G_K, K_K, C_num, pres_K, blkmod1_K, blkmod2_K
-        integer                                                :: i, k, l, q
+        type(scalar_field), dimension(:), intent(in) :: nc_iface_vel_y_vf
+        type(scalar_field), dimension(:), intent(in) :: nc_iface_vel_y_hatR_vf
+        real(wp) :: rho_K, G_K, K_K, C_num, pres_K, blkmod1_K, blkmod2_K, alpha_K, alpha_rho_K
+        integer :: i, k, l, q
 
-        $:GPU_PARALLEL_LOOP(collapse=3, private='[rho_K, G_K, K_K, C_num, pres_K, blkmod1_K, blkmod2_K]')
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[rho_K, G_K, K_K, C_num, pres_K, blkmod1_K, blkmod2_K, alpha_K, alpha_rho_K]')
         do q = 0, p
             do l = 0, n
                 do k = 0, m
@@ -627,10 +627,12 @@ contains
                         ! Same two-component K as the HLLD anchor state (see m_riemann_solver_hypo_hlld.fpp), including the
                         ! verysmall denominator regularization
                         pres_K = q_prim_vf(eqn_idx%E)%sf(k, l, q)
-                        call s_phase_bulk_modulus(pres_K, q_prim_vf(eqn_idx%adv%beg)%sf(k, l, q), &
-                                                  & q_prim_vf(eqn_idx%cont%beg)%sf(k, l, q), 1, blkmod1_K)
-                        call s_phase_bulk_modulus(pres_K, q_prim_vf(eqn_idx%adv%end)%sf(k, l, q), &
-                                                  & q_prim_vf(eqn_idx%cont%end)%sf(k, l, q), 2, blkmod2_K)
+                        alpha_K = q_prim_vf(eqn_idx%adv%beg)%sf(k, l, q)
+                        alpha_rho_K = q_prim_vf(eqn_idx%cont%beg)%sf(k, l, q)
+                        call s_phase_bulk_modulus(pres_K, alpha_K, alpha_rho_K, 1, blkmod1_K)
+                        alpha_K = q_prim_vf(eqn_idx%adv%end)%sf(k, l, q)
+                        alpha_rho_K = q_prim_vf(eqn_idx%cont%end)%sf(k, l, q)
+                        call s_phase_bulk_modulus(pres_K, alpha_K, alpha_rho_K, 2, blkmod2_K)
                         blkmod1_K = blkmod1_K + (4._wp/3._wp)*Gs_hypo(1)
                         blkmod2_K = blkmod2_K + (4._wp/3._wp)*Gs_hypo(2)
                         K_K = q_prim_vf(eqn_idx%adv%beg)%sf(k, l, q)*q_prim_vf(eqn_idx%adv%end)%sf(k, l, &

src/simulation/m_ibm.fpp

@@ -171,6 +171,7 @@ contains
             real(wp), dimension(nb*nnode)    :: presb_IP, massv_IP
             real(wp), dimension(num_species) :: Ys_IP
         #:endif
+        real(wp) :: alpha_q, alpha_rho_q, e_q
         real(wp) :: T_IP, mw_IP, e_IP  !< Image-point temperature, mixture MW, and mass-specific internal energy (chemistry)
         real(wp) :: v_blow_eff         !< Effective surface blowing speed (after any pressure-coupled burn-rate scaling)
         ! Primitive variables at the image point associated with a ghost point, interpolated from surrounding fluid cells.
@@ -223,7 +224,7 @@ contains
             $:GPU_PARALLEL_LOOP(private='[i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, vel_g, vel_norm_IP, &
                                 & r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, &
                                 & innerp, norm, buf, radial_vector, rotation_velocity, j, k, l, q, qv_K, c_IP, nbub, patch_id, &
-                                & Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff, vel_sum_g, E_ghost]')
+                                & Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff, vel_sum_g, E_ghost, alpha_q, alpha_rho_q, e_q]')
             do i = 1, num_gps
                 gp = ghost_points(i)
                 j = gp%loc(1)
@@ -444,9 +445,10 @@ contains
                 if (model_eqns == model_eqns_6eq) then
                     $:GPU_LOOP(parallelism='[seq]')
                     do q = eqn_idx%int_en%beg, eqn_idx%int_en%end
-                        call s_phase_internal_energy(pres_IP, alpha_IP(q - eqn_idx%int_en%beg + 1), &
-                                                     & alpha_rho_IP(q - eqn_idx%int_en%beg + 1), q - eqn_idx%int_en%beg + 1, &
-                                                     & q_cons_vf(q)%sf(j, k, l))
+                        alpha_q = alpha_IP(q - eqn_idx%int_en%beg + 1)
+                        alpha_rho_q = alpha_rho_IP(q - eqn_idx%int_en%beg + 1)
+                        call s_phase_internal_energy(pres_IP, alpha_q, alpha_rho_q, q - eqn_idx%int_en%beg + 1, e_q)
+                        q_cons_vf(q)%sf(j, k, l) = e_q
                     end do
                 end if
             end do

src/simulation/m_pressure_relaxation.fpp

@@ -123,7 +123,7 @@ contains
         #:else
             real(wp), dimension(num_fluids) :: pres_K_init, rho_K_init, rho_K_s
         #:endif
-        real(wp)           :: gamma_K, pi_inf_K, dpi_K, dgamma_K, c2_K
+        real(wp)           :: gamma_K, pi_inf_K, dpi_K, dgamma_K, c2_K, alpha_i, alpha_rho_i, rho_i, p_i, rho_s_i
         integer, parameter :: MAX_ITER = 50
         ! Pressure relaxation convergence tolerance
         real(wp), parameter :: TOLERANCE = 1.e-10_wp
@@ -137,9 +137,10 @@ contains
                 ! Phasic internal energy carries the formation energy: alpha_rho_k*qv_k must be
                 ! removed before inverting the stiffened-gas EOS, or a nonzero qv inflates the
                 ! phasic pressure by rho_k*qv_k/gamma_k (this is what breaks the reactive burn).
-                call s_phase_coefficients(q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l), &
-                                          & q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l), i, rho_K_init(i), gamma_K, pi_inf_K, &
-                                          & dpi_K, dgamma_K)
+                alpha_rho_i = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)
+                alpha_i = q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l)
+                call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_K, pi_inf_K, dpi_K, dgamma_K)
+                rho_K_init(i) = rho_i
                 pres_K_init(i) = ((q_cons_vf(i + eqn_idx%int_en%beg - 1)%sf(j, k, l) - q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, &
                             & k, l)*qvs(i))/q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l) - pi_inf_K)/gamma_K
                 if (.not. f_is_state_dependent(i)) then
@@ -174,7 +175,10 @@ contains
                 $:GPU_LOOP(parallelism='[seq]')
                 do i = 1, num_fluids
                     if (q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l) > sgm_eps .and. any_state_dependent_eos) then
-                        call s_phase_density_on_isentrope(i, rho_K_init(i), pres_K_init(i), pres_relax, rho_K_s(i), c2_K)
+                        rho_i = rho_K_init(i)
+                        p_i = pres_K_init(i)
+                        call s_phase_density_on_isentrope(i, rho_i, p_i, pres_relax, rho_s_i, c2_K)
+                        rho_K_s(i) = rho_s_i
                         f_pres = f_pres + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)/rho_K_s(i)
                         df_pres = df_pres - q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)/(rho_K_s(i)**2*c2_K)
                     else if (q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l) > sgm_eps) then
@@ -207,7 +211,7 @@ contains
         type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
         integer, intent(in)                                    :: j, k, l
         real(wp), intent(in)                                   :: rho, gamma, pi_inf, qv_mix
-        real(wp)                                               :: dyn_pres, pres_relax
+        real(wp)                                               :: dyn_pres, pres_relax, alpha_i, alpha_rho_i, e_i
         integer                                                :: i
 
         dyn_pres = 0._wp
@@ -220,9 +224,10 @@ contains
 
         $:GPU_LOOP(parallelism='[seq]')
         do i = 1, num_fluids
-            call s_phase_internal_energy(pres_relax, q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l), &
-                                         & q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l), i, &
-                                         & q_cons_vf(i + eqn_idx%int_en%beg - 1)%sf(j, k, l))
+            alpha_i = q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l)
+            alpha_rho_i = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)
+            call s_phase_internal_energy(pres_relax, alpha_i, alpha_rho_i, i, e_i)
+            q_cons_vf(i + eqn_idx%int_en%beg - 1)%sf(j, k, l) = e_i
         end do
 
     end subroutine s_correct_internal_energies