Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly
Appearance settings

Make CBC not so terrible #870

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Jun 9, 2025
Merged

Make CBC not so terrible #870

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d11643f
refac cbc
sbryngelson Jun 8, 2025
8b680c0
refac cbc
sbryngelson Jun 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
300 changes: 112 additions & 188 deletions 300 src/simulation/m_compute_cbc.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,11 +1,9 @@
!>
!! @file m_compute_cbc.f90
!! @brief Contains module m_compute_cbc
!! @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, &
Expand All @@ -18,11 +16,72 @@
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)))

Check warning on line 26 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L26 

Added line #L26 was not covered by tests
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)

Check warning on line 39 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L39 

Added line #L39 was not covered by tests
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))

Check warning on line 52 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L52 

Added line #L52 was not covered by tests
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)

Check warning on line 65 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L65 

Added line #L65 was not covered by tests
end do
end subroutine s_fill_advection_L

!> Fill chemistry L variables
pure subroutine s_fill_chemistry_L(L, lambda_factor, lambda2, dYs_ds)

Check warning on line 70 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L70 

Added line #L70 was not covered by tests
!$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)

Check warning on line 80 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L79-L80 

Added lines #L79 - L80 were not covered by tests
end do
end subroutine s_fill_chemistry_L

Check warning on line 82 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L82 

Added line #L82 was not covered by tests

!> 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
Expand All @@ -31,26 +90,16 @@
#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

Check warning on line 98 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L97-L98 

Added lines #L97 - L98 were not covered by tests
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
Expand All @@ -65,42 +114,22 @@
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)))

Check warning on line 120 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L120 

Added line #L120 was not covered by tests

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)

Check warning on line 126 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L126 

Added line #L126 was not covered by tests

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)))

Check warning on line 129 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L129 

Added line #L129 was not covered by tests
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
Expand All @@ -109,30 +138,15 @@
#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

Check warning on line 145 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L144-L145 

Added lines #L144 - L145 were not covered by tests
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
Expand All @@ -148,40 +162,15 @@
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)

Check warning on line 165 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L165 

Added line #L165 was not covered by tests
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)

Check warning on line 169 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L169 

Added line #L169 was not covered by tests
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
Expand All @@ -196,30 +185,14 @@
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)

Check warning on line 188 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L188 

Added line #L188 was not covered by tests
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
Expand All @@ -234,57 +207,26 @@
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)

Check warning on line 210 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L210 

Added line #L210 was not covered by tests
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
#else
!$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

Check warning on line 225 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L225 

Added line #L225 was not covered by tests
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
Expand All @@ -299,30 +241,12 @@
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)

Check warning on line 245 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L245 

Added line #L245 was not covered by tests
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)

Check warning on line 249 in src/simulation/m_compute_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_compute_cbc.fpp#L249 

Added line #L249 was not covered by tests
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
Loading
Morty Proxy This is a proxified and sanitized view of the page, visit original site.