diff --git a/src/simulation/m_hb_function.fpp b/src/simulation/m_hb_function.fpp index bd0d19da7..1b7d1a838 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 aeb34c8a9..8949d4861 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 37afd8988..07ef40051 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)