From dfb5fced9f99c361f8b227ca5f7f41a833cfa0a8 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Mon, 26 May 2025 19:16:52 -0400 Subject: [PATCH 1/8] refac that routine --- src/simulation/m_riemann_solvers.fpp | 614 ++++++++------------------- 1 file changed, 166 insertions(+), 448 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index df388c479a..f942ef3749 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -4467,22 +4467,22 @@ contains end subroutine s_compute_cylindrical_viscous_source_flux - !> The goal of this subroutine is to evaluate and account - !! for the contribution of viscous stresses in the source - !! flux for the momentum and energy. - !! @param velL_vf Left, WENO reconstructed, cell-boundary values of the velocity - !! @param velR_vf Right, WENO reconstructed, cell-boundary values of the velocity - !! @param dvelL_dx_vf Left, WENO reconstructed cell-avg. x-dir derivative of the velocity - !! @param dvelL_dy_vf Left, WENO reconstructed cell-avg. y-dir derivative of the velocity - !! @param dvelL_dz_vf Left, WENO reconstructed cell-avg. z-dir derivative of the velocity - !! @param dvelR_dx_vf Right, WENO reconstructed cell-avg. x-dir derivative of the velocity - !! @param dvelR_dy_vf Right, WENO reconstructed cell-avg. y-dir derivative of the velocity - !! @param dvelR_dz_vf Right, WENO reconstructed cell-avg. z-dir derivative of the velocity - !! @param flux_src_vf Intercell flux - !! @param norm_dir Dimensional splitting coordinate direction - !! @param ix Index bounds in first coordinate direction - !! @param iy Index bounds in second coordinate direction - !! @param iz Index bounds in third coordinate direction + !> @brief Computes Cartesian viscous source flux contributions for momentum and energy. + !! Calculates averaged velocity gradients, gets Re and interface velocities, + !! calls helpers for shear/bulk stress, then updates `flux_src_vf`. + !! @param[in] velL_vf Left boundary velocity (num_dims scalar_field). + !! @param[in] dvelL_dx_vf Left boundary d(vel)/dx (num_dims scalar_field). + !! @param[in] dvelL_dy_vf Left boundary d(vel)/dy (num_dims scalar_field). + !! @param[in] dvelL_dz_vf Left boundary d(vel)/dz (num_dims scalar_field). + !! @param[in] velR_vf Right boundary velocity (num_dims scalar_field). + !! @param[in] dvelR_dx_vf Right boundary d(vel)/dx (num_dims scalar_field). + !! @param[in] dvelR_dy_vf Right boundary d(vel)/dy (num_dims scalar_field). + !! @param[in] dvelR_dz_vf Right boundary d(vel)/dz (num_dims scalar_field). + !! @param[inout] flux_src_vf Intercell source flux array to update (sys_size scalar_field). + !! @param[in] norm_dir Interface normal direction (1=x, 2=y, 3=z). + !! @param[in] ix X-direction loop bounds (int_bounds_info). + !! @param[in] iy Y-direction loop bounds (int_bounds_info). + !! @param[in] iz Z-direction loop bounds (int_bounds_info). subroutine s_compute_cartesian_viscous_source_flux(velL_vf, & dvelL_dx_vf, & dvelL_dy_vf, & @@ -4495,469 +4495,187 @@ contains norm_dir, & ix, iy, iz) - type(scalar_field), & - dimension(num_dims), & - intent(in) :: velL_vf, velR_vf, & - dvelL_dx_vf, dvelR_dx_vf, & - dvelL_dy_vf, dvelR_dy_vf, & - dvelL_dz_vf, dvelR_dz_vf - - type(scalar_field), & - dimension(sys_size), & - intent(inout) :: flux_src_vf - + ! Arguments + type(scalar_field), dimension(num_dims), intent(in) :: velL_vf, velR_vf + type(scalar_field), dimension(num_dims), intent(in) :: dvelL_dx_vf, dvelR_dx_vf + type(scalar_field), dimension(num_dims), intent(in) :: dvelL_dy_vf, dvelR_dy_vf + type(scalar_field), dimension(num_dims), intent(in) :: dvelL_dz_vf, dvelR_dz_vf + type(scalar_field), dimension(sys_size), intent(inout) :: flux_src_vf integer, intent(in) :: norm_dir type(int_bounds_info), intent(in) :: ix, iy, iz - ! Arithmetic mean of the left and right, WENO-reconstructed, cell- - ! boundary values of cell-average first-order spatial derivatives - ! of velocity - real(wp), dimension(num_dims) :: dvel_avg_dx - real(wp), dimension(num_dims) :: dvel_avg_dy - real(wp), dimension(num_dims) :: dvel_avg_dz - - real(wp), dimension(num_dims, num_dims) :: tau_Re !< Viscous stress tensor - - integer :: i, j, k, l !< Generic loop iterators - - ! Viscous Stresses in x-direction - if (norm_dir == 1) then - - if (shear_stress) then ! Shear stresses - !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) - - tau_Re(1, 1) = (4._wp/3._wp)*dvel_avg_dx(1)/ & - Re_avg_rsx_vf(j, k, l, 1) - - flux_src_vf(momxb)%sf(j, k, l) = & - flux_src_vf(momxb)%sf(j, k, l) - & - tau_Re(1, 1) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, 1)* & - tau_Re(1, 1) - - end do - end do - end do - end if - - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) - - tau_Re(1, 1) = dvel_avg_dx(1)/ & - Re_avg_rsx_vf(j, k, l, 2) - - flux_src_vf(momxb)%sf(j, k, l) = & - flux_src_vf(momxb)%sf(j, k, l) - & - tau_Re(1, 1) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, 1)* & - tau_Re(1, 1) - - end do - end do - end do - end if - - if (n == 0) return - - if (shear_stress) then ! Shear stresses - !$acc parallel loop collapse(3) gang vector default(present) private(dvel_avg_dx, dvel_avg_dy, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - !$acc loop seq - do i = 1, 2 - dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j + 1, k, l)) - end do - - dvel_avg_dx(2) = 5e-1_wp*(dvelL_dx_vf(2)%sf(j, k, l) & - + dvelR_dx_vf(2)%sf(j + 1, k, l)) - - tau_Re(1, 1) = -(2._wp/3._wp)*dvel_avg_dy(2)/ & - Re_avg_rsx_vf(j, k, l, 1) - - tau_Re(1, 2) = (dvel_avg_dy(1) + dvel_avg_dx(2))/ & - Re_avg_rsx_vf(j, k, l, 1) - - !$acc loop seq - do i = 1, 2 - - flux_src_vf(contxe + i)%sf(j, k, l) = & - flux_src_vf(contxe + i)%sf(j, k, l) - & - tau_Re(1, i) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, i)* & - tau_Re(1, i) - - end do - - end do - end do - end do - end if - - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dy, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j + 1, k, l)) - - tau_Re(1, 1) = dvel_avg_dy(2)/ & - Re_avg_rsx_vf(j, k, l, 2) - - flux_src_vf(momxb)%sf(j, k, l) = & - flux_src_vf(momxb)%sf(j, k, l) - & - tau_Re(1, 1) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, 1)* & - tau_Re(1, 1) - - end do + ! Local variables + real(wp), dimension(num_dims, num_dims) :: vel_grad_avg !< Averaged velocity gradient tensor `d(vel_i)/d(coord_j)`. + real(wp), dimension(num_dims, num_dims) :: current_tau_shear !< Current shear stress tensor. + 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. + 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. + + !$acc parallel loop collapse(3) gang vector default(present) & + !$acc private(idx_right_phys, vel_grad_avg, & + !$acc current_tau_shear, current_tau_bulk, vel_src_at_interface, & + !$acc Re_shear, Re_bulk, divergence_v, i_dim, vel_comp_idx) + do l_loop = isz%beg, isz%end + do k_loop = isy%beg, isy%end + do j_loop = isx%beg, isx%end + + idx_right_phys(1) = j_loop + idx_right_phys(2) = k_loop + idx_right_phys(3) = l_loop + idx_right_phys(norm_dir) = idx_right_phys(norm_dir) + 1 + + vel_grad_avg = 0.0_wp + do vel_comp_idx = 1, num_dims + vel_grad_avg(vel_comp_idx, 1) = 0.5_wp*(dvelL_dx_vf(vel_comp_idx)%sf(j_loop, k_loop, l_loop) + & + dvelR_dx_vf(vel_comp_idx)%sf(idx_right_phys(1), idx_right_phys(2), idx_right_phys(3))) + if (num_dims > 1) then + vel_grad_avg(vel_comp_idx, 2) = 0.5_wp*(dvelL_dy_vf(vel_comp_idx)%sf(j_loop, k_loop, l_loop) + & + dvelR_dy_vf(vel_comp_idx)%sf(idx_right_phys(1), idx_right_phys(2), idx_right_phys(3))) + end if + if (num_dims > 2) then + vel_grad_avg(vel_comp_idx, 3) = 0.5_wp*(dvelL_dz_vf(vel_comp_idx)%sf(j_loop, k_loop, l_loop) + & + dvelR_dz_vf(vel_comp_idx)%sf(idx_right_phys(1), idx_right_phys(2), idx_right_phys(3))) + end if end do - end do - end if - - if (p == 0) return - - if (shear_stress) then ! Shear stresses - !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - !$acc loop seq - do i = 1, 3, 2 - dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j + 1, k, l)) - end do - - dvel_avg_dx(3) = 5e-1_wp*(dvelL_dx_vf(3)%sf(j, k, l) & - + dvelR_dx_vf(3)%sf(j + 1, k, l)) - - tau_Re(1, 1) = -(2._wp/3._wp)*dvel_avg_dz(3)/ & - Re_avg_rsx_vf(j, k, l, 1) - - tau_Re(1, 3) = (dvel_avg_dz(1) + dvel_avg_dx(3))/ & - Re_avg_rsx_vf(j, k, l, 1) - - !$acc loop seq - do i = 1, 3, 2 - flux_src_vf(contxe + i)%sf(j, k, l) = & - flux_src_vf(contxe + i)%sf(j, k, l) - & - tau_Re(1, i) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, i)* & - tau_Re(1, i) - end do - - end do + divergence_v = 0.0_wp + do i_dim = 1, num_dims + divergence_v = divergence_v + vel_grad_avg(i_dim, i_dim) end do - end do - end if - - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j + 1, k, l)) - - tau_Re(1, 1) = dvel_avg_dz(3)/ & - Re_avg_rsx_vf(j, k, l, 2) - - flux_src_vf(momxb)%sf(j, k, l) = & - flux_src_vf(momxb)%sf(j, k, l) - & - tau_Re(1, 1) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, 1)* & - tau_Re(1, 1) + vel_src_at_interface = 0.0_wp + if (norm_dir == 1) then + Re_shear = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 1) + Re_bulk = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 2) + do i_dim = 1, num_dims + vel_src_at_interface(i_dim) = vel_src_rsx_vf(j_loop, k_loop, l_loop, i_dim) end do - end do - end do - end if - ! END: Viscous Stresses in x-direction - - ! Viscous Stresses in y-direction - elseif (norm_dir == 2) then - - if (shear_stress) then ! Shear stresses - !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - !$acc loop seq - do i = 1, 2 - - dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k + 1, l)) - - dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k + 1, l)) - - end do - - tau_Re(2, 1) = (dvel_avg_dy(1) + dvel_avg_dx(2))/ & - Re_avg_rsy_vf(k, j, l, 1) - - tau_Re(2, 2) = (4._wp*dvel_avg_dy(2) & - - 2._wp*dvel_avg_dx(1))/ & - (3._wp*Re_avg_rsy_vf(k, j, l, 1)) - - !$acc loop seq - do i = 1, 2 - - flux_src_vf(contxe + i)%sf(j, k, l) = & - flux_src_vf(contxe + i)%sf(j, k, l) - & - tau_Re(2, i) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsy_vf(k, j, l, i)* & - tau_Re(2, i) - - end do - + else if (norm_dir == 2) then + Re_shear = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 1) + Re_bulk = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 2) + do i_dim = 1, num_dims + vel_src_at_interface(i_dim) = vel_src_rsy_vf(k_loop, j_loop, l_loop, i_dim) end do - end do - end do - end if - - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k + 1, l)) - - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k + 1, l)) - - tau_Re(2, 2) = (dvel_avg_dx(1) + dvel_avg_dy(2))/ & - Re_avg_rsy_vf(k, j, l, 2) - - flux_src_vf(momxb + 1)%sf(j, k, l) = & - flux_src_vf(momxb + 1)%sf(j, k, l) - & - tau_Re(2, 2) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsy_vf(k, j, l, 2)* & - tau_Re(2, 2) - + else + Re_shear = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 1) + Re_bulk = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 2) + do i_dim = 1, num_dims + vel_src_at_interface(i_dim) = vel_src_rsz_vf(l_loop, k_loop, j_loop, i_dim) end do - end do - end do - end if - - if (p == 0) return - - if (shear_stress) then ! Shear stresses - !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dy, dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - !$acc loop seq - do i = 2, 3 - dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k + 1, l)) - end do - - dvel_avg_dy(3) = 5e-1_wp*(dvelL_dy_vf(3)%sf(j, k, l) & - + dvelR_dy_vf(3)%sf(j, k + 1, l)) - - tau_Re(2, 2) = -(2._wp/3._wp)*dvel_avg_dz(3)/ & - Re_avg_rsy_vf(k, j, l, 1) - - tau_Re(2, 3) = (dvel_avg_dz(2) + dvel_avg_dy(3))/ & - Re_avg_rsy_vf(k, j, l, 1) - - !$acc loop seq - do i = 2, 3 + end if - flux_src_vf(contxe + i)%sf(j, k, l) = & - flux_src_vf(contxe + i)%sf(j, k, l) - & - tau_Re(2, i) + if (shear_stress) then + current_tau_shear = 0.0_wp + call s_calculate_shear_stress_tensor(vel_grad_avg, Re_shear, divergence_v, current_tau_shear) - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsy_vf(k, j, l, i)* & - tau_Re(2, i) - - end do + do i_dim = 1, num_dims + flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) = & + flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) - current_tau_shear(norm_dir, i_dim) + flux_src_vf(E_idx)%sf(j_loop, k_loop, l_loop) = & + flux_src_vf(E_idx)%sf(j_loop, k_loop, l_loop) - & + vel_src_at_interface(i_dim)*current_tau_shear(norm_dir, i_dim) end do - end do - end do - end if - - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end + end if - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k + 1, l)) + if (bulk_stress) then + current_tau_bulk = 0.0_wp + call s_calculate_bulk_stress_tensor(Re_bulk, divergence_v, current_tau_bulk) - tau_Re(2, 2) = dvel_avg_dz(3)/ & - Re_avg_rsy_vf(k, j, l, 2) - - flux_src_vf(momxb + 1)%sf(j, k, l) = & - flux_src_vf(momxb + 1)%sf(j, k, l) - & - tau_Re(2, 2) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsy_vf(k, j, l, 2)* & - tau_Re(2, 2) + do i_dim = 1, num_dims + flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) = & + flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) - current_tau_bulk(norm_dir, i_dim) + flux_src_vf(E_idx)%sf(j_loop, k_loop, l_loop) = & + flux_src_vf(E_idx)%sf(j_loop, k_loop, l_loop) - & + vel_src_at_interface(i_dim)*current_tau_bulk(norm_dir, i_dim) end do - end do - end do - end if - ! END: Viscous Stresses in y-direction - - ! Viscous Stresses in z-direction - else - - if (shear_stress) then ! Shear stresses - !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - !$acc loop seq - do i = 1, 3, 2 - dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k, l + 1)) - end do - - !$acc loop seq - do i = 2, 3 - dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k, l + 1)) - end do - - !$acc loop seq - do i = 1, 3 - dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k, l + 1)) - end do - - tau_Re(3, 1) = (dvel_avg_dz(1) + dvel_avg_dx(3))/ & - Re_avg_rsz_vf(l, k, j, 1) - - tau_Re(3, 2) = (dvel_avg_dz(2) + dvel_avg_dy(3))/ & - Re_avg_rsz_vf(l, k, j, 1) - - tau_Re(3, 3) = (4._wp*dvel_avg_dz(3) & - - 2._wp*dvel_avg_dx(1) & - - 2._wp*dvel_avg_dy(2))/ & - (3._wp*Re_avg_rsz_vf(l, k, j, 1)) - - !$acc loop seq - do i = 1, 3 - - flux_src_vf(contxe + i)%sf(j, k, l) = & - flux_src_vf(contxe + i)%sf(j, k, l) - & - tau_Re(3, i) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsz_vf(l, k, j, i)* & - tau_Re(3, i) - - end do + end if - end do - end do end do - end if + end do + end do + !$acc end parallel loop - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end + end subroutine s_compute_cartesian_viscous_source_flux - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k, l + 1)) + !> @brief Calculates shear stress tensor components. + !! tau_ij_shear = ( (dui/dxj + duj/dxi) - (2/3)*(div_v)*delta_ij ) / Re_shear + !! @param[in] vel_grad_avg Averaged velocity gradient tensor (d(vel_i)/d(coord_j)). + !! @param[in] Re_shear Shear Reynolds number. + !! @param[in] divergence_v Velocity divergence (du/dx + dv/dy + dw/dz). + !! @param[out] tau_shear_out Calculated shear stress tensor (stress on i-face, j-direction). + subroutine s_calculate_shear_stress_tensor(vel_grad_avg, Re_shear, divergence_v, tau_shear_out) + use m_global_parameters, only: wp, num_dims, sgm_eps + implicit none + + ! Arguments + real(wp), dimension(num_dims, num_dims), intent(in) :: vel_grad_avg + real(wp), intent(in) :: Re_shear + real(wp), intent(in) :: divergence_v + real(wp), dimension(num_dims, num_dims), intent(out) :: tau_shear_out + + ! Local variables + integer :: i_dim !< Loop iterator for face normal. + integer :: j_dim !< Loop iterator for force component direction. + + tau_shear_out = 0.0_wp + + if (abs(Re_shear) < sgm_eps) then + return + end if - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k, l + 1)) + do i_dim = 1, num_dims + do j_dim = 1, num_dims + tau_shear_out(i_dim, j_dim) = (vel_grad_avg(j_dim, i_dim) + vel_grad_avg(i_dim, j_dim))/Re_shear + if (i_dim == j_dim) then + tau_shear_out(i_dim, j_dim) = tau_shear_out(i_dim, j_dim) - & + (2.0_wp/3.0_wp)*divergence_v/Re_shear + end if + end do + end do - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k, l + 1)) + end subroutine s_calculate_shear_stress_tensor - tau_Re(3, 3) = (dvel_avg_dx(1) & - + dvel_avg_dy(2) & - + dvel_avg_dz(3))/ & - Re_avg_rsz_vf(l, k, j, 2) + !> @brief Calculates bulk stress tensor components (diagonal only). + !! tau_ii_bulk = (div_v) / Re_bulk. Off-diagonals are zero. + !! @param[in] Re_bulk Bulk Reynolds number. + !! @param[in] divergence_v Velocity divergence (du/dx + dv/dy + dw/dz). + !! @param[out] tau_bulk_out Calculated bulk stress tensor (stress on i-face, i-direction). + subroutine s_calculate_bulk_stress_tensor(Re_bulk, divergence_v, tau_bulk_out) + use m_global_parameters, only: wp, num_dims, sgm_eps + implicit none - flux_src_vf(momxe)%sf(j, k, l) = & - flux_src_vf(momxe)%sf(j, k, l) - & - tau_Re(3, 3) + ! Arguments + real(wp), intent(in) :: Re_bulk + real(wp), intent(in) :: divergence_v + real(wp), dimension(num_dims, num_dims), intent(out) :: tau_bulk_out - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsz_vf(l, k, j, 3)* & - tau_Re(3, 3) + ! Local variables + integer :: i_dim !< Loop iterator for diagonal components. - end do - end do - end do - end if + tau_bulk_out = 0.0_wp + if (abs(Re_bulk) < sgm_eps) then + return end if - ! END: Viscous Stresses in z-direction - end subroutine s_compute_cartesian_viscous_source_flux + do i_dim = 1, num_dims + tau_bulk_out(i_dim, i_dim) = divergence_v/Re_bulk + end do + + end subroutine s_calculate_bulk_stress_tensor !> Deallocation and/or disassociation procedures that are !! needed to finalize the selected Riemann problem solver From bf64aa52995139c374629a7114342fd8ab1903fb Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 27 May 2025 14:02:44 -0400 Subject: [PATCH 2/8] get the cylindrical part --- src/simulation/m_riemann_solvers.fpp | 814 +++++++-------------------- 1 file changed, 210 insertions(+), 604 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index f942ef3749..b6376fe72a 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -3940,560 +3940,181 @@ contains end subroutine s_initialize_riemann_solver - !> The goal of this subroutine is to evaluate and account - !! for the contribution of viscous stresses in the source - !! flux for the momentum and energy. - !! @param velL_vf Left, WENO reconstructed, cell-boundary values of the velocity - !! @param velR_vf Right, WENO reconstructed, cell-boundary values of the velocity - !! @param dvelL_dx_vf Left, WENO reconstructed cell-avg. x-dir derivative of the velocity - !! @param dvelL_dy_vf Left, WENO reconstructed cell-avg. y-dir derivative of the velocity - !! @param dvelL_dz_vf Left, WENO reconstructed cell-avg. z-dir derivative of the velocity - !! @param dvelR_dx_vf Right, WENO reconstructed cell-avg. x-dir derivative of the velocity - !! @param dvelR_dy_vf Right, WENO reconstructed cell-avg. y-dir derivative of the velocity - !! @param dvelR_dz_vf Right, WENO reconstructed cell-avg. z-dir derivative of the velocity - !! @param flux_src_vf Intercell flux - !! @param norm_dir Dimensional splitting coordinate direction - !! @param ix Index bounds in first coordinate direction - !! @param iy Index bounds in second coordinate direction - !! @param iz Index bounds in third coordinate direction + !> @brief Computes cylindrical viscous source flux contributions for momentum and energy. + !! Calculates Cartesian components of the stress tensor using averaged velocity derivatives + !! and cylindrical geometric factors, then updates `flux_src_vf`. + !! Assumes x-dir is axial (z_cyl), y-dir is radial (r_cyl), z-dir is azimuthal (theta_cyl for derivatives). + !! @param[in] velL_vf Left boundary velocity ($v_x, v_y, v_z$) (num_dims scalar_field). + !! @param[in] dvelL_dx_vf Left boundary $\partial v_i/\partial x$ (num_dims scalar_field). + !! @param[in] dvelL_dy_vf Left boundary $\partial v_i/\partial y$ (num_dims scalar_field). + !! @param[in] dvelL_dz_vf Left boundary $\partial v_i/\partial z$ (num_dims scalar_field). + !! @param[in] velR_vf Right boundary velocity ($v_x, v_y, v_z$) (num_dims scalar_field). + !! @param[in] dvelR_dx_vf Right boundary $\partial v_i/\partial x$ (num_dims scalar_field). + !! @param[in] dvelR_dy_vf Right boundary $\partial v_i/\partial y$ (num_dims scalar_field). + !! @param[in] dvelR_dz_vf Right boundary $\partial v_i/\partial z$ (num_dims scalar_field). + !! @param[inout] flux_src_vf Intercell source flux array to update (sys_size scalar_field). + !! @param[in] norm_dir Interface normal direction (1=x-face, 2=y-face, 3=z-face). + !! @param[in] ix Global X-direction loop bounds (int_bounds_info). + !! @param[in] iy Global Y-direction loop bounds (int_bounds_info). + !! @param[in] iz Global Z-direction loop bounds (int_bounds_info). subroutine s_compute_cylindrical_viscous_source_flux(velL_vf, & - dvelL_dx_vf, & - dvelL_dy_vf, & - dvelL_dz_vf, & + dvelL_dx_vf, dvelL_dy_vf, dvelL_dz_vf, & velR_vf, & - dvelR_dx_vf, & - dvelR_dy_vf, & - dvelR_dz_vf, & - flux_src_vf, & - norm_dir, & - ix, iy, iz) - - type(scalar_field), & - dimension(num_dims), & - intent(in) :: velL_vf, velR_vf, & - dvelL_dx_vf, dvelR_dx_vf, & - dvelL_dy_vf, dvelR_dy_vf, & - dvelL_dz_vf, dvelR_dz_vf - - type(scalar_field), & - dimension(sys_size), & - intent(inout) :: flux_src_vf + dvelR_dx_vf, dvelR_dy_vf, dvelR_dz_vf, & + flux_src_vf, norm_dir, ix, iy, iz) + type(scalar_field), dimension(num_dims), intent(in) :: velL_vf, velR_vf + type(scalar_field), dimension(num_dims), intent(in) :: dvelL_dx_vf, dvelR_dx_vf + type(scalar_field), dimension(num_dims), intent(in) :: dvelL_dy_vf, dvelR_dy_vf + type(scalar_field), dimension(num_dims), intent(in) :: dvelL_dz_vf, dvelR_dz_vf + type(scalar_field), dimension(sys_size), intent(inout) :: flux_src_vf integer, intent(in) :: norm_dir - type(int_bounds_info), intent(in) :: ix, iy, iz - ! Arithmetic mean of the left and right, WENO-reconstructed, cell- - ! boundary values of cell-average first-order spatial derivatives - ! of velocity - real(wp), dimension(num_dims) :: avg_vel - real(wp), dimension(num_dims) :: dvel_avg_dx - real(wp), dimension(num_dims) :: dvel_avg_dy - real(wp), dimension(num_dims) :: dvel_avg_dz - - ! Viscous stress tensor - real(wp), dimension(num_dims, num_dims) :: tau_Re - - ! Generic loop iterators - integer :: i, j, k, l - - ! Viscous Stresses in z-direction - if (norm_dir == 1) then - if (shear_stress) then ! Shear stresses - !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) - - tau_Re(1, 1) = (4._wp/3._wp)*dvel_avg_dx(1)/ & - Re_avg_rsx_vf(j, k, l, 1) - - flux_src_vf(momxb)%sf(j, k, l) = & - flux_src_vf(momxb)%sf(j, k, l) - & - tau_Re(1, 1) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, 1)* & - tau_Re(1, 1) - - end do - end do - end do - end if - - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) - - tau_Re(1, 1) = dvel_avg_dx(1)/ & - Re_avg_rsx_vf(j, k, l, 2) - - flux_src_vf(momxb)%sf(j, k, l) = & - flux_src_vf(momxb)%sf(j, k, l) - & - tau_Re(1, 1) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, 1)* & - tau_Re(1, 1) - - end do - end do - end do - end if - - if (n == 0) return - - if (shear_stress) then ! Shear stresses - !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j + 1, k, l)) - - !$acc loop seq - do i = 1, 2 - dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j + 1, k, l)) - end do - - dvel_avg_dx(2) = 5e-1_wp*(dvelL_dx_vf(2)%sf(j, k, l) & - + dvelR_dx_vf(2)%sf(j + 1, k, l)) - - tau_Re(1, 1) = -(2._wp/3._wp)*(dvel_avg_dy(2) + & - avg_vel(2)/y_cc(k))/ & - Re_avg_rsx_vf(j, k, l, 1) - - tau_Re(1, 2) = (dvel_avg_dy(1) + dvel_avg_dx(2))/ & - Re_avg_rsx_vf(j, k, l, 1) - - !$acc loop seq - do i = 1, 2 - flux_src_vf(contxe + i)%sf(j, k, l) = & - flux_src_vf(contxe + i)%sf(j, k, l) - & - tau_Re(1, i) - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, i)* & - tau_Re(1, i) - end do - - end do - end do - end do - end if - - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dy, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j + 1, k, l)) - - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j + 1, k, l)) - - tau_Re(1, 1) = (dvel_avg_dy(2) + & - avg_vel(2)/y_cc(k))/ & - Re_avg_rsx_vf(j, k, l, 2) - - flux_src_vf(momxb)%sf(j, k, l) = & - flux_src_vf(momxb)%sf(j, k, l) - & - tau_Re(1, 1) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, 1)* & - tau_Re(1, 1) - - end do - end do - end do - end if - - if (p == 0) return - - if (shear_stress) then ! Shear stresses - !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - !$acc loop seq - do i = 1, 3, 2 - dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j + 1, k, l)) - end do - - dvel_avg_dx(3) = 5e-1_wp*(dvelL_dx_vf(3)%sf(j, k, l) & - + dvelR_dx_vf(3)%sf(j + 1, k, l)) - - tau_Re(1, 1) = -(2._wp/3._wp)*dvel_avg_dz(3)/y_cc(k)/ & - Re_avg_rsx_vf(j, k, l, 1) - - tau_Re(1, 3) = (dvel_avg_dz(1)/y_cc(k) + dvel_avg_dx(3))/ & - Re_avg_rsx_vf(j, k, l, 1) - - !$acc loop seq - do i = 1, 3, 2 - - flux_src_vf(contxe + i)%sf(j, k, l) = & - flux_src_vf(contxe + i)%sf(j, k, l) - & - tau_Re(1, i) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, i)* & - tau_Re(1, i) - - end do - - end do - end do - end do - end if - - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private( avg_vel, dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j + 1, k, l)) - - tau_Re(1, 1) = dvel_avg_dz(3)/y_cc(k)/ & - Re_avg_rsx_vf(j, k, l, 2) - - flux_src_vf(momxb)%sf(j, k, l) = & - flux_src_vf(momxb)%sf(j, k, l) - & - tau_Re(1, 1) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsx_vf(j, k, l, 1)* & - tau_Re(1, 1) - - end do - end do - end do - end if - ! END: Viscous Stresses in z-direction - - ! Viscous Stresses in r-direction - elseif (norm_dir == 2) then - - if (shear_stress) then ! Shear stresses - - !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j, k + 1, l)) - - !$acc loop seq - do i = 1, 2 - - dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k + 1, l)) - - dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k + 1, l)) - - end do - - tau_Re(2, 1) = (dvel_avg_dy(1) + dvel_avg_dx(2))/ & - Re_avg_rsy_vf(k, j, l, 1) - - tau_Re(2, 2) = (4._wp*dvel_avg_dy(2) & - - 2._wp*dvel_avg_dx(1) & - - 2._wp*avg_vel(2)/y_cb(k))/ & - (3._wp*Re_avg_rsy_vf(k, j, l, 1)) - - !$acc loop seq - do i = 1, 2 - - flux_src_vf(contxe + i)%sf(j, k, l) = & - flux_src_vf(contxe + i)%sf(j, k, l) - & - tau_Re(2, i) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsy_vf(k, j, l, i)* & - tau_Re(2, i) - - end do - - end do - end do - end do - end if - - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j, k + 1, l)) - - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k + 1, l)) - - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k + 1, l)) - - tau_Re(2, 2) = (dvel_avg_dx(1) + dvel_avg_dy(2) + & - avg_vel(2)/y_cb(k))/ & - Re_avg_rsy_vf(k, j, l, 2) - - flux_src_vf(momxb + 1)%sf(j, k, l) = & - flux_src_vf(momxb + 1)%sf(j, k, l) - & - tau_Re(2, 2) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsy_vf(k, j, l, 2)* & - tau_Re(2, 2) - - end do - end do - end do - end if - - if (p == 0) return - - if (shear_stress) then ! Shear stresses - !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dy, dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - avg_vel(3) = 5e-1_wp*(velL_vf(3)%sf(j, k, l) & - + velR_vf(3)%sf(j, k + 1, l)) - - !$acc loop seq - do i = 2, 3 - dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k + 1, l)) - end do - - dvel_avg_dy(3) = 5e-1_wp*(dvelL_dy_vf(3)%sf(j, k, l) & - + dvelR_dy_vf(3)%sf(j, k + 1, l)) - - tau_Re(2, 2) = -(2._wp/3._wp)*dvel_avg_dz(3)/y_cb(k)/ & - Re_avg_rsy_vf(k, j, l, 1) - - tau_Re(2, 3) = ((dvel_avg_dz(2) - avg_vel(3))/ & - y_cb(k) + dvel_avg_dy(3))/ & - Re_avg_rsy_vf(k, j, l, 1) - - !$acc loop seq - do i = 2, 3 - - flux_src_vf(contxe + i)%sf(j, k, l) = & - flux_src_vf(contxe + i)%sf(j, k, l) - & - tau_Re(2, i) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsy_vf(k, j, l, i)* & - tau_Re(2, i) - - end do - - end do - end do - end do - end if - - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k + 1, l)) + ! Local variables + real(wp), dimension(num_dims) :: avg_v_int !!< Averaged interface velocity $(v_x, v_y, v_z)$ (grid directions). + real(wp), dimension(num_dims) :: avg_dvdx_int !!< Averaged interface $\partial v_i/\partial x$ (grid dir 1). + real(wp), dimension(num_dims) :: avg_dvdy_int !!< Averaged interface $\partial v_i/\partial y$ (grid dir 2). + real(wp), dimension(num_dims) :: avg_dvdz_int !!< Averaged interface $\partial v_i/\partial z$ (grid dir 3). - tau_Re(2, 2) = dvel_avg_dz(3)/y_cb(k)/ & - Re_avg_rsy_vf(k, j, l, 2) + real(wp), dimension(num_dims) :: stress_vector_shear !!< Shear stress vector $(\sigma_{N1}, \sigma_{N2}, \sigma_{N3})$ on N-face (grid directions). + real(wp) :: stress_normal_bulk !!< Normal bulk stress component $\sigma_{NN}$ on N-face. - flux_src_vf(momxb + 1)%sf(j, k, l) = & - flux_src_vf(momxb + 1)%sf(j, k, l) - & - tau_Re(2, 2) + real(wp) :: Re_s, Re_b !!< Effective interface shear and bulk Reynolds numbers. + real(wp), dimension(num_dims) :: vel_src_int !!< Interface velocity $(v_1,v_2,v_3)$ (grid directions) for viscous work. + real(wp) :: r_eff !!< Effective radius at interface for cylindrical terms. + real(wp) :: div_v_term_const !!< Common term $-(2/3)(\nabla \cdot \mathbf{v}) / \text{Re}_s$ for shear stress diagonal. + real(wp) :: divergence_cyl !!< Full divergence $\nabla \cdot \mathbf{v}$ in cylindrical coordinates. - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsy_vf(k, j, l, 2)* & - tau_Re(2, 2) + integer :: j, k, l !!< Loop iterators for $x, y, z$ grid directions. + integer :: i_vel !!< Loop iterator for velocity components. + integer :: idx_rp(3) !!< Indices $(j,k,l)$ of 'right' point for averaging. - end do + !$acc parallel loop collapse(3) gang vector default(present) & + !$acc private(idx_rp, avg_v_int, avg_dvdx_int, avg_dvdy_int, avg_dvdz_int, & + !$acc Re_s, Re_b, vel_src_int, r_eff, divergence_cyl, & + !$acc stress_vector_shear, stress_normal_bulk, div_v_term_const) + do l = iz%beg, iz%end + do k = iy%beg, iy%end + do j = ix%beg, ix%end + + ! Determine indices for the 'right' state for averaging across the interface + idx_rp = [j, k, l] + idx_rp(norm_dir) = idx_rp(norm_dir) + 1 + + ! Average velocities and their derivatives at the interface + ! For cylindrical: x-dir ~ axial (z_cyl), y-dir ~ radial (r_cyl), z-dir ~ azimuthal (theta_cyl) + do i_vel = 1, num_dims + avg_v_int(i_vel) = 0.5_wp*(velL_vf(i_vel)%sf(j, k, l) + velR_vf(i_vel)%sf(idx_rp(1), idx_rp(2), idx_rp(3))) + + avg_dvdx_int(i_vel) = 0.5_wp*(dvelL_dx_vf(i_vel)%sf(j, k, l) + & + dvelR_dx_vf(i_vel)%sf(idx_rp(1), idx_rp(2), idx_rp(3))) + if (num_dims > 1) then + avg_dvdy_int(i_vel) = 0.5_wp*(dvelL_dy_vf(i_vel)%sf(j, k, l) + & + dvelR_dy_vf(i_vel)%sf(idx_rp(1), idx_rp(2), idx_rp(3))) + else + avg_dvdy_int(i_vel) = 0.0_wp + end if + if (num_dims > 2) then + avg_dvdz_int(i_vel) = 0.5_wp*(dvelL_dz_vf(i_vel)%sf(j, k, l) + & + dvelR_dz_vf(i_vel)%sf(idx_rp(1), idx_rp(2), idx_rp(3))) + else + avg_dvdz_int(i_vel) = 0.0_wp + end if end do - end do - end if - ! END: Viscous Stresses in r-direction - - ! Viscous Stresses in theta-direction - else - - if (shear_stress) then ! Shear stresses - !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - !$acc loop seq - do i = 2, 3 - avg_vel(i) = 5e-1_wp*(velL_vf(i)%sf(j, k, l) & - + velR_vf(i)%sf(j, k, l + 1)) - end do - - !$acc loop seq - do i = 1, 3, 2 - dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k, l + 1)) - end do - do i = 2, 3 - dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k, l + 1)) - end do - - !$acc loop seq - do i = 1, 3 - dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k, l + 1)) - end do - - tau_Re(3, 1) = (dvel_avg_dz(1)/y_cc(k) + dvel_avg_dx(3))/ & - Re_avg_rsz_vf(l, k, j, 1)/ & - y_cc(k) + ! Get Re numbers and interface velocity for viscous work + select case (norm_dir) + case (1) ! x-face (axial face in z_cyl direction) + Re_s = Re_avg_rsx_vf(j, k, l, 1) + Re_b = Re_avg_rsx_vf(j, k, l, 2) + vel_src_int = vel_src_rsx_vf(j, k, l, 1:num_dims) + r_eff = y_cc(k) + case (2) ! y-face (radial face in r_cyl direction) + Re_s = Re_avg_rsy_vf(k, j, l, 1) + Re_b = Re_avg_rsy_vf(k, j, l, 2) + vel_src_int = vel_src_rsy_vf(k, j, l, 1:num_dims) + r_eff = y_cb(k) + case (3) ! z-face (azimuthal face in theta_cyl direction) + Re_s = Re_avg_rsz_vf(l, k, j, 1) + Re_b = Re_avg_rsz_vf(l, k, j, 2) + vel_src_int = vel_src_rsz_vf(l, k, j, 1:num_dims) + r_eff = y_cc(k) + end select + + ! Divergence in cylindrical coordinates (vx=vz_cyl, vy=vr_cyl, vz=vtheta_cyl) + divergence_cyl = avg_dvdx_int(1) + avg_dvdy_int(2) + avg_v_int(2)/r_eff + if (num_dims > 2) then + divergence_cyl = divergence_cyl + avg_dvdz_int(3)/r_eff + end if - tau_Re(3, 2) = ((dvel_avg_dz(2) - avg_vel(3))/ & - y_cc(k) + dvel_avg_dy(3))/ & - Re_avg_rsz_vf(l, k, j, 1)/ & - y_cc(k) + stress_vector_shear = 0.0_wp + stress_normal_bulk = 0.0_wp - tau_Re(3, 3) = (4._wp*dvel_avg_dz(3)/y_cc(k) & - - 2._wp*dvel_avg_dx(1) & - - 2._wp*dvel_avg_dy(2) & - + 4._wp*avg_vel(2)/y_cc(k))/ & - (3._wp*Re_avg_rsz_vf(l, k, j, 1))/ & - y_cc(k) + if (shear_stress) then + div_v_term_const = -(2.0_wp/3.0_wp)*divergence_cyl/Re_s - !$acc loop seq - do i = 1, 3 - flux_src_vf(contxe + i)%sf(j, k, l) = & - flux_src_vf(contxe + i)%sf(j, k, l) - & - tau_Re(3, i) - - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsz_vf(l, k, j, i)* & - tau_Re(3, i) - end do + select case (norm_dir) + case (1) ! X-face (axial normal, z_cyl) + stress_vector_shear(1) = (2.0_wp*avg_dvdx_int(1))/Re_s + div_v_term_const + if (num_dims > 1) then + stress_vector_shear(2) = (avg_dvdy_int(1) + avg_dvdx_int(2))/Re_s + end if + if (num_dims > 2) then + stress_vector_shear(3) = (avg_dvdz_int(1)/r_eff + avg_dvdx_int(3))/Re_s + end if + case (2) ! Y-face (radial normal, r_cyl) + if (num_dims > 1) then + stress_vector_shear(1) = (avg_dvdy_int(1) + avg_dvdx_int(2))/Re_s + stress_vector_shear(2) = (2.0_wp*avg_dvdy_int(2))/Re_s + div_v_term_const + if (num_dims > 2) then + stress_vector_shear(3) = (avg_dvdz_int(2)/r_eff - avg_v_int(3)/r_eff + avg_dvdy_int(3))/Re_s + end if + else + stress_vector_shear(1) = (2.0_wp*avg_dvdx_int(1))/Re_s + div_v_term_const + end if + case (3) ! Z-face (azimuthal normal, theta_cyl) + if (num_dims > 2) then + stress_vector_shear(1) = (avg_dvdz_int(1)/r_eff + avg_dvdx_int(3))/Re_s + stress_vector_shear(2) = (avg_dvdz_int(2)/r_eff - avg_v_int(3)/r_eff + avg_dvdy_int(3))/Re_s + stress_vector_shear(3) = (2.0_wp*(avg_dvdz_int(3)/r_eff + avg_v_int(2)/r_eff))/Re_s + div_v_term_const + end if + end select + do i_vel = 1, num_dims + flux_src_vf(momxb + i_vel - 1)%sf(j, k, l) = flux_src_vf(momxb + i_vel - 1)%sf(j, k, l) - stress_vector_shear(i_vel) + flux_src_vf(E_idx)%sf(j, k, l) = flux_src_vf(E_idx)%sf(j, k, l) - vel_src_int(i_vel)*stress_vector_shear(i_vel) end do - end do - end do - end if - - if (bulk_stress) then ! Bulk stresses - !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) - do l = isz%beg, isz%end - do k = isy%beg, isy%end - do j = isx%beg, isx%end - - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j, k, l + 1)) - - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k, l + 1)) - - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k, l + 1)) - - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k, l + 1)) - - tau_Re(3, 3) = (dvel_avg_dx(1) & - + dvel_avg_dy(2) & - + dvel_avg_dz(3)/y_cc(k) & - + avg_vel(2)/y_cc(k))/ & - Re_avg_rsz_vf(l, k, j, 2)/ & - y_cc(k) + end if - flux_src_vf(momxe)%sf(j, k, l) = & - flux_src_vf(momxe)%sf(j, k, l) - & - tau_Re(3, 3) + if (bulk_stress) then + stress_normal_bulk = divergence_cyl/Re_b - flux_src_vf(E_idx)%sf(j, k, l) = & - flux_src_vf(E_idx)%sf(j, k, l) - & - vel_src_rsz_vf(l, k, j, 3)* & - tau_Re(3, 3) + flux_src_vf(momxb + norm_dir - 1)%sf(j, k, l) = flux_src_vf(momxb + norm_dir - 1)%sf(j, k, l) - stress_normal_bulk + flux_src_vf(E_idx)%sf(j, k, l) = flux_src_vf(E_idx)%sf(j, k, l) - vel_src_int(norm_dir)*stress_normal_bulk + end if - end do - end do end do - end if - - end if - ! END: Viscous Stresses in theta-direction + end do + end do + !$acc end parallel loop end subroutine s_compute_cylindrical_viscous_source_flux - !> @brief Computes Cartesian viscous source flux contributions for momentum and energy. - !! Calculates averaged velocity gradients, gets Re and interface velocities, - !! calls helpers for shear/bulk stress, then updates `flux_src_vf`. - !! @param[in] velL_vf Left boundary velocity (num_dims scalar_field). - !! @param[in] dvelL_dx_vf Left boundary d(vel)/dx (num_dims scalar_field). - !! @param[in] dvelL_dy_vf Left boundary d(vel)/dy (num_dims scalar_field). - !! @param[in] dvelL_dz_vf Left boundary d(vel)/dz (num_dims scalar_field). - !! @param[in] velR_vf Right boundary velocity (num_dims scalar_field). - !! @param[in] dvelR_dx_vf Right boundary d(vel)/dx (num_dims scalar_field). - !! @param[in] dvelR_dy_vf Right boundary d(vel)/dy (num_dims scalar_field). - !! @param[in] dvelR_dz_vf Right boundary d(vel)/dz (num_dims scalar_field). - !! @param[inout] flux_src_vf Intercell source flux array to update (sys_size scalar_field). + !> @brief Computes Cartesian viscous source flux contributions. + !! Calls helpers for shear/bulk stress after averaging gradients. + !! @param[in] velL_vf,velR_vf Left/Right boundary velocities. + !! @param[in] dvelL_dx_vf,dvelR_dx_vf Left/Right boundary $\partial v_i/\partial x$. + !! @param[in] dvelL_dy_vf,dvelR_dy_vf Left/Right boundary $\partial v_i/\partial y$. + !! @param[in] dvelL_dz_vf,dvelR_dz_vf Left/Right boundary $\partial v_i/\partial z$. + !! @param[inout] flux_src_vf Intercell source flux array to update. !! @param[in] norm_dir Interface normal direction (1=x, 2=y, 3=z). - !! @param[in] ix X-direction loop bounds (int_bounds_info). - !! @param[in] iy Y-direction loop bounds (int_bounds_info). - !! @param[in] iz Z-direction loop bounds (int_bounds_info). + !! @param[in] ix,iy,iz Global X,Y,Z-direction loop bounds. subroutine s_compute_cartesian_viscous_source_flux(velL_vf, & - dvelL_dx_vf, & - dvelL_dy_vf, & - dvelL_dz_vf, & + dvelL_dx_vf, dvelL_dy_vf, dvelL_dz_vf, & velR_vf, & - dvelR_dx_vf, & - dvelR_dy_vf, & - dvelR_dz_vf, & - flux_src_vf, & - norm_dir, & - ix, iy, iz) + dvelR_dx_vf, dvelR_dy_vf, dvelR_dz_vf, & + flux_src_vf, norm_dir, ix, iy, iz) ! Arguments type(scalar_field), dimension(num_dims), intent(in) :: velL_vf, velR_vf @@ -4505,30 +4126,28 @@ contains type(int_bounds_info), intent(in) :: ix, iy, iz ! Local variables - real(wp), dimension(num_dims, num_dims) :: vel_grad_avg !< Averaged velocity gradient tensor `d(vel_i)/d(coord_j)`. - real(wp), dimension(num_dims, num_dims) :: current_tau_shear !< Current shear stress tensor. - 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. - 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), dimension(num_dims, num_dims) :: vel_grad_avg !!< Averaged velocity gradient tensor $L_{ij} = \partial v_i/\partial x_j$. + real(wp), dimension(num_dims, num_dims) :: current_tau_shear !!< Current shear stress tensor $\tau_{ij}^{(s)}$. + real(wp), dimension(num_dims, num_dims) :: current_tau_bulk !!< Current bulk stress tensor $\tau_{ij}^{(b)}$ (diagonal). + real(wp), dimension(num_dims) :: vel_src_at_interface !!< Interface velocities $(v_x,v_y,v_z)$ for viscous work. + integer, dimension(3) :: idx_right_phys !!< Physical $(j,k,l)$ indices for 'right' state. + + real(wp) :: Re_shear_int !!< Interface shear Reynolds number. + real(wp) :: Re_bulk_int !!< 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 iterator for flux updates. + integer :: vel_comp_idx !!< Iterator for velocity components. + real(wp) :: divergence_v !!< Velocity divergence $\nabla \cdot \mathbf{v}$ at interface. !$acc parallel loop collapse(3) gang vector default(present) & - !$acc private(idx_right_phys, vel_grad_avg, & - !$acc current_tau_shear, current_tau_bulk, vel_src_at_interface, & - !$acc Re_shear, Re_bulk, divergence_v, i_dim, vel_comp_idx) - do l_loop = isz%beg, isz%end - do k_loop = isy%beg, isy%end - do j_loop = isx%beg, isx%end + !$acc private(idx_right_phys, vel_grad_avg, current_tau_shear, current_tau_bulk, & + !$acc vel_src_at_interface, Re_shear_int, Re_bulk_int, divergence_v) + do l_loop = iz%beg, iz%end + do k_loop = iy%beg, iy%end + do j_loop = ix%beg, ix%end idx_right_phys(1) = j_loop idx_right_phys(2) = k_loop @@ -4554,31 +4173,28 @@ contains divergence_v = divergence_v + vel_grad_avg(i_dim, i_dim) end do - vel_src_at_interface = 0.0_wp - if (norm_dir == 1) then - Re_shear = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 1) - Re_bulk = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 2) - do i_dim = 1, num_dims - vel_src_at_interface(i_dim) = vel_src_rsx_vf(j_loop, k_loop, l_loop, i_dim) - end do - else if (norm_dir == 2) then - Re_shear = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 1) - Re_bulk = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 2) - do i_dim = 1, num_dims - vel_src_at_interface(i_dim) = vel_src_rsy_vf(k_loop, j_loop, l_loop, i_dim) - end do - else - Re_shear = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 1) - Re_bulk = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 2) - do i_dim = 1, num_dims - vel_src_at_interface(i_dim) = vel_src_rsz_vf(l_loop, k_loop, j_loop, i_dim) - end do - end if + select case (norm_dir) + case (1) + Re_shear_int = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 1) + Re_bulk_int = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 2) + vel_src_at_interface(1:num_dims) = vel_src_rsx_vf(j_loop, k_loop, l_loop, 1:num_dims) + case (2) + Re_shear_int = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 1) + Re_bulk_int = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 2) + vel_src_at_interface(1:num_dims) = vel_src_rsy_vf(k_loop, j_loop, l_loop, 1:num_dims) + case (3) + Re_shear_int = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 1) + Re_bulk_int = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 2) + vel_src_at_interface(1:num_dims) = vel_src_rsz_vf(l_loop, k_loop, j_loop, 1:num_dims) + case default + cycle + end select + + current_tau_shear = 0.0_wp + current_tau_bulk = 0.0_wp if (shear_stress) then - current_tau_shear = 0.0_wp - call s_calculate_shear_stress_tensor(vel_grad_avg, Re_shear, divergence_v, current_tau_shear) - + call s_calculate_shear_stress_tensor(vel_grad_avg, Re_shear_int, divergence_v, current_tau_shear) do i_dim = 1, num_dims flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) = & flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) - current_tau_shear(norm_dir, i_dim) @@ -4590,9 +4206,7 @@ contains end if if (bulk_stress) then - current_tau_bulk = 0.0_wp - call s_calculate_bulk_stress_tensor(Re_bulk, divergence_v, current_tau_bulk) - + call s_calculate_bulk_stress_tensor(Re_bulk_int, divergence_v, current_tau_bulk) do i_dim = 1, num_dims flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) = & flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) - current_tau_bulk(norm_dir, i_dim) @@ -4610,16 +4224,15 @@ contains end subroutine s_compute_cartesian_viscous_source_flux - !> @brief Calculates shear stress tensor components. - !! tau_ij_shear = ( (dui/dxj + duj/dxi) - (2/3)*(div_v)*delta_ij ) / Re_shear - !! @param[in] vel_grad_avg Averaged velocity gradient tensor (d(vel_i)/d(coord_j)). + !> @brief Calculates Cartesian shear stress tensor components. + !! @param[in] vel_grad_avg Averaged velocity gradient tensor $L_{ij} = \partial v_i/\partial x_j$. !! @param[in] Re_shear Shear Reynolds number. - !! @param[in] divergence_v Velocity divergence (du/dx + dv/dy + dw/dz). - !! @param[out] tau_shear_out Calculated shear stress tensor (stress on i-face, j-direction). + !! @param[in] divergence_v Velocity divergence $\nabla \cdot \mathbf{v}$. + !! @param[out] tau_shear_out Calculated shear stress tensor $\tau_{ij}$. subroutine s_calculate_shear_stress_tensor(vel_grad_avg, Re_shear, divergence_v, tau_shear_out) - use m_global_parameters, only: wp, num_dims, sgm_eps - implicit none + !$acc routine seq + implicit none ! Arguments real(wp), dimension(num_dims, num_dims), intent(in) :: vel_grad_avg real(wp), intent(in) :: Re_shear @@ -4627,34 +4240,31 @@ contains real(wp), dimension(num_dims, num_dims), intent(out) :: tau_shear_out ! Local variables - integer :: i_dim !< Loop iterator for face normal. - integer :: j_dim !< Loop iterator for force component direction. + integer :: i_face_norm !!< Tensor row index $i$ (face normal). + integer :: j_force_dir !!< Tensor col index $j$ (force component direction). tau_shear_out = 0.0_wp - if (abs(Re_shear) < sgm_eps) then - return - end if - - do i_dim = 1, num_dims - do j_dim = 1, num_dims - tau_shear_out(i_dim, j_dim) = (vel_grad_avg(j_dim, i_dim) + vel_grad_avg(i_dim, j_dim))/Re_shear - if (i_dim == j_dim) then - tau_shear_out(i_dim, j_dim) = tau_shear_out(i_dim, j_dim) - & - (2.0_wp/3.0_wp)*divergence_v/Re_shear + do i_face_norm = 1, num_dims + do j_force_dir = 1, num_dims + tau_shear_out(i_face_norm, j_force_dir) = (vel_grad_avg(j_force_dir, i_face_norm) + & + vel_grad_avg(i_face_norm, j_force_dir))/Re_shear + if (i_face_norm == j_force_dir) then + tau_shear_out(i_face_norm, j_force_dir) = tau_shear_out(i_face_norm, j_force_dir) - & + (2.0_wp/3.0_wp)*divergence_v/Re_shear end if end do end do end subroutine s_calculate_shear_stress_tensor - !> @brief Calculates bulk stress tensor components (diagonal only). - !! tau_ii_bulk = (div_v) / Re_bulk. Off-diagonals are zero. + !> @brief Calculates Cartesian bulk stress tensor components (diagonal only). !! @param[in] Re_bulk Bulk Reynolds number. - !! @param[in] divergence_v Velocity divergence (du/dx + dv/dy + dw/dz). - !! @param[out] tau_bulk_out Calculated bulk stress tensor (stress on i-face, i-direction). + !! @param[in] divergence_v Velocity divergence $\nabla \cdot \mathbf{v}$. + !! @param[out] tau_bulk_out Calculated bulk stress tensor $\tau_{ii}$. subroutine s_calculate_bulk_stress_tensor(Re_bulk, divergence_v, tau_bulk_out) - use m_global_parameters, only: wp, num_dims, sgm_eps + !$acc routine seq + implicit none ! Arguments @@ -4663,14 +4273,10 @@ contains real(wp), dimension(num_dims, num_dims), intent(out) :: tau_bulk_out ! Local variables - integer :: i_dim !< Loop iterator for diagonal components. + integer :: i_dim !!< Loop iterator for diagonal components. tau_bulk_out = 0.0_wp - if (abs(Re_bulk) < sgm_eps) then - return - end if - do i_dim = 1, num_dims tau_bulk_out(i_dim, i_dim) = divergence_v/Re_bulk end do From 77c7c28bebf2ceab1d7d1ab3b0920629ab92967d Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 27 May 2025 15:50:27 -0400 Subject: [PATCH 3/8] fix --- src/simulation/m_riemann_solvers.fpp | 164 +++++++++++++++------------ 1 file changed, 93 insertions(+), 71 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index b6376fe72a..8524a5ea08 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -4101,20 +4101,33 @@ contains end subroutine s_compute_cylindrical_viscous_source_flux - !> @brief Computes Cartesian viscous source flux contributions. - !! Calls helpers for shear/bulk stress after averaging gradients. - !! @param[in] velL_vf,velR_vf Left/Right boundary velocities. - !! @param[in] dvelL_dx_vf,dvelR_dx_vf Left/Right boundary $\partial v_i/\partial x$. - !! @param[in] dvelL_dy_vf,dvelR_dy_vf Left/Right boundary $\partial v_i/\partial y$. - !! @param[in] dvelL_dz_vf,dvelR_dz_vf Left/Right boundary $\partial v_i/\partial z$. - !! @param[inout] flux_src_vf Intercell source flux array to update. + !> @brief Computes Cartesian viscous source flux contributions for momentum and energy. + !! Calculates averaged velocity gradients, gets Re and interface velocities, + !! calls helpers for shear/bulk stress, then updates `flux_src_vf`. + !! @param[in] velL_vf Left boundary velocity (num_dims scalar_field). + !! @param[in] dvelL_dx_vf Left boundary d(vel)/dx (num_dims scalar_field). + !! @param[in] dvelL_dy_vf Left boundary d(vel)/dy (num_dims scalar_field). + !! @param[in] dvelL_dz_vf Left boundary d(vel)/dz (num_dims scalar_field). + !! @param[in] velR_vf Right boundary velocity (num_dims scalar_field). + !! @param[in] dvelR_dx_vf Right boundary d(vel)/dx (num_dims scalar_field). + !! @param[in] dvelR_dy_vf Right boundary d(vel)/dy (num_dims scalar_field). + !! @param[in] dvelR_dz_vf Right boundary d(vel)/dz (num_dims scalar_field). + !! @param[inout] flux_src_vf Intercell source flux array to update (sys_size scalar_field). !! @param[in] norm_dir Interface normal direction (1=x, 2=y, 3=z). - !! @param[in] ix,iy,iz Global X,Y,Z-direction loop bounds. + !! @param[in] ix X-direction loop bounds (int_bounds_info). + !! @param[in] iy Y-direction loop bounds (int_bounds_info). + !! @param[in] iz Z-direction loop bounds (int_bounds_info). subroutine s_compute_cartesian_viscous_source_flux(velL_vf, & - dvelL_dx_vf, dvelL_dy_vf, dvelL_dz_vf, & + dvelL_dx_vf, & + dvelL_dy_vf, & + dvelL_dz_vf, & velR_vf, & - dvelR_dx_vf, dvelR_dy_vf, dvelR_dz_vf, & - flux_src_vf, norm_dir, ix, iy, iz) + dvelR_dx_vf, & + dvelR_dy_vf, & + dvelR_dz_vf, & + flux_src_vf, & + norm_dir, & + ix, iy, iz) ! Arguments type(scalar_field), dimension(num_dims), intent(in) :: velL_vf, velR_vf @@ -4126,28 +4139,30 @@ contains type(int_bounds_info), intent(in) :: ix, iy, iz ! Local variables - real(wp), dimension(num_dims, num_dims) :: vel_grad_avg !!< Averaged velocity gradient tensor $L_{ij} = \partial v_i/\partial x_j$. - real(wp), dimension(num_dims, num_dims) :: current_tau_shear !!< Current shear stress tensor $\tau_{ij}^{(s)}$. - real(wp), dimension(num_dims, num_dims) :: current_tau_bulk !!< Current bulk stress tensor $\tau_{ij}^{(b)}$ (diagonal). - real(wp), dimension(num_dims) :: vel_src_at_interface !!< Interface velocities $(v_x,v_y,v_z)$ for viscous work. - integer, dimension(3) :: idx_right_phys !!< Physical $(j,k,l)$ indices for 'right' state. - - real(wp) :: Re_shear_int !!< Interface shear Reynolds number. - real(wp) :: Re_bulk_int !!< 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 iterator for flux updates. - integer :: vel_comp_idx !!< Iterator for velocity components. - real(wp) :: divergence_v !!< Velocity divergence $\nabla \cdot \mathbf{v}$ at interface. + real(wp), dimension(num_dims, num_dims) :: vel_grad_avg !< Averaged velocity gradient tensor `d(vel_i)/d(coord_j)`. + real(wp), dimension(num_dims, num_dims) :: current_tau_shear !< Current shear stress tensor. + 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. + 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. !$acc parallel loop collapse(3) gang vector default(present) & - !$acc private(idx_right_phys, vel_grad_avg, current_tau_shear, current_tau_bulk, & - !$acc vel_src_at_interface, Re_shear_int, Re_bulk_int, divergence_v) - do l_loop = iz%beg, iz%end - do k_loop = iy%beg, iy%end - do j_loop = ix%beg, ix%end + !$acc private(idx_right_phys, vel_grad_avg, & + !$acc current_tau_shear, current_tau_bulk, vel_src_at_interface, & + !$acc Re_shear, Re_bulk, divergence_v, i_dim, vel_comp_idx) + do l_loop = isz%beg, isz%end + do k_loop = isy%beg, isy%end + do j_loop = isx%beg, isx%end idx_right_phys(1) = j_loop idx_right_phys(2) = k_loop @@ -4173,28 +4188,31 @@ contains divergence_v = divergence_v + vel_grad_avg(i_dim, i_dim) end do - select case (norm_dir) - case (1) - Re_shear_int = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 1) - Re_bulk_int = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 2) - vel_src_at_interface(1:num_dims) = vel_src_rsx_vf(j_loop, k_loop, l_loop, 1:num_dims) - case (2) - Re_shear_int = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 1) - Re_bulk_int = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 2) - vel_src_at_interface(1:num_dims) = vel_src_rsy_vf(k_loop, j_loop, l_loop, 1:num_dims) - case (3) - Re_shear_int = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 1) - Re_bulk_int = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 2) - vel_src_at_interface(1:num_dims) = vel_src_rsz_vf(l_loop, k_loop, j_loop, 1:num_dims) - case default - cycle - end select - - current_tau_shear = 0.0_wp - current_tau_bulk = 0.0_wp + vel_src_at_interface = 0.0_wp + if (norm_dir == 1) then + Re_shear = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 1) + Re_bulk = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 2) + do i_dim = 1, num_dims + vel_src_at_interface(i_dim) = vel_src_rsx_vf(j_loop, k_loop, l_loop, i_dim) + end do + else if (norm_dir == 2) then + Re_shear = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 1) + Re_bulk = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 2) + do i_dim = 1, num_dims + vel_src_at_interface(i_dim) = vel_src_rsy_vf(k_loop, j_loop, l_loop, i_dim) + end do + else + Re_shear = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 1) + Re_bulk = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 2) + do i_dim = 1, num_dims + vel_src_at_interface(i_dim) = vel_src_rsz_vf(l_loop, k_loop, j_loop, i_dim) + end do + end if if (shear_stress) then - call s_calculate_shear_stress_tensor(vel_grad_avg, Re_shear_int, divergence_v, current_tau_shear) + current_tau_shear = 0.0_wp + call s_calculate_shear_stress_tensor(vel_grad_avg, Re_shear, divergence_v, current_tau_shear) + do i_dim = 1, num_dims flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) = & flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) - current_tau_shear(norm_dir, i_dim) @@ -4206,7 +4224,9 @@ contains end if if (bulk_stress) then - call s_calculate_bulk_stress_tensor(Re_bulk_int, divergence_v, current_tau_bulk) + current_tau_bulk = 0.0_wp + call s_calculate_bulk_stress_tensor(Re_bulk, divergence_v, current_tau_bulk) + do i_dim = 1, num_dims flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) = & flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) - current_tau_bulk(norm_dir, i_dim) @@ -4224,15 +4244,17 @@ contains end subroutine s_compute_cartesian_viscous_source_flux - !> @brief Calculates Cartesian shear stress tensor components. - !! @param[in] vel_grad_avg Averaged velocity gradient tensor $L_{ij} = \partial v_i/\partial x_j$. + !> @brief Calculates shear stress tensor components. + !! tau_ij_shear = ( (dui/dxj + duj/dxi) - (2/3)*(div_v)*delta_ij ) / Re_shear + !! @param[in] vel_grad_avg Averaged velocity gradient tensor (d(vel_i)/d(coord_j)). !! @param[in] Re_shear Shear Reynolds number. - !! @param[in] divergence_v Velocity divergence $\nabla \cdot \mathbf{v}$. - !! @param[out] tau_shear_out Calculated shear stress tensor $\tau_{ij}$. + !! @param[in] divergence_v Velocity divergence (du/dx + dv/dy + dw/dz). + !! @param[out] tau_shear_out Calculated shear stress tensor (stress on i-face, j-direction). subroutine s_calculate_shear_stress_tensor(vel_grad_avg, Re_shear, divergence_v, tau_shear_out) - !$acc routine seq + !$acc routine seq implicit none + ! Arguments real(wp), dimension(num_dims, num_dims), intent(in) :: vel_grad_avg real(wp), intent(in) :: Re_shear @@ -4240,30 +4262,30 @@ contains real(wp), dimension(num_dims, num_dims), intent(out) :: tau_shear_out ! Local variables - integer :: i_face_norm !!< Tensor row index $i$ (face normal). - integer :: j_force_dir !!< Tensor col index $j$ (force component direction). + integer :: i_dim !< Loop iterator for face normal. + integer :: j_dim !< Loop iterator for force component direction. tau_shear_out = 0.0_wp - do i_face_norm = 1, num_dims - do j_force_dir = 1, num_dims - tau_shear_out(i_face_norm, j_force_dir) = (vel_grad_avg(j_force_dir, i_face_norm) + & - vel_grad_avg(i_face_norm, j_force_dir))/Re_shear - if (i_face_norm == j_force_dir) then - tau_shear_out(i_face_norm, j_force_dir) = tau_shear_out(i_face_norm, j_force_dir) - & - (2.0_wp/3.0_wp)*divergence_v/Re_shear + do i_dim = 1, num_dims + do j_dim = 1, num_dims + tau_shear_out(i_dim, j_dim) = (vel_grad_avg(j_dim, i_dim) + vel_grad_avg(i_dim, j_dim))/Re_shear + if (i_dim == j_dim) then + tau_shear_out(i_dim, j_dim) = tau_shear_out(i_dim, j_dim) - & + (2.0_wp/3.0_wp)*divergence_v/Re_shear end if end do end do end subroutine s_calculate_shear_stress_tensor - !> @brief Calculates Cartesian bulk stress tensor components (diagonal only). + !> @brief Calculates bulk stress tensor components (diagonal only). + !! tau_ii_bulk = (div_v) / Re_bulk. Off-diagonals are zero. !! @param[in] Re_bulk Bulk Reynolds number. - !! @param[in] divergence_v Velocity divergence $\nabla \cdot \mathbf{v}$. - !! @param[out] tau_bulk_out Calculated bulk stress tensor $\tau_{ii}$. + !! @param[in] divergence_v Velocity divergence (du/dx + dv/dy + dw/dz). + !! @param[out] tau_bulk_out Calculated bulk stress tensor (stress on i-face, i-direction). subroutine s_calculate_bulk_stress_tensor(Re_bulk, divergence_v, tau_bulk_out) - !$acc routine seq + !$acc routine seq implicit none @@ -4273,7 +4295,7 @@ contains real(wp), dimension(num_dims, num_dims), intent(out) :: tau_bulk_out ! Local variables - integer :: i_dim !!< Loop iterator for diagonal components. + integer :: i_dim !< Loop iterator for diagonal components. tau_bulk_out = 0.0_wp From 7b657ed869ce35abc2dd4ead03284f064ea7fcd1 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 27 May 2025 15:54:17 -0400 Subject: [PATCH 4/8] pretty; --- src/simulation/m_riemann_solvers.fpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 8524a5ea08..2d257eaf8f 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -4251,7 +4251,7 @@ contains !! @param[in] divergence_v Velocity divergence (du/dx + dv/dy + dw/dz). !! @param[out] tau_shear_out Calculated shear stress tensor (stress on i-face, j-direction). subroutine s_calculate_shear_stress_tensor(vel_grad_avg, Re_shear, divergence_v, tau_shear_out) - !$acc routine seq + !$acc routine seq implicit none @@ -4285,7 +4285,7 @@ contains !! @param[in] divergence_v Velocity divergence (du/dx + dv/dy + dw/dz). !! @param[out] tau_bulk_out Calculated bulk stress tensor (stress on i-face, i-direction). subroutine s_calculate_bulk_stress_tensor(Re_bulk, divergence_v, tau_bulk_out) - !$acc routine seq + !$acc routine seq implicit none From e63dc1e1676486f9e3e92bc7a4d363651fb0ea4c Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 27 May 2025 21:55:44 -0400 Subject: [PATCH 5/8] make it work --- src/simulation/m_riemann_solvers.fpp | 435 +++++++++++++++++++-------- 1 file changed, 303 insertions(+), 132 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 2d257eaf8f..8c1622ecbc 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -3940,164 +3940,338 @@ contains end subroutine s_initialize_riemann_solver - !> @brief Computes cylindrical viscous source flux contributions for momentum and energy. - !! Calculates Cartesian components of the stress tensor using averaged velocity derivatives - !! and cylindrical geometric factors, then updates `flux_src_vf`. - !! Assumes x-dir is axial (z_cyl), y-dir is radial (r_cyl), z-dir is azimuthal (theta_cyl for derivatives). - !! @param[in] velL_vf Left boundary velocity ($v_x, v_y, v_z$) (num_dims scalar_field). - !! @param[in] dvelL_dx_vf Left boundary $\partial v_i/\partial x$ (num_dims scalar_field). - !! @param[in] dvelL_dy_vf Left boundary $\partial v_i/\partial y$ (num_dims scalar_field). - !! @param[in] dvelL_dz_vf Left boundary $\partial v_i/\partial z$ (num_dims scalar_field). - !! @param[in] velR_vf Right boundary velocity ($v_x, v_y, v_z$) (num_dims scalar_field). - !! @param[in] dvelR_dx_vf Right boundary $\partial v_i/\partial x$ (num_dims scalar_field). - !! @param[in] dvelR_dy_vf Right boundary $\partial v_i/\partial y$ (num_dims scalar_field). - !! @param[in] dvelR_dz_vf Right boundary $\partial v_i/\partial z$ (num_dims scalar_field). - !! @param[inout] flux_src_vf Intercell source flux array to update (sys_size scalar_field). - !! @param[in] norm_dir Interface normal direction (1=x-face, 2=y-face, 3=z-face). - !! @param[in] ix Global X-direction loop bounds (int_bounds_info). - !! @param[in] iy Global Y-direction loop bounds (int_bounds_info). - !! @param[in] iz Global Z-direction loop bounds (int_bounds_info). + !> @brief Computes viscous source fluxes for cylindrical coordinates. + !! This subroutine evaluates the contribution of viscous stresses (both shear and bulk) + !! to the source terms for momentum and energy equations in a cylindrical coordinate system. + !! It handles calculations based on the normal direction of the flux (axial, radial, or azimuthal). + !! + !! @param velL_vf Left WENO-reconstructed cell-boundary values of the velocity components. + !! @param dvelL_dx_vf Left WENO-reconstructed cell-boundary values of the x-derivative of velocity components. + !! @param dvelL_dy_vf Left WENO-reconstructed cell-boundary values of the y-derivative of velocity components. + !! @param dvelL_dz_vf Left WENO-reconstructed cell-boundary values of the z-derivative of velocity components. + !! @param velR_vf Right WENO-reconstructed cell-boundary values of the velocity components. + !! @param dvelR_dx_vf Right WENO-reconstructed cell-boundary values of the x-derivative of velocity components. + !! @param dvelR_dy_vf Right WENO-reconstructed cell-boundary values of the y-derivative of velocity components. + !! @param dvelR_dz_vf Right WENO-reconstructed cell-boundary values of the z-derivative of velocity components. + !! @param flux_src_vf Array to which the computed viscous source fluxes are added. + !! @param norm_dir Integer indicating the normal direction for flux calculation (1 for axial, 2 for radial, 3 for azimuthal). + !! @param ix Index bounds in the first coordinate direction for loops. + !! @param iy Index bounds in the second coordinate direction for loops. + !! @param iz Index bounds in the third coordinate direction for loops. subroutine s_compute_cylindrical_viscous_source_flux(velL_vf, & - dvelL_dx_vf, dvelL_dy_vf, dvelL_dz_vf, & + dvelL_dx_vf, & + dvelL_dy_vf, & + dvelL_dz_vf, & velR_vf, & - dvelR_dx_vf, dvelR_dy_vf, dvelR_dz_vf, & - flux_src_vf, norm_dir, ix, iy, iz) + dvelR_dx_vf, & + dvelR_dy_vf, & + dvelR_dz_vf, & + flux_src_vf, & + norm_dir, & + ix, iy, iz) + + implicit none type(scalar_field), dimension(num_dims), intent(in) :: velL_vf, velR_vf type(scalar_field), dimension(num_dims), intent(in) :: dvelL_dx_vf, dvelR_dx_vf type(scalar_field), dimension(num_dims), intent(in) :: dvelL_dy_vf, dvelR_dy_vf type(scalar_field), dimension(num_dims), intent(in) :: dvelL_dz_vf, dvelR_dz_vf + type(scalar_field), dimension(sys_size), intent(inout) :: flux_src_vf + integer, intent(in) :: norm_dir type(int_bounds_info), intent(in) :: ix, iy, iz - ! Local variables - real(wp), dimension(num_dims) :: avg_v_int !!< Averaged interface velocity $(v_x, v_y, v_z)$ (grid directions). - real(wp), dimension(num_dims) :: avg_dvdx_int !!< Averaged interface $\partial v_i/\partial x$ (grid dir 1). - real(wp), dimension(num_dims) :: avg_dvdy_int !!< Averaged interface $\partial v_i/\partial y$ (grid dir 2). - real(wp), dimension(num_dims) :: avg_dvdz_int !!< Averaged interface $\partial v_i/\partial z$ (grid dir 3). + real(wp), dimension(num_dims) :: avg_vel + real(wp), dimension(num_dims) :: dvel_avg_dx + real(wp), dimension(num_dims) :: dvel_avg_dy + real(wp), dimension(num_dims) :: dvel_avg_dz - real(wp), dimension(num_dims) :: stress_vector_shear !!< Shear stress vector $(\sigma_{N1}, \sigma_{N2}, \sigma_{N3})$ on N-face (grid directions). - real(wp) :: stress_normal_bulk !!< Normal bulk stress component $\sigma_{NN}$ on N-face. + real(wp), dimension(num_dims, num_dims) :: tau_Re ! Viscous stress tensor components - real(wp) :: Re_s, Re_b !!< Effective interface shear and bulk Reynolds numbers. - real(wp), dimension(num_dims) :: vel_src_int !!< Interface velocity $(v_1,v_2,v_3)$ (grid directions) for viscous work. - real(wp) :: r_eff !!< Effective radius at interface for cylindrical terms. - real(wp) :: div_v_term_const !!< Common term $-(2/3)(\nabla \cdot \mathbf{v}) / \text{Re}_s$ for shear stress diagonal. - real(wp) :: divergence_cyl !!< Full divergence $\nabla \cdot \mathbf{v}$ in cylindrical coordinates. + integer :: i_loop ! Loop iterator for components + integer :: j_idx, k_idx, l_idx ! Spatial loop iterators - integer :: j, k, l !!< Loop iterators for $x, y, z$ grid directions. - integer :: i_vel !!< Loop iterator for velocity components. - integer :: idx_rp(3) !!< Indices $(j,k,l)$ of 'right' point for averaging. + ! norm_dir = 1: Axial direction flux (e.g., z-flux if grid x-axis is solver's z-axis) + ! In this block, for derivatives: dx is axial, dy is radial, dz is azimuthal component of gradient. + if (norm_dir == 1) then + ! Part 1: Axial derivative contributions (1D, 2D, 3D cylindrical) + !$acc parallel loop collapse(3) gang vector default(present) & + !$acc private(dvel_avg_dx, tau_Re) + do l_idx = iz%beg, iz%end + do k_idx = iy%beg, iy%end + do j_idx = ix%beg, ix%end + + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j_idx, k_idx, l_idx) + & + dvelR_dx_vf(1)%sf(j_idx + 1, k_idx, l_idx)) + + if (shear_stress) then + ! tau_zz_shear (normal-normal stress component in axial direction) + tau_Re(1, 1) = (4._wp/3._wp)*dvel_avg_dx(1)/Re_avg_rsx_vf(j_idx, k_idx, l_idx, 1) + flux_src_vf(momxb)%sf(j_idx, k_idx, l_idx) = flux_src_vf(momxb)%sf(j_idx, k_idx, l_idx) - tau_Re(1, 1) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsx_vf(j_idx, k_idx, l_idx, 1)*tau_Re(1, 1) + end if - !$acc parallel loop collapse(3) gang vector default(present) & - !$acc private(idx_rp, avg_v_int, avg_dvdx_int, avg_dvdy_int, avg_dvdz_int, & - !$acc Re_s, Re_b, vel_src_int, r_eff, divergence_cyl, & - !$acc stress_vector_shear, stress_normal_bulk, div_v_term_const) - do l = iz%beg, iz%end - do k = iy%beg, iy%end - do j = ix%beg, ix%end - - ! Determine indices for the 'right' state for averaging across the interface - idx_rp = [j, k, l] - idx_rp(norm_dir) = idx_rp(norm_dir) + 1 - - ! Average velocities and their derivatives at the interface - ! For cylindrical: x-dir ~ axial (z_cyl), y-dir ~ radial (r_cyl), z-dir ~ azimuthal (theta_cyl) - do i_vel = 1, num_dims - avg_v_int(i_vel) = 0.5_wp*(velL_vf(i_vel)%sf(j, k, l) + velR_vf(i_vel)%sf(idx_rp(1), idx_rp(2), idx_rp(3))) - - avg_dvdx_int(i_vel) = 0.5_wp*(dvelL_dx_vf(i_vel)%sf(j, k, l) + & - dvelR_dx_vf(i_vel)%sf(idx_rp(1), idx_rp(2), idx_rp(3))) - if (num_dims > 1) then - avg_dvdy_int(i_vel) = 0.5_wp*(dvelL_dy_vf(i_vel)%sf(j, k, l) + & - dvelR_dy_vf(i_vel)%sf(idx_rp(1), idx_rp(2), idx_rp(3))) - else - avg_dvdy_int(i_vel) = 0.0_wp + if (bulk_stress) then + ! tau_zz_bulk + tau_Re(1, 1) = dvel_avg_dx(1)/Re_avg_rsx_vf(j_idx, k_idx, l_idx, 2) + flux_src_vf(momxb)%sf(j_idx, k_idx, l_idx) = flux_src_vf(momxb)%sf(j_idx, k_idx, l_idx) - tau_Re(1, 1) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsx_vf(j_idx, k_idx, l_idx, 1)*tau_Re(1, 1) end if - if (num_dims > 2) then - avg_dvdz_int(i_vel) = 0.5_wp*(dvelL_dz_vf(i_vel)%sf(j, k, l) + & - dvelR_dz_vf(i_vel)%sf(idx_rp(1), idx_rp(2), idx_rp(3))) - else - avg_dvdz_int(i_vel) = 0.0_wp + end do + end do + end do + + if (n == 0) return ! No radial extent + + ! Part 2: Radial derivative contributions to axial flux (2D, 3D cylindrical) + !$acc parallel loop collapse(3) gang vector default(present) & + !$acc private(avg_vel, dvel_avg_dx, dvel_avg_dy, tau_Re, i_loop) + do l_idx = iz%beg, iz%end + do k_idx = iy%beg, iy%end + do j_idx = ix%beg, ix%end + + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j_idx, k_idx, l_idx) + velR_vf(2)%sf(j_idx + 1, k_idx, l_idx)) ! v_r + + do i_loop = 1, 2 ! u_axial (1), u_radial (2) + dvel_avg_dy(i_loop) = 0.5_wp*(dvelL_dy_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & + dvelR_dy_vf(i_loop)%sf(j_idx + 1, k_idx, l_idx)) ! du_i/dr + end do + dvel_avg_dx(2) = 0.5_wp*(dvelL_dx_vf(2)%sf(j_idx, k_idx, l_idx) + & + dvelR_dx_vf(2)%sf(j_idx + 1, k_idx, l_idx)) ! du_r/dx_axial (dx_axial is normal here) + + if (shear_stress) then + ! tau_zz component from radial terms: -(2/3)mu (dv_r/dr + v_r/r) + tau_Re(1, 1) = -(2._wp/3._wp)*(dvel_avg_dy(2) + avg_vel(2)/y_cc(k_idx))/ & + Re_avg_rsx_vf(j_idx, k_idx, l_idx, 1) + ! tau_zr shear stress: mu (dv_axial/dr + dv_r/dx_axial) + tau_Re(1, 2) = (dvel_avg_dy(1) + dvel_avg_dx(2))/Re_avg_rsx_vf(j_idx, k_idx, l_idx, 1) + + do i_loop = 1, 2 ! Axial (1) and radial (2) momentum components + flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) = & + flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) - tau_Re(1, i_loop) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = & + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsx_vf(j_idx, k_idx, l_idx, i_loop)*tau_Re(1, i_loop) + end do + end if + + if (bulk_stress) then + ! tau_zz_bulk from radial terms: mu_bulk (dv_r/dr + v_r/r) + tau_Re(1, 1) = (dvel_avg_dy(2) + avg_vel(2)/y_cc(k_idx))/Re_avg_rsx_vf(j_idx, k_idx, l_idx, 2) + flux_src_vf(momxb)%sf(j_idx, k_idx, l_idx) = flux_src_vf(momxb)%sf(j_idx, k_idx, l_idx) - tau_Re(1, 1) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsx_vf(j_idx, k_idx, l_idx, 1)*tau_Re(1, 1) end if end do + end do + end do - ! Get Re numbers and interface velocity for viscous work - select case (norm_dir) - case (1) ! x-face (axial face in z_cyl direction) - Re_s = Re_avg_rsx_vf(j, k, l, 1) - Re_b = Re_avg_rsx_vf(j, k, l, 2) - vel_src_int = vel_src_rsx_vf(j, k, l, 1:num_dims) - r_eff = y_cc(k) - case (2) ! y-face (radial face in r_cyl direction) - Re_s = Re_avg_rsy_vf(k, j, l, 1) - Re_b = Re_avg_rsy_vf(k, j, l, 2) - vel_src_int = vel_src_rsy_vf(k, j, l, 1:num_dims) - r_eff = y_cb(k) - case (3) ! z-face (azimuthal face in theta_cyl direction) - Re_s = Re_avg_rsz_vf(l, k, j, 1) - Re_b = Re_avg_rsz_vf(l, k, j, 2) - vel_src_int = vel_src_rsz_vf(l, k, j, 1:num_dims) - r_eff = y_cc(k) - end select - - ! Divergence in cylindrical coordinates (vx=vz_cyl, vy=vr_cyl, vz=vtheta_cyl) - divergence_cyl = avg_dvdx_int(1) + avg_dvdy_int(2) + avg_v_int(2)/r_eff - if (num_dims > 2) then - divergence_cyl = divergence_cyl + avg_dvdz_int(3)/r_eff - end if + if (p == 0) return ! No azimuthal extent - stress_vector_shear = 0.0_wp - stress_normal_bulk = 0.0_wp + ! Part 3: Azimuthal derivative contributions to axial flux (3D cylindrical) + !$acc parallel loop collapse(3) gang vector default(present) & + !$acc private(avg_vel, dvel_avg_dx, dvel_avg_dz, tau_Re, i_loop) + do l_idx = iz%beg, iz%end + do k_idx = iy%beg, iy%end + do j_idx = ix%beg, ix%end - if (shear_stress) then - div_v_term_const = -(2.0_wp/3.0_wp)*divergence_cyl/Re_s + do i_loop = 1, 3, 2 ! u_axial (1), u_azimuthal (3) + dvel_avg_dz(i_loop) = 0.5_wp*(dvelL_dz_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & + dvelR_dz_vf(i_loop)%sf(j_idx + 1, k_idx, l_idx)) ! du_i/(r dtheta) + end do + dvel_avg_dx(3) = 0.5_wp*(dvelL_dx_vf(3)%sf(j_idx, k_idx, l_idx) + & + dvelR_dx_vf(3)%sf(j_idx + 1, k_idx, l_idx)) ! du_theta/dx_axial + + if (shear_stress) then + ! tau_zz component from azimuthal terms: -(2/3)mu (1/r dv_theta/dtheta) + tau_Re(1, 1) = -(2._wp/3._wp)*(dvel_avg_dz(3)/y_cc(k_idx))/Re_avg_rsx_vf(j_idx, k_idx, l_idx, 1) + ! tau_z_theta shear stress: mu ( (1/r)dv_axial/dtheta + dv_theta/dx_axial ) + tau_Re(1, 3) = (dvel_avg_dz(1)/y_cc(k_idx) + dvel_avg_dx(3))/Re_avg_rsx_vf(j_idx, k_idx, l_idx, 1) + + do i_loop = 1, 3, 2 ! Axial (1) and azimuthal (3) momentum components + flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) = & + flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) - tau_Re(1, i_loop) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = & + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsx_vf(j_idx, k_idx, l_idx, i_loop)*tau_Re(1, i_loop) + end do + end if - select case (norm_dir) - case (1) ! X-face (axial normal, z_cyl) - stress_vector_shear(1) = (2.0_wp*avg_dvdx_int(1))/Re_s + div_v_term_const - if (num_dims > 1) then - stress_vector_shear(2) = (avg_dvdy_int(1) + avg_dvdx_int(2))/Re_s - end if - if (num_dims > 2) then - stress_vector_shear(3) = (avg_dvdz_int(1)/r_eff + avg_dvdx_int(3))/Re_s - end if - case (2) ! Y-face (radial normal, r_cyl) - if (num_dims > 1) then - stress_vector_shear(1) = (avg_dvdy_int(1) + avg_dvdx_int(2))/Re_s - stress_vector_shear(2) = (2.0_wp*avg_dvdy_int(2))/Re_s + div_v_term_const - if (num_dims > 2) then - stress_vector_shear(3) = (avg_dvdz_int(2)/r_eff - avg_v_int(3)/r_eff + avg_dvdy_int(3))/Re_s - end if - else - stress_vector_shear(1) = (2.0_wp*avg_dvdx_int(1))/Re_s + div_v_term_const - end if - case (3) ! Z-face (azimuthal normal, theta_cyl) - if (num_dims > 2) then - stress_vector_shear(1) = (avg_dvdz_int(1)/r_eff + avg_dvdx_int(3))/Re_s - stress_vector_shear(2) = (avg_dvdz_int(2)/r_eff - avg_v_int(3)/r_eff + avg_dvdy_int(3))/Re_s - stress_vector_shear(3) = (2.0_wp*(avg_dvdz_int(3)/r_eff + avg_v_int(2)/r_eff))/Re_s + div_v_term_const - end if - end select + if (bulk_stress) then + ! tau_zz_bulk from azimuthal terms: mu_bulk (1/r dv_theta/dtheta) + tau_Re(1, 1) = (dvel_avg_dz(3)/y_cc(k_idx))/Re_avg_rsx_vf(j_idx, k_idx, l_idx, 2) + flux_src_vf(momxb)%sf(j_idx, k_idx, l_idx) = flux_src_vf(momxb)%sf(j_idx, k_idx, l_idx) - tau_Re(1, 1) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsx_vf(j_idx, k_idx, l_idx, 1)*tau_Re(1, 1) + end if + end do + end do + end do - do i_vel = 1, num_dims - flux_src_vf(momxb + i_vel - 1)%sf(j, k, l) = flux_src_vf(momxb + i_vel - 1)%sf(j, k, l) - stress_vector_shear(i_vel) - flux_src_vf(E_idx)%sf(j, k, l) = flux_src_vf(E_idx)%sf(j, k, l) - vel_src_int(i_vel)*stress_vector_shear(i_vel) + ! norm_dir = 2: Radial direction flux (e.g., r-flux if grid y-axis is solver's r-axis) + ! In this block, for derivatives: dy is radial (normal), dx is axial, dz is azimuthal component of gradient. + elseif (norm_dir == 2) then + ! Part 1: Axial and radial derivative contributions to radial flux + !$acc parallel loop collapse(3) gang vector default(present) & + !$acc private(avg_vel, dvel_avg_dx, dvel_avg_dy, tau_Re, i_loop) + do l_idx = iz%beg, iz%end + do k_idx = iy%beg, iy%end ! Radial flux index + do j_idx = ix%beg, ix%end + + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j_idx, k_idx, l_idx) + & + velR_vf(2)%sf(j_idx, k_idx + 1, l_idx)) ! v_r at interface k_idx+1/2 + + do i_loop = 1, 2 ! u_axial (1), u_radial (2) + dvel_avg_dx(i_loop) = 0.5_wp*(dvelL_dx_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & + dvelR_dx_vf(i_loop)%sf(j_idx, k_idx + 1, l_idx)) ! du_i/dx_axial + dvel_avg_dy(i_loop) = 0.5_wp*(dvelL_dy_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & + dvelR_dy_vf(i_loop)%sf(j_idx, k_idx + 1, l_idx)) ! du_i/dr_normal end do - end if - if (bulk_stress) then - stress_normal_bulk = divergence_cyl/Re_b + if (shear_stress) then + ! tau_rx_axial shear stress (r normal, x axial) + tau_Re(2, 1) = (dvel_avg_dy(1) + dvel_avg_dx(2))/Re_avg_rsy_vf(k_idx, j_idx, l_idx, 1) + ! tau_rr_shear (normal-normal stress): (4/3)mu dv_r/dr - (2/3)mu (dv_axial/dx_axial + v_r/r + (1/r)dv_theta/dtheta) + ! This specific form is (4/3)mu dv_r/dr - (2/3)mu (dv_axial/dx_axial + v_r/r) + tau_Re(2, 2) = (4._wp*dvel_avg_dy(2) - 2._wp*dvel_avg_dx(1) - & + 2._wp*avg_vel(2)/y_cb(k_idx))/(3._wp*Re_avg_rsy_vf(k_idx, j_idx, l_idx, 1)) + + do i_loop = 1, 2 ! Axial (1) and radial (2) momentum components + flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) = & + flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) - tau_Re(2, i_loop) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = & + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsy_vf(k_idx, j_idx, l_idx, i_loop)*tau_Re(2, i_loop) + end do + end if - flux_src_vf(momxb + norm_dir - 1)%sf(j, k, l) = flux_src_vf(momxb + norm_dir - 1)%sf(j, k, l) - stress_normal_bulk - flux_src_vf(E_idx)%sf(j, k, l) = flux_src_vf(E_idx)%sf(j, k, l) - vel_src_int(norm_dir)*stress_normal_bulk - end if + if (bulk_stress) then + ! tau_rr_bulk: mu_bulk (div V) = mu_bulk (dv_axial/dx_axial + dv_r/dr + v_r/r + (1/r)dv_theta/dtheta) + ! This specific form includes axial and radial terms + tau_Re(2, 2) = (dvel_avg_dx(1) + dvel_avg_dy(2) + avg_vel(2)/y_cb(k_idx))/Re_avg_rsy_vf(k_idx, j_idx, l_idx, 2) + flux_src_vf(momxb + 1)%sf(j_idx, k_idx, l_idx) = flux_src_vf(momxb + 1)%sf(j_idx, k_idx, l_idx) - tau_Re(2, 2) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsy_vf(k_idx, j_idx, l_idx, 2)*tau_Re(2, 2) + end if + end do + end do + end do + + if (p == 0) return ! No azimuthal extent + + ! Part 2: Azimuthal derivative contributions to radial flux (3D cylindrical) + !$acc parallel loop collapse(3) gang vector default(present) & + !$acc private(avg_vel, dvel_avg_dy, dvel_avg_dz, tau_Re, i_loop) + do l_idx = iz%beg, iz%end + do k_idx = iy%beg, iy%end ! Radial flux index + do j_idx = ix%beg, ix%end + + avg_vel(3) = 0.5_wp*(velL_vf(3)%sf(j_idx, k_idx, l_idx) + & + velR_vf(3)%sf(j_idx, k_idx + 1, l_idx)) ! v_theta at interface + do i_loop = 2, 3 ! u_radial (2), u_azimuthal (3) + dvel_avg_dz(i_loop) = 0.5_wp*(dvelL_dz_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & + dvelR_dz_vf(i_loop)%sf(j_idx, k_idx + 1, l_idx)) ! du_i/(r dtheta) + end do + dvel_avg_dy(3) = 0.5_wp*(dvelL_dy_vf(3)%sf(j_idx, k_idx, l_idx) + & + dvelR_dy_vf(3)%sf(j_idx, k_idx + 1, l_idx)) ! du_theta/dr_normal + + if (shear_stress) then + ! tau_rr component from azimuthal terms: -(2/3)mu ( (1/r)dv_theta/dtheta ) + tau_Re(2, 2) = -(2._wp/3._wp)*(dvel_avg_dz(3)/y_cb(k_idx))/Re_avg_rsy_vf(k_idx, j_idx, l_idx, 1) + ! tau_r_theta shear stress: mu ( (1/r)dv_r/dtheta - v_theta/r + dv_theta/dr ) + tau_Re(2, 3) = ((dvel_avg_dz(2) - avg_vel(3))/y_cb(k_idx) + dvel_avg_dy(3))/Re_avg_rsy_vf(k_idx, j_idx, l_idx, 1) + + do i_loop = 2, 3 ! Radial (2) and azimuthal (3) momentum components + flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) = & + flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) - tau_Re(2, i_loop) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = & + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsy_vf(k_idx, j_idx, l_idx, i_loop)*tau_Re(2, i_loop) + end do + end if + + if (bulk_stress) then + ! tau_rr_bulk from azimuthal terms: mu_bulk ( (1/r)dv_theta/dtheta ) + tau_Re(2, 2) = (dvel_avg_dz(3)/y_cb(k_idx))/Re_avg_rsy_vf(k_idx, j_idx, l_idx, 2) + flux_src_vf(momxb + 1)%sf(j_idx, k_idx, l_idx) = flux_src_vf(momxb + 1)%sf(j_idx, k_idx, l_idx) - tau_Re(2, 2) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsy_vf(k_idx, j_idx, l_idx, 2)*tau_Re(2, 2) + end if + end do end do end do - end do - !$acc end parallel loop + + ! norm_dir = 3: Azimuthal direction flux (e.g., theta-flux if grid z-axis is solver's theta-axis) + ! In this block, for derivatives: dz is azimuthal (normal), dx is axial, dy is radial component of gradient. + else + !$acc parallel loop collapse(3) gang vector default(present) & + !$acc private(avg_vel, dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re, i_loop) + do l_idx = iz%beg, iz%end ! Azimuthal flux index + do k_idx = iy%beg, iy%end + do j_idx = ix%beg, ix%end + + do i_loop = 2, 3 ! u_radial (2), u_azimuthal (3) + avg_vel(i_loop) = 0.5_wp*(velL_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & + velR_vf(i_loop)%sf(j_idx, k_idx, l_idx + 1)) + end do + + do i_loop = 1, 3, 2 ! u_axial (1), u_azimuthal (3) + dvel_avg_dx(i_loop) = 0.5_wp*(dvelL_dx_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & + dvelR_dx_vf(i_loop)%sf(j_idx, k_idx, l_idx + 1)) ! du_i/dx_axial + end do + do i_loop = 2, 3 ! u_radial (2), u_azimuthal (3) + dvel_avg_dy(i_loop) = 0.5_wp*(dvelL_dy_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & + dvelR_dy_vf(i_loop)%sf(j_idx, k_idx, l_idx + 1)) ! du_i/dr_radial + end do + do i_loop = 1, 3 ! All components + dvel_avg_dz(i_loop) = 0.5_wp*(dvelL_dz_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & + dvelR_dz_vf(i_loop)%sf(j_idx, k_idx, l_idx + 1)) ! du_i/(r dtheta)_normal + end do + + if (shear_stress) then + ! tau_theta_axial shear stress + tau_Re(3, 1) = (dvel_avg_dz(1)/y_cc(k_idx) + dvel_avg_dx(3))/Re_avg_rsz_vf(l_idx, k_idx, j_idx, 1) + ! tau_theta_radial shear stress + tau_Re(3, 2) = ((dvel_avg_dz(2) - avg_vel(3))/y_cc(k_idx) + dvel_avg_dy(3))/Re_avg_rsz_vf(l_idx, k_idx, j_idx, 1) + ! tau_theta_theta_shear (normal-normal stress) + tau_Re(3, 3) = (4._wp*dvel_avg_dz(3)/y_cc(k_idx) - 2._wp*dvel_avg_dx(1) - & + 2._wp*dvel_avg_dy(2) + 4._wp*avg_vel(2)/y_cc(k_idx))/ & + (3._wp*Re_avg_rsz_vf(l_idx, k_idx, j_idx, 1)) + + ! Final scaling by 1/r for azimuthal flux (original code implies this for tau_Re in this context) + tau_Re(3, 1) = tau_Re(3, 1)/y_cc(k_idx) + tau_Re(3, 2) = tau_Re(3, 2)/y_cc(k_idx) + tau_Re(3, 3) = tau_Re(3, 3)/y_cc(k_idx) + + do i_loop = 1, 3 ! Axial, radial, azimuthal momentum + flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) = & + flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) - tau_Re(3, i_loop) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = & + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsz_vf(l_idx, k_idx, j_idx, i_loop)*tau_Re(3, i_loop) + end do + end if + + if (bulk_stress) then + ! tau_theta_theta_bulk + tau_Re(3, 3) = (dvel_avg_dx(1) + dvel_avg_dy(2) + dvel_avg_dz(3)/y_cc(k_idx) + avg_vel(2)/y_cc(k_idx))/ & + Re_avg_rsz_vf(l_idx, k_idx, j_idx, 2) + + ! Final scaling by 1/r + tau_Re(3, 3) = tau_Re(3, 3)/y_cc(k_idx) + + flux_src_vf(momxe)%sf(j_idx, k_idx, l_idx) = flux_src_vf(momxe)%sf(j_idx, k_idx, l_idx) - tau_Re(3, 3) + flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) = flux_src_vf(E_idx)%sf(j_idx, k_idx, l_idx) - & + vel_src_rsz_vf(l_idx, k_idx, j_idx, 3)*tau_Re(3, 3) + end if + end do + end do + end do + end if end subroutine s_compute_cylindrical_viscous_source_flux @@ -4114,9 +4288,6 @@ contains !! @param[in] dvelR_dz_vf Right boundary d(vel)/dz (num_dims scalar_field). !! @param[inout] flux_src_vf Intercell source flux array to update (sys_size scalar_field). !! @param[in] norm_dir Interface normal direction (1=x, 2=y, 3=z). - !! @param[in] ix X-direction loop bounds (int_bounds_info). - !! @param[in] iy Y-direction loop bounds (int_bounds_info). - !! @param[in] iz Z-direction loop bounds (int_bounds_info). subroutine s_compute_cartesian_viscous_source_flux(velL_vf, & dvelL_dx_vf, & dvelL_dy_vf, & From 0867211e6b62447a468a06af0d3d5f94ea2486a6 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 27 May 2025 22:16:03 -0400 Subject: [PATCH 6/8] add loop seq --- src/simulation/m_riemann_solvers.fpp | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 8c1622ecbc..ac4c6079b0 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -4035,6 +4035,7 @@ contains avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j_idx, k_idx, l_idx) + velR_vf(2)%sf(j_idx + 1, k_idx, l_idx)) ! v_r + !$acc loop seq do i_loop = 1, 2 ! u_axial (1), u_radial (2) dvel_avg_dy(i_loop) = 0.5_wp*(dvelL_dy_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & dvelR_dy_vf(i_loop)%sf(j_idx + 1, k_idx, l_idx)) ! du_i/dr @@ -4049,6 +4050,7 @@ contains ! tau_zr shear stress: mu (dv_axial/dr + dv_r/dx_axial) tau_Re(1, 2) = (dvel_avg_dy(1) + dvel_avg_dx(2))/Re_avg_rsx_vf(j_idx, k_idx, l_idx, 1) + !$acc loop seq do i_loop = 1, 2 ! Axial (1) and radial (2) momentum components flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) = & flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) - tau_Re(1, i_loop) @@ -4078,6 +4080,7 @@ contains do k_idx = iy%beg, iy%end do j_idx = ix%beg, ix%end + !$acc loop seq do i_loop = 1, 3, 2 ! u_axial (1), u_azimuthal (3) dvel_avg_dz(i_loop) = 0.5_wp*(dvelL_dz_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & dvelR_dz_vf(i_loop)%sf(j_idx + 1, k_idx, l_idx)) ! du_i/(r dtheta) @@ -4091,6 +4094,7 @@ contains ! tau_z_theta shear stress: mu ( (1/r)dv_axial/dtheta + dv_theta/dx_axial ) tau_Re(1, 3) = (dvel_avg_dz(1)/y_cc(k_idx) + dvel_avg_dx(3))/Re_avg_rsx_vf(j_idx, k_idx, l_idx, 1) + !$acc loop seq do i_loop = 1, 3, 2 ! Axial (1) and azimuthal (3) momentum components flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) = & flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) - tau_Re(1, i_loop) @@ -4124,6 +4128,7 @@ contains avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j_idx, k_idx, l_idx) + & velR_vf(2)%sf(j_idx, k_idx + 1, l_idx)) ! v_r at interface k_idx+1/2 + !$acc loop seq do i_loop = 1, 2 ! u_axial (1), u_radial (2) dvel_avg_dx(i_loop) = 0.5_wp*(dvelL_dx_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & dvelR_dx_vf(i_loop)%sf(j_idx, k_idx + 1, l_idx)) ! du_i/dx_axial @@ -4172,6 +4177,7 @@ contains avg_vel(3) = 0.5_wp*(velL_vf(3)%sf(j_idx, k_idx, l_idx) + & velR_vf(3)%sf(j_idx, k_idx + 1, l_idx)) ! v_theta at interface + !$acc loop seq do i_loop = 2, 3 ! u_radial (2), u_azimuthal (3) dvel_avg_dz(i_loop) = 0.5_wp*(dvelL_dz_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & dvelR_dz_vf(i_loop)%sf(j_idx, k_idx + 1, l_idx)) ! du_i/(r dtheta) @@ -4185,6 +4191,7 @@ contains ! tau_r_theta shear stress: mu ( (1/r)dv_r/dtheta - v_theta/r + dv_theta/dr ) tau_Re(2, 3) = ((dvel_avg_dz(2) - avg_vel(3))/y_cb(k_idx) + dvel_avg_dy(3))/Re_avg_rsy_vf(k_idx, j_idx, l_idx, 1) + !$acc loop seq do i_loop = 2, 3 ! Radial (2) and azimuthal (3) momentum components flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) = & flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) - tau_Re(2, i_loop) @@ -4214,19 +4221,25 @@ contains do k_idx = iy%beg, iy%end do j_idx = ix%beg, ix%end + !$acc loop seq do i_loop = 2, 3 ! u_radial (2), u_azimuthal (3) avg_vel(i_loop) = 0.5_wp*(velL_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & velR_vf(i_loop)%sf(j_idx, k_idx, l_idx + 1)) end do + !$acc loop seq do i_loop = 1, 3, 2 ! u_axial (1), u_azimuthal (3) dvel_avg_dx(i_loop) = 0.5_wp*(dvelL_dx_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & dvelR_dx_vf(i_loop)%sf(j_idx, k_idx, l_idx + 1)) ! du_i/dx_axial end do + + !$acc loop seq do i_loop = 2, 3 ! u_radial (2), u_azimuthal (3) dvel_avg_dy(i_loop) = 0.5_wp*(dvelL_dy_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & dvelR_dy_vf(i_loop)%sf(j_idx, k_idx, l_idx + 1)) ! du_i/dr_radial end do + + !$acc loop seq do i_loop = 1, 3 ! All components dvel_avg_dz(i_loop) = 0.5_wp*(dvelL_dz_vf(i_loop)%sf(j_idx, k_idx, l_idx) + & dvelR_dz_vf(i_loop)%sf(j_idx, k_idx, l_idx + 1)) ! du_i/(r dtheta)_normal @@ -4247,6 +4260,7 @@ contains tau_Re(3, 2) = tau_Re(3, 2)/y_cc(k_idx) tau_Re(3, 3) = tau_Re(3, 3)/y_cc(k_idx) + !$acc loop seq do i_loop = 1, 3 ! Axial, radial, azimuthal momentum flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) = & flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) - tau_Re(3, i_loop) @@ -4341,6 +4355,8 @@ contains idx_right_phys(norm_dir) = idx_right_phys(norm_dir) + 1 vel_grad_avg = 0.0_wp + + !$acc loop seq do vel_comp_idx = 1, num_dims vel_grad_avg(vel_comp_idx, 1) = 0.5_wp*(dvelL_dx_vf(vel_comp_idx)%sf(j_loop, k_loop, l_loop) + & dvelR_dx_vf(vel_comp_idx)%sf(idx_right_phys(1), idx_right_phys(2), idx_right_phys(3))) @@ -4355,6 +4371,7 @@ contains end do divergence_v = 0.0_wp + !$acc loop seq do i_dim = 1, num_dims divergence_v = divergence_v + vel_grad_avg(i_dim, i_dim) end do @@ -4363,18 +4380,21 @@ contains if (norm_dir == 1) then Re_shear = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 1) Re_bulk = Re_avg_rsx_vf(j_loop, k_loop, l_loop, 2) + !$acc loop seq do i_dim = 1, num_dims vel_src_at_interface(i_dim) = vel_src_rsx_vf(j_loop, k_loop, l_loop, i_dim) end do else if (norm_dir == 2) then Re_shear = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 1) Re_bulk = Re_avg_rsy_vf(k_loop, j_loop, l_loop, 2) + !$acc loop seq do i_dim = 1, num_dims vel_src_at_interface(i_dim) = vel_src_rsy_vf(k_loop, j_loop, l_loop, i_dim) end do else Re_shear = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 1) Re_bulk = Re_avg_rsz_vf(l_loop, k_loop, j_loop, 2) + !$acc loop seq do i_dim = 1, num_dims vel_src_at_interface(i_dim) = vel_src_rsz_vf(l_loop, k_loop, j_loop, i_dim) end do @@ -4384,6 +4404,7 @@ contains current_tau_shear = 0.0_wp call s_calculate_shear_stress_tensor(vel_grad_avg, Re_shear, divergence_v, current_tau_shear) + !$acc loop seq do i_dim = 1, num_dims flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) = & flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) - current_tau_shear(norm_dir, i_dim) @@ -4398,6 +4419,7 @@ contains current_tau_bulk = 0.0_wp call s_calculate_bulk_stress_tensor(Re_bulk, divergence_v, current_tau_bulk) + !$acc loop seq do i_dim = 1, num_dims flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) = & flux_src_vf(momxb + i_dim - 1)%sf(j_loop, k_loop, l_loop) - current_tau_bulk(norm_dir, i_dim) From 25e16156071e120487b634ec9193d2a740e159b2 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 28 May 2025 08:09:50 -0400 Subject: [PATCH 7/8] Update m_riemann_solvers.fpp --- src/simulation/m_riemann_solvers.fpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index ac4c6079b0..36539c4b4d 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -4144,6 +4144,7 @@ contains tau_Re(2, 2) = (4._wp*dvel_avg_dy(2) - 2._wp*dvel_avg_dx(1) - & 2._wp*avg_vel(2)/y_cb(k_idx))/(3._wp*Re_avg_rsy_vf(k_idx, j_idx, l_idx, 1)) + !$acc loop seq do i_loop = 1, 2 ! Axial (1) and radial (2) momentum components flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) = & flux_src_vf(contxe + i_loop)%sf(j_idx, k_idx, l_idx) - tau_Re(2, i_loop) From 0db9b67198580672eeb5483d5c3c05e70df4c03a Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Thu, 29 May 2025 21:55:43 -0400 Subject: [PATCH 8/8] Update case.py --- toolchain/mfc/test/case.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/toolchain/mfc/test/case.py b/toolchain/mfc/test/case.py index 3720aaee49..97d4fd121a 100644 --- a/toolchain/mfc/test/case.py +++ b/toolchain/mfc/test/case.py @@ -255,6 +255,8 @@ def compute_tolerance(self) -> float: if "Example" in self.trace.split(" -> "): tolerance = 1e-3 + elif "Cylindrical" in self.trace.split(" -> "): + tolerance = 1e-9 elif self.params.get("hypoelasticity", 'F') == 'T': tolerance = 1e-7 elif self.params.get("mixlayer_perturb", 'F') == 'T':