Skip to content

Host-capture Re_size in viscous and non-Newtonian device code #2001

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

Draft
sbryngelson wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

Host-capture Re_size in viscous and non-Newtonian device code #2001

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
7 changes: 4 additions & 3 deletions src/simulation/m_hb_function.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,23 +57,24 @@ 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]')

real(wp), dimension(*), intent(in) :: alpha
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

! Plain serial loops: this is already a seq device routine, so no GPU_LOOP
! 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
Expand Down
30 changes: 17 additions & 13 deletions src/simulation/m_riemann_state.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
33 changes: 20 additions & 13 deletions src/simulation/m_viscous.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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

Expand All @@ -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

Expand All @@ -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]')

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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, &
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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)
Expand Down
Loading