From d11643fbd880ce7f9a66e77c01ef19253f004bc4 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sun, 8 Jun 2025 18:04:13 -0400 Subject: [PATCH 1/2] refac cbc --- src/simulation/m_compute_cbc.fpp | 302 ++++++++++++------------------- 1 file changed, 113 insertions(+), 189 deletions(-) diff --git a/src/simulation/m_compute_cbc.fpp b/src/simulation/m_compute_cbc.fpp index 7bab7d3a67..28180aaf71 100644 --- a/src/simulation/m_compute_cbc.fpp +++ b/src/simulation/m_compute_cbc.fpp @@ -1,11 +1,9 @@ !> -!! @file m_compute_cbc.f90 -!! @brief Contains module m_compute_cbc +!! @file m_compute_cbc.fpp +!! @brief CBC computation module module m_compute_cbc - - use m_global_parameters !< Definitions of the global parameters - + use m_global_parameters implicit none private; public :: s_compute_slip_wall_L, & @@ -18,11 +16,72 @@ module m_compute_cbc s_compute_supersonic_outflow_L contains + !> Base L1 calculation + pure function f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) result(L1) + !$acc routine seq + real(wp), dimension(3), intent(in) :: lambda + real(wp), intent(in) :: rho, c, dpres_ds + real(wp), dimension(num_dims), intent(in) :: dvel_ds + real(wp) :: L1 + L1 = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1))) + end function f_base_L1 + + !> Fill density L variables + pure subroutine s_fill_density_L(L, lambda_factor, lambda2, c, mf, dalpha_rho_ds, dpres_ds) + !$acc routine seq + real(wp), dimension(sys_size), intent(inout) :: L + real(wp), intent(in) :: lambda_factor, lambda2, c + real(wp), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds + real(wp), intent(in) :: dpres_ds + integer :: i + + do i = 2, momxb + L(i) = lambda_factor*lambda2*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds) + end do + end subroutine s_fill_density_L - !> The L variables for the slip wall CBC, see pg. 451 of - !! Thompson (1990). At the slip wall (frictionless wall), - !! the normal component of velocity is zero at all times, - !! while the transverse velocities may be nonzero. + !> Fill velocity L variables + pure subroutine s_fill_velocity_L(L, lambda_factor, lambda2, dvel_ds) + !$acc routine seq + real(wp), dimension(sys_size), intent(inout) :: L + real(wp), intent(in) :: lambda_factor, lambda2 + real(wp), dimension(num_dims), intent(in) :: dvel_ds + integer :: i + + do i = momxb + 1, momxe + L(i) = lambda_factor*lambda2*dvel_ds(dir_idx(i - contxe)) + end do + end subroutine s_fill_velocity_L + + !> Fill advection L variables + pure subroutine s_fill_advection_L(L, lambda_factor, lambda2, dadv_ds) + !$acc routine seq + real(wp), dimension(sys_size), intent(inout) :: L + real(wp), intent(in) :: lambda_factor, lambda2 + real(wp), dimension(num_fluids), intent(in) :: dadv_ds + integer :: i + + do i = E_idx, advxe - 1 + L(i) = lambda_factor*lambda2*dadv_ds(i - momxe) + end do + end subroutine s_fill_advection_L + + !> Fill chemistry L variables + pure subroutine s_fill_chemistry_L(L, lambda_factor, lambda2, dYs_ds) + !$acc routine seq + real(wp), dimension(sys_size), intent(inout) :: L + real(wp), intent(in) :: lambda_factor, lambda2 + real(wp), dimension(num_species), intent(in) :: dYs_ds + integer :: i + + if (.not. chemistry) return + + do i = chemxb, chemxe + L(i) = lambda_factor*lambda2*dYs_ds(i - chemxb + 1) + end do + end subroutine s_fill_chemistry_L + + !> Slip wall CBC (Thompson 1990, pg. 451) pure subroutine s_compute_slip_wall_L(lambda, L, rho, c, dpres_ds, dvel_ds) #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_slip_wall_L @@ -31,26 +90,16 @@ contains #endif real(wp), dimension(3), intent(in) :: lambda real(wp), dimension(sys_size), intent(inout) :: L - real(wp), intent(in) :: rho, c - real(wp), intent(in) :: dpres_ds + real(wp), intent(in) :: rho, c, dpres_ds real(wp), dimension(num_dims), intent(in) :: dvel_ds - integer :: i - L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1))) - - do i = 2, advxe - L(i) = 0._wp - end do - + L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) + L(2:advxe - 1) = 0._wp L(advxe) = L(1) - end subroutine s_compute_slip_wall_L - !> The L variables for the nonreflecting subsonic buffer CBC - !! see pg. 13 of Thompson (1987). The nonreflecting subsonic - !! buffer reduces the amplitude of any reflections caused by - !! outgoing waves. + !> Nonreflecting subsonic buffer CBC (Thompson 1987, pg. 13) pure subroutine s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds) #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_buffer_L @@ -65,42 +114,22 @@ contains real(wp), dimension(num_dims), intent(in) :: dvel_ds real(wp), dimension(num_fluids), intent(in) :: dadv_ds real(wp), dimension(num_species), intent(in) :: dYs_ds + real(wp) :: lambda_factor - integer :: i !< Generic loop iterator - - L(1) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(1)))*lambda(1) & - *(dpres_ds - rho*c*dvel_ds(dir_idx(1))) + lambda_factor = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(1))) + L(1) = lambda_factor*lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1))) - do i = 2, momxb - L(i) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(2)))*lambda(2) & - *(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds) - end do - - do i = momxb + 1, momxe - L(i) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(2)))*lambda(2) & - *(dvel_ds(dir_idx(i - contxe))) - end do - - do i = E_idx, advxe - 1 - L(i) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(2)))*lambda(2) & - *(dadv_ds(i - momxe)) - end do - - L(advxe) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(3)))*lambda(3) & - *(dpres_ds + rho*c*dvel_ds(dir_idx(1))) - - if (chemistry) then - do i = chemxb, chemxe - L(i) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(2)))*lambda(2) & - *(dYs_ds(i - chemxb + 1)) - end do - end if + lambda_factor = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(2))) + call s_fill_density_L(L, lambda_factor, lambda(2), c, mf, dalpha_rho_ds, dpres_ds) + call s_fill_velocity_L(L, lambda_factor, lambda(2), dvel_ds) + call s_fill_advection_L(L, lambda_factor, lambda(2), dadv_ds) + call s_fill_chemistry_L(L, lambda_factor, lambda(2), dYs_ds) + lambda_factor = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(3))) + L(advxe) = lambda_factor*lambda(3)*(dpres_ds + rho*c*dvel_ds(dir_idx(1))) end subroutine s_compute_nonreflecting_subsonic_buffer_L - !> The L variables for the nonreflecting subsonic inflow CBC - !! see pg. 455, Thompson (1990). This nonreflecting subsonic - !! CBC assumes an incoming flow and reduces the amplitude of - !! any reflections caused by outgoing waves. + + !> Nonreflecting subsonic inflow CBC (Thompson 1990, pg. 455) pure subroutine s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, dpres_ds, dvel_ds) #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_inflow_L @@ -109,30 +138,15 @@ contains #endif real(wp), dimension(3), intent(in) :: lambda real(wp), dimension(sys_size), intent(inout) :: L - real(wp), intent(in) :: rho, c - real(wp), intent(in) :: dpres_ds + real(wp), intent(in) :: rho, c, dpres_ds real(wp), dimension(num_dims), intent(in) :: dvel_ds - integer :: i - - L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1))) - - do i = 2, advxe - L(i) = 0._wp - end do - - if (chemistry) then - do i = chemxb, chemxe - L(i) = 0._wp - end do - end if - + L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) + L(2:advxe) = 0._wp + if (chemistry) L(chemxb:chemxe) = 0._wp end subroutine s_compute_nonreflecting_subsonic_inflow_L - !> The L variables for the nonreflecting subsonic outflow - !! CBC see pg. 454 of Thompson (1990). This nonreflecting - !! subsonic CBC presumes an outgoing flow and reduces the - !! amplitude of any reflections caused by outgoing waves. + !> Nonreflecting subsonic outflow CBC (Thompson 1990, pg. 454) pure subroutine s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds) #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_outflow_L @@ -148,40 +162,15 @@ contains real(wp), dimension(num_fluids), intent(in) :: dadv_ds real(wp), dimension(num_species), intent(in) :: dYs_ds - integer :: i !> Generic loop iterator - - L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1))) - - do i = 2, momxb - L(i) = lambda(2)*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds) - end do - - do i = momxb + 1, momxe - L(i) = lambda(2)*(dvel_ds(dir_idx(i - contxe))) - end do - - do i = E_idx, advxe - 1 - L(i) = lambda(2)*(dadv_ds(i - momxe)) - end do - - ! bubble index + L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) + call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds) + call s_fill_velocity_L(L, 1._wp, lambda(2), dvel_ds) + call s_fill_advection_L(L, 1._wp, lambda(2), dadv_ds) + call s_fill_chemistry_L(L, 1._wp, lambda(2), dYs_ds) L(advxe) = 0._wp - - if (chemistry) then - do i = chemxb, chemxe - L(i) = lambda(2)*dYs_ds(i - chemxb + 1) - end do - end if - end subroutine s_compute_nonreflecting_subsonic_outflow_L - !> The L variables for the force-free subsonic outflow CBC, - !! see pg. 454 of Thompson (1990). The force-free subsonic - !! outflow sets to zero the sum of all of the forces which - !! are acting on a fluid element for the normal coordinate - !! direction to the boundary. As a result, a fluid element - !! at the boundary is simply advected outward at the fluid - !! velocity. + !> Force-free subsonic outflow CBC (Thompson 1990, pg. 454) pure subroutine s_compute_force_free_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_force_free_subsonic_outflow_L @@ -196,30 +185,14 @@ contains real(wp), dimension(num_dims), intent(in) :: dvel_ds real(wp), dimension(num_fluids), intent(in) :: dadv_ds - integer :: i !> Generic loop iterator - - L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1))) - - do i = 2, momxb - L(i) = lambda(2)*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds) - end do - - do i = momxb + 1, momxe - L(i) = lambda(2)*(dvel_ds(dir_idx(i - contxe))) - end do - - do i = E_idx, advxe - 1 - L(i) = lambda(2)*(dadv_ds(i - momxe)) - end do - + L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) + call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds) + call s_fill_velocity_L(L, 1._wp, lambda(2), dvel_ds) + call s_fill_advection_L(L, 1._wp, lambda(2), dadv_ds) L(advxe) = L(1) + 2._wp*rho*c*lambda(2)*dvel_ds(dir_idx(1)) - end subroutine s_compute_force_free_subsonic_outflow_L - !> L variables for the constant pressure subsonic outflow - !! CBC see pg. 455 Thompson (1990). The constant pressure - !! subsonic outflow maintains a fixed pressure at the CBC - !! boundary in absence of any transverse effects. + !> Constant pressure subsonic outflow CBC (Thompson 1990, pg. 455) pure subroutine s_compute_constant_pressure_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_constant_pressure_subsonic_outflow_L @@ -234,31 +207,14 @@ contains real(wp), dimension(num_dims), intent(in) :: dvel_ds real(wp), dimension(num_fluids), intent(in) :: dadv_ds - integer :: i !> Generic loop iterator - - L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1))) - - do i = 2, momxb - L(i) = lambda(2)*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds) - end do - - do i = momxb + 1, momxe - L(i) = lambda(2)*(dvel_ds(dir_idx(i - contxe))) - end do - - do i = E_idx, advxe - 1 - L(i) = lambda(2)*(dadv_ds(i - momxe)) - end do - + L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) + call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds) + call s_fill_velocity_L(L, 1._wp, lambda(2), dvel_ds) + call s_fill_advection_L(L, 1._wp, lambda(2), dadv_ds) L(advxe) = -L(1) - end subroutine s_compute_constant_pressure_subsonic_outflow_L - !> L variables for the supersonic inflow CBC, see pg. 453 - !! Thompson (1990). The supersonic inflow CBC is a steady - !! state, or nearly a steady state, CBC in which only the - !! transverse terms may generate a time dependence at the - !! inflow boundary. + !> Supersonic inflow CBC (Thompson 1990, pg. 453) pure subroutine s_compute_supersonic_inflow_L(L) #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_supersonic_inflow_L @@ -266,25 +222,11 @@ contains !$acc routine seq #endif real(wp), dimension(sys_size), intent(inout) :: L - - integer :: i - - do i = 1, advxe - L(i) = 0._wp - end do - - if (chemistry) then - do i = chemxb, chemxe - L(i) = 0._wp - end do - end if - + L(1:advxe) = 0._wp + if (chemistry) L(chemxb:chemxe) = 0._wp end subroutine s_compute_supersonic_inflow_L - !> L variables for the supersonic outflow CBC, see pg. 453 - !! of Thompson (1990). For the supersonic outflow CBC, the - !! flow evolution at the boundary is determined completely - !! by the interior data. + !> Supersonic outflow CBC (Thompson 1990, pg. 453) pure subroutine s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds) #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_supersonic_outflow_L @@ -299,30 +241,12 @@ contains real(wp), dimension(num_dims), intent(in) :: dvel_ds real(wp), dimension(num_fluids), intent(in) :: dadv_ds real(wp), dimension(num_species), intent(in) :: dYs_ds - integer :: i !< Generic loop iterator - - L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1))) - - do i = 2, momxb - L(i) = lambda(2)*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds) - end do - - do i = momxb + 1, momxe - L(i) = lambda(2)*(dvel_ds(dir_idx(i - contxe))) - end do - - do i = E_idx, advxe - 1 - L(i) = lambda(2)*(dadv_ds(i - momxe)) - end do + L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) + call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds) + call s_fill_velocity_L(L, 1._wp, lambda(2), dvel_ds) + call s_fill_advection_L(L, 1._wp, lambda(2), dadv_ds) + call s_fill_chemistry_L(L, 1._wp, lambda(2), dYs_ds) L(advxe) = lambda(3)*(dpres_ds + rho*c*dvel_ds(dir_idx(1))) - - if (chemistry) then - do i = chemxb, chemxe - L(i) = lambda(2)*dYs_ds(i - chemxb + 1) - end do - end if - end subroutine s_compute_supersonic_outflow_L - end module m_compute_cbc From 8b680c00d471ce5f4cb55fdeb8c31c417bdfbbb3 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sun, 8 Jun 2025 18:04:50 -0400 Subject: [PATCH 2/2] refac cbc --- src/simulation/m_compute_cbc.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_compute_cbc.fpp b/src/simulation/m_compute_cbc.fpp index 28180aaf71..47b3f9ce81 100644 --- a/src/simulation/m_compute_cbc.fpp +++ b/src/simulation/m_compute_cbc.fpp @@ -1,5 +1,5 @@ !> -!! @file m_compute_cbc.fpp +!! @file m_compute_cbc.f90 !! @brief CBC computation module module m_compute_cbc