diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index df388c479a..36539c4b4d 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -3940,22 +3940,24 @@ 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 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, & @@ -3968,521 +3970,339 @@ 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 + implicit none - type(scalar_field), & - dimension(sys_size), & - intent(inout) :: flux_src_vf + 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 - integer, intent(in) :: norm_dir + 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 + real(wp), dimension(num_dims, num_dims) :: tau_Re ! Viscous stress tensor components - ! Generic loop iterators - integer :: i, j, k, l + integer :: i_loop ! Loop iterator for components + integer :: j_idx, k_idx, l_idx ! Spatial loop iterators - ! Viscous Stresses in z-direction + ! 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 - 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) + ! 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 - end do + 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 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)) + end do - tau_Re(1, 1) = dvel_avg_dx(1)/ & - Re_avg_rsx_vf(j, k, l, 2) + if (n == 0) return ! No radial extent - flux_src_vf(momxb)%sf(j, k, l) = & - flux_src_vf(momxb)%sf(j, k, l) - & - tau_Re(1, 1) + ! 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 - 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) + 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 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) + 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) + + !$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) + 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 - end do + 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 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) + end do - flux_src_vf(momxb)%sf(j, k, l) = & - flux_src_vf(momxb)%sf(j, k, l) - & - tau_Re(1, 1) + if (p == 0) return ! No azimuthal extent - 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) + ! 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 + !$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) 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) - + 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) + + !$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) + 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 - 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 + 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 if - ! END: Viscous Stresses in z-direction + end do - ! Viscous Stresses in r-direction + ! 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 + + !$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 + 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 - 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)) + 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)) !$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) - + 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 - end do + 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 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)) + end do - 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) + if (p == 0) return ! No azimuthal extent - flux_src_vf(momxb + 1)%sf(j, k, l) = & - flux_src_vf(momxb + 1)%sf(j, k, l) - & - tau_Re(2, 2) + ! 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 - 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) + 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) 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) - + 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) + + !$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) + 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 - 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)) - - tau_Re(2, 2) = dvel_avg_dz(3)/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 + 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 if - ! END: Viscous Stresses in r-direction + end do - ! Viscous Stresses in theta-direction + ! 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 - - 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) - - 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) - - 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) - - !$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 - + !$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 + + !$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 - 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)) + !$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 - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k, l + 1)) + !$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 - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k, l + 1)) + !$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 + end do - 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) + 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) + + !$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) + 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 - flux_src_vf(momxe)%sf(j, k, l) = & - flux_src_vf(momxe)%sf(j, k, l) - & - tau_Re(3, 3) + 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) - 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) + ! Final scaling by 1/r + tau_Re(3, 3) = tau_Re(3, 3)/y_cc(k_idx) - end do + 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 if - + end do end if - ! END: Viscous Stresses in theta-direction 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). subroutine s_compute_cartesian_viscous_source_flux(velL_vf, & dvelL_dx_vf, & dvelL_dy_vf, & @@ -4495,469 +4315,189 @@ 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 + ! 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 + + !$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))) + 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 (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 + 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 - 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 + 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) + !$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 - 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) - + 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 - 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 - + 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 - 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) - - 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 - - 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) - - 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 - - flux_src_vf(contxe + i)%sf(j, k, l) = & - flux_src_vf(contxe + i)%sf(j, k, l) - & - tau_Re(2, i) + end if - 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) + 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) - end do + !$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) + 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 - - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k + 1, l)) - - tau_Re(2, 2) = dvel_avg_dz(3)/ & - Re_avg_rsy_vf(k, j, l, 2) + end if - flux_src_vf(momxb + 1)%sf(j, k, l) = & - flux_src_vf(momxb + 1)%sf(j, k, l) - & - tau_Re(2, 2) + if (bulk_stress) then + current_tau_bulk = 0.0_wp + call s_calculate_bulk_stress_tensor(Re_bulk, divergence_v, current_tau_bulk) - 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) + !$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) + 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) + !$acc routine seq + + 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 + + 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_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k, l + 1)) + end subroutine s_calculate_shear_stress_tensor - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k, l + 1)) + !> @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) + !$acc routine seq - tau_Re(3, 3) = (dvel_avg_dx(1) & - + dvel_avg_dy(2) & - + dvel_avg_dz(3))/ & - Re_avg_rsz_vf(l, k, j, 2) + 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 - end if - ! END: Viscous Stresses in z-direction + do i_dim = 1, num_dims + tau_bulk_out(i_dim, i_dim) = divergence_v/Re_bulk + end do - end subroutine s_compute_cartesian_viscous_source_flux + end subroutine s_calculate_bulk_stress_tensor !> Deallocation and/or disassociation procedures that are !! needed to finalize the selected Riemann problem solver 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':