From 437af6f51073badce7177df81350a6ad8263c253 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 9 Oct 2026 23:18:07 -0400 Subject: [PATCH] Host-capture Re_size in viscous and non-Newtonian device code #1588 worked around amdflang reading the declare-target Re_size stale (as 0) across translation units by copying it into firstprivate locals in the Riemann solvers. Device code in other modules still read the module variable directly, so on AMD flang it could silently drop viscosity the same way: - s_compute_axis_inv_re (m_viscous), called from the four cylindrical axis stress kernels - s_compute_mixture_inv_re (m_hb_function), called from the axis routine and from the cylindrical and Cartesian viscous source-flux kernels in m_riemann_state Apply the same pattern: each enclosing kernel takes Re_size_loc1/2 as firstprivate host copies and passes them down, and the device routines use merge(Re_size_loc1, Re_size_loc2, i == 1) in place of Re_size(i), as s_compute_interface_reynolds already does. Values are identical to Re_size, so results do not change on CPU or NVIDIA: viscous, non-Newtonian and cylindrical tests pass unchanged and the NVHPC OpenACC build compiles. Not tested on AMD hardware. --- src/simulation/m_hb_function.fpp | 7 ++++--- src/simulation/m_riemann_state.fpp | 30 +++++++++++++++------------ src/simulation/m_viscous.fpp | 33 ++++++++++++++++++------------ 3 files changed, 41 insertions(+), 29 deletions(-) diff --git a/src/simulation/m_hb_function.fpp b/src/simulation/m_hb_function.fpp index bd0d19da7e..1b7d1a838f 100644 --- a/src/simulation/m_hb_function.fpp +++ b/src/simulation/m_hb_function.fpp @@ -57,7 +57,7 @@ contains !> Mixture inverse Reynolds (= 1/mu_mix) per direction (1=shear, 2=bulk) at one state. Newtonian fluids accumulate the legacy !! per-fluid alpha/Res term. For non-Newtonian fluids only the shear direction (i==1) uses the Herschel-Bulkley viscosity at the !! given shear rate; the bulk direction (i==2) contributes zero, since non-Newtonian bulk viscosity is not supported. - subroutine s_compute_mixture_inv_re(alpha, shear_rate, Res, Re_out) + subroutine s_compute_mixture_inv_re(alpha, shear_rate, Res, Re_out, Re_size_loc1, Re_size_loc2) $:GPU_ROUTINE(parallelism='[seq]') @@ -65,6 +65,7 @@ contains real(wp), intent(in) :: shear_rate real(wp), dimension(2,*), intent(in) :: Res real(wp), dimension(2), intent(out) :: Re_out + integer, intent(in) :: Re_size_loc1, Re_size_loc2 integer :: i, q, fl real(wp) :: mu_q @@ -72,8 +73,8 @@ contains ! directives (they emit empty on Cray/AMD and are redundant here). do i = 1, 2 Re_out(i) = dflt_real - if (Re_size(i) > 0) Re_out(i) = 0._wp - do q = 1, Re_size(i) + if (merge(Re_size_loc1, Re_size_loc2, i == 1) > 0) Re_out(i) = 0._wp + do q = 1, merge(Re_size_loc1, Re_size_loc2, i == 1) fl = Re_idx(i, q) if (any_non_newtonian .and. is_non_newtonian(fl)) then if (i == 1) then diff --git a/src/simulation/m_riemann_state.fpp b/src/simulation/m_riemann_state.fpp index 1c589d4f7f..de7ff88396 100644 --- a/src/simulation/m_riemann_state.fpp +++ b/src/simulation/m_riemann_state.fpp @@ -810,10 +810,12 @@ contains real(wp), dimension(2) :: Re_nn real(wp), dimension(num_fluids) :: alpha_avg integer :: fl + integer :: Re_size_loc1, Re_size_loc2 !< host copies; amdflang reads Re_size stale cross-TU + Re_size_loc1 = Re_size(1); Re_size_loc2 = Re_size(2) $:GPU_PARALLEL_LOOP(collapse=3, private='[idx_rp, avg_v_int, avg_dvdx_int, avg_dvdy_int, avg_dvdz_int, Re_s, Re_b, & & vel_src_int, r_eff, divergence_cyl, stress_vector_shear, stress_normal_bulk, div_v_term_const, & - & gamma_dot, D_xx, D_yy, D_zz, D_xy, D_xz, D_yz, Re_nn, alpha_avg, fl]') + & gamma_dot, D_xx, D_yy, D_zz, D_xy, D_xz, D_yz, Re_nn, alpha_avg, fl]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = iz%beg, iz%end do k = iy%beg, iy%end do j = ix%beg, ix%end @@ -870,7 +872,7 @@ contains ! Raw cell-centered alphas can under/overshoot near interfaces; clamp to [0,1] alpha_avg(fl) = min(max(alpha_avg(fl), 0._wp), 1._wp) end do - call s_compute_mixture_inv_re(alpha_avg, gamma_dot, Res_gs, Re_nn) + call s_compute_mixture_inv_re(alpha_avg, gamma_dot, Res_gs, Re_nn, Re_size_loc1, Re_size_loc2) end if ! Get Re numbers and interface velocity for viscous work @@ -1011,23 +1013,25 @@ contains real(wp), dimension(num_dims, num_dims) :: current_tau_bulk !< Current bulk stress tensor. real(wp), dimension(num_dims) :: vel_src_at_interface !< Interface velocities (u,v,w) for viscous work. #:endif - integer, dimension(3) :: idx_right_phys !< Physical (j,k,l) indices for right state. - real(wp) :: Re_shear !< Interface shear Reynolds number. - real(wp) :: Re_bulk !< Interface bulk Reynolds number. - integer :: j_loop !< Physical x-index loop iterator. - integer :: k_loop !< Physical y-index loop iterator. - integer :: l_loop !< Physical z-index loop iterator. - integer :: i_dim !< Generic dimension/component iterator. - integer :: vel_comp_idx !< Velocity component iterator (1=u, 2=v, 3=w). - real(wp) :: divergence_v !< Velocity divergence at interface. + integer, dimension(3) :: idx_right_phys !< Physical (j,k,l) indices for right state. + real(wp) :: Re_shear !< Interface shear Reynolds number. + real(wp) :: Re_bulk !< Interface bulk Reynolds number. + integer :: j_loop !< Physical x-index loop iterator. + integer :: k_loop !< Physical y-index loop iterator. + integer :: l_loop !< Physical z-index loop iterator. + integer :: i_dim !< Generic dimension/component iterator. + integer :: vel_comp_idx !< Velocity component iterator (1=u, 2=v, 3=w). + real(wp) :: divergence_v !< Velocity divergence at interface. real(wp) :: gamma_dot, D_xx, D_yy, D_zz, D_xy, D_xz, D_yz real(wp), dimension(2) :: Re_nn real(wp), dimension(num_fluids) :: alpha_avg integer :: fl + integer :: Re_size_loc1, Re_size_loc2 !< host copies; amdflang reads Re_size stale cross-TU + Re_size_loc1 = Re_size(1); Re_size_loc2 = Re_size(2) $:GPU_PARALLEL_LOOP(collapse=3, private='[idx_right_phys, vel_grad_avg, current_tau_shear, current_tau_bulk, & & vel_src_at_interface, Re_shear, Re_bulk, divergence_v, i_dim, vel_comp_idx, gamma_dot, D_xx, D_yy, & - & D_zz, D_xy, D_xz, D_yz, Re_nn, alpha_avg, fl]') + & D_zz, D_xy, D_xz, D_yz, Re_nn, alpha_avg, fl]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l_loop = isz%beg, isz%end do k_loop = isy%beg, isy%end do j_loop = isx%beg, isx%end @@ -1081,7 +1085,7 @@ contains ! Raw cell-centered alphas can under/overshoot near interfaces; clamp to [0,1] alpha_avg(fl) = min(max(alpha_avg(fl), 0._wp), 1._wp) end do - call s_compute_mixture_inv_re(alpha_avg, gamma_dot, Res_gs, Re_nn) + call s_compute_mixture_inv_re(alpha_avg, gamma_dot, Res_gs, Re_nn, Re_size_loc1, Re_size_loc2) end if divergence_v = 0.0_wp diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp index 37afd8988e..07ef40051c 100644 --- a/src/simulation/m_viscous.fpp +++ b/src/simulation/m_viscous.fpp @@ -54,7 +54,7 @@ contains !> Compute viscous stress tensor near cylindrical axis, avoiding 1/r singularity at y_cb(-1)=0 !> Inverse Reynolds numbers of the mixture at one cell. Non-Newtonian fluids need the local shear rate, which comes from the !! velocity gradients; Newtonian ones reduce to the volume-fraction-weighted harmonic mean. - subroutine s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc) + subroutine s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc, Re_size_loc1, Re_size_loc2) $:GPU_ROUTINE(function_name='s_compute_axis_inv_re', parallelism='[seq]', cray_inline=True) @@ -66,6 +66,7 @@ contains #:endif integer, intent(in) :: j, k, l real(wp), dimension(2), intent(out) :: Re_visc + integer, intent(in) :: Re_size_loc1, Re_size_loc2 real(wp) :: gamma_dot_c integer :: i, q @@ -80,15 +81,15 @@ contains & l) + grad_y_vf(3)%sf(j, k, l))) end if #:endif - call s_compute_mixture_inv_re(alpha_visc, gamma_dot_c, Res_viscous, Re_visc) + call s_compute_mixture_inv_re(alpha_visc, gamma_dot_c, Res_viscous, Re_visc, Re_size_loc1, Re_size_loc2) else $:GPU_LOOP(parallelism='[seq]') do i = 1, 2 Re_visc(i) = dflt_real - if (Re_size(i) > 0) Re_visc(i) = 0._wp + if (merge(Re_size_loc1, Re_size_loc2, i == 1) > 0) Re_visc(i) = 0._wp $:GPU_LOOP(parallelism='[seq]') - do q = 1, Re_size(i) + do q = 1, merge(Re_size_loc1, Re_size_loc2, i == 1) Re_visc(i) = alpha_visc(Re_idx(i, q))/Res_viscous(i, q) + Re_visc(i) end do @@ -115,9 +116,11 @@ contains real(wp), dimension(num_dims, num_dims) :: tau_Re #:endif - integer :: i, j, k, l, q !< Generic loop iterator + integer :: i, j, k, l, q !< Generic loop iterator + integer :: Re_size_loc1, Re_size_loc2 !< host copies; amdflang reads Re_size stale cross-TU is1_viscous = ix; is2_viscous = iy; is3_viscous = iz + Re_size_loc1 = Re_size(1); Re_size_loc2 = Re_size(2) $:GPU_UPDATE(device='[is1_viscous, is2_viscous, is3_viscous]') @@ -137,7 +140,7 @@ contains #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 if (shear_stress) then ! Shear stresses $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, q, rho_visc, gamma_visc, pi_inf_visc, qv_visc, & - & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]') + & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 do j = is1_viscous%beg, is1_viscous%end @@ -164,7 +167,8 @@ contains call s_compute_mixture_coefficients(alpha_rho_visc, alpha_visc, rho_visc, gamma_visc, pi_inf_visc, & & qv_visc) if (viscous) then - call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc) + call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc, & + & Re_size_loc1, Re_size_loc2) end if ! Shear stress near cylindrical axis: includes v/r hoop term @@ -191,7 +195,7 @@ contains #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 if (bulk_stress) then ! Bulk stresses $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, q, rho_visc, gamma_visc, pi_inf_visc, qv_visc, & - & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]') + & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 do j = is1_viscous%beg, is1_viscous%end @@ -218,7 +222,8 @@ contains call s_compute_mixture_coefficients(alpha_rho_visc, alpha_visc, rho_visc, gamma_visc, pi_inf_visc, & & qv_visc) if (viscous) then - call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc) + call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc, & + & Re_size_loc1, Re_size_loc2) end if tau_Re(2, 2) = (grad_x_vf(1)%sf(j, k, l) + grad_y_vf(2)%sf(j, k, & @@ -239,7 +244,7 @@ contains #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 if (shear_stress) then ! Shear stresses $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, q, rho_visc, gamma_visc, pi_inf_visc, qv_visc, & - & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]') + & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 do j = is1_viscous%beg, is1_viscous%end @@ -266,7 +271,8 @@ contains call s_compute_mixture_coefficients(alpha_rho_visc, alpha_visc, rho_visc, gamma_visc, pi_inf_visc, & & qv_visc) if (viscous) then - call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc) + call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc, & + & Re_size_loc1, Re_size_loc2) end if tau_Re(2, 2) = -(2._wp/3._wp)*grad_z_vf(3)%sf(j, k, l)/y_cc(k)/Re_visc(1) @@ -290,7 +296,7 @@ contains if (bulk_stress) then ! Bulk stresses $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, q, rho_visc, gamma_visc, pi_inf_visc, qv_visc, & - & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]') + & alpha_visc_sum, alpha_visc, alpha_rho_visc, Re_visc, tau_Re]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 do j = is1_viscous%beg, is1_viscous%end @@ -317,7 +323,8 @@ contains call s_compute_mixture_coefficients(alpha_rho_visc, alpha_visc, rho_visc, gamma_visc, pi_inf_visc, & & qv_visc) if (viscous) then - call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc) + call s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, j, k, l, Re_visc, & + & Re_size_loc1, Re_size_loc2) end if tau_Re(2, 2) = grad_z_vf(3)%sf(j, k, l)/y_cc(k)/Re_visc(2)