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

Commit db44da1

Browse filesBrowse files
sbryngelsonsbryngelson
authored
Make CBC not so terrible (#870)
1 parent 9879eb7 commit db44da1
Copy full SHA for db44da1

File tree

Expand file treeCollapse file tree

1 file changed

+112
-188
lines changed
Filter options
Expand file treeCollapse file tree

1 file changed

+112
-188
lines changed

‎src/simulation/m_compute_cbc.fpp

Copy file name to clipboard
+112-188Lines changed: 112 additions & 188 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,9 @@
11
!>
22
!! @file m_compute_cbc.f90
3-
!! @brief Contains module m_compute_cbc
3+
!! @brief CBC computation module
44

55
module m_compute_cbc
6-
7-
use m_global_parameters !< Definitions of the global parameters
8-
6+
use m_global_parameters
97
implicit none
108

119
private; public :: s_compute_slip_wall_L, &
@@ -18,11 +16,72 @@ module m_compute_cbc
1816
s_compute_supersonic_outflow_L
1917

2018
contains
19+
!> Base L1 calculation
20+
pure function f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) result(L1)
21+
!$acc routine seq
22+
real(wp), dimension(3), intent(in) :: lambda
23+
real(wp), intent(in) :: rho, c, dpres_ds
24+
real(wp), dimension(num_dims), intent(in) :: dvel_ds
25+
real(wp) :: L1
26+
L1 = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
27+
end function f_base_L1
28+
29+
!> Fill density L variables
30+
pure subroutine s_fill_density_L(L, lambda_factor, lambda2, c, mf, dalpha_rho_ds, dpres_ds)
31+
!$acc routine seq
32+
real(wp), dimension(sys_size), intent(inout) :: L
33+
real(wp), intent(in) :: lambda_factor, lambda2, c
34+
real(wp), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds
35+
real(wp), intent(in) :: dpres_ds
36+
integer :: i
37+
38+
do i = 2, momxb
39+
L(i) = lambda_factor*lambda2*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds)
40+
end do
41+
end subroutine s_fill_density_L
2142

22-
!> The L variables for the slip wall CBC, see pg. 451 of
23-
!! Thompson (1990). At the slip wall (frictionless wall),
24-
!! the normal component of velocity is zero at all times,
25-
!! while the transverse velocities may be nonzero.
43+
!> Fill velocity L variables
44+
pure subroutine s_fill_velocity_L(L, lambda_factor, lambda2, dvel_ds)
45+
!$acc routine seq
46+
real(wp), dimension(sys_size), intent(inout) :: L
47+
real(wp), intent(in) :: lambda_factor, lambda2
48+
real(wp), dimension(num_dims), intent(in) :: dvel_ds
49+
integer :: i
50+
51+
do i = momxb + 1, momxe
52+
L(i) = lambda_factor*lambda2*dvel_ds(dir_idx(i - contxe))
53+
end do
54+
end subroutine s_fill_velocity_L
55+
56+
!> Fill advection L variables
57+
pure subroutine s_fill_advection_L(L, lambda_factor, lambda2, dadv_ds)
58+
!$acc routine seq
59+
real(wp), dimension(sys_size), intent(inout) :: L
60+
real(wp), intent(in) :: lambda_factor, lambda2
61+
real(wp), dimension(num_fluids), intent(in) :: dadv_ds
62+
integer :: i
63+
64+
do i = E_idx, advxe - 1
65+
L(i) = lambda_factor*lambda2*dadv_ds(i - momxe)
66+
end do
67+
end subroutine s_fill_advection_L
68+
69+
!> Fill chemistry L variables
70+
pure subroutine s_fill_chemistry_L(L, lambda_factor, lambda2, dYs_ds)
71+
!$acc routine seq
72+
real(wp), dimension(sys_size), intent(inout) :: L
73+
real(wp), intent(in) :: lambda_factor, lambda2
74+
real(wp), dimension(num_species), intent(in) :: dYs_ds
75+
integer :: i
76+
77+
if (.not. chemistry) return
78+
79+
do i = chemxb, chemxe
80+
L(i) = lambda_factor*lambda2*dYs_ds(i - chemxb + 1)
81+
end do
82+
end subroutine s_fill_chemistry_L
83+
84+
!> Slip wall CBC (Thompson 1990, pg. 451)
2685
pure subroutine s_compute_slip_wall_L(lambda, L, rho, c, dpres_ds, dvel_ds)
2786
#ifdef _CRAYFTN
2887
!DIR$ INLINEALWAYS s_compute_slip_wall_L
@@ -31,26 +90,16 @@ contains
3190
#endif
3291
real(wp), dimension(3), intent(in) :: lambda
3392
real(wp), dimension(sys_size), intent(inout) :: L
34-
real(wp), intent(in) :: rho, c
35-
real(wp), intent(in) :: dpres_ds
93+
real(wp), intent(in) :: rho, c, dpres_ds
3694
real(wp), dimension(num_dims), intent(in) :: dvel_ds
37-
3895
integer :: i
3996

40-
L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
41-
42-
do i = 2, advxe
43-
L(i) = 0._wp
44-
end do
45-
97+
L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds)
98+
L(2:advxe - 1) = 0._wp
4699
L(advxe) = L(1)
47-
48100
end subroutine s_compute_slip_wall_L
49101

50-
!> The L variables for the nonreflecting subsonic buffer CBC
51-
!! see pg. 13 of Thompson (1987). The nonreflecting subsonic
52-
!! buffer reduces the amplitude of any reflections caused by
53-
!! outgoing waves.
102+
!> Nonreflecting subsonic buffer CBC (Thompson 1987, pg. 13)
54103
pure subroutine s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds)
55104
#ifdef _CRAYFTN
56105
!DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_buffer_L
@@ -65,42 +114,22 @@ contains
65114
real(wp), dimension(num_dims), intent(in) :: dvel_ds
66115
real(wp), dimension(num_fluids), intent(in) :: dadv_ds
67116
real(wp), dimension(num_species), intent(in) :: dYs_ds
117+
real(wp) :: lambda_factor
68118

69-
integer :: i !< Generic loop iterator
70-
71-
L(1) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(1)))*lambda(1) &
72-
*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
119+
lambda_factor = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(1)))
120+
L(1) = lambda_factor*lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
73121

74-
do i = 2, momxb
75-
L(i) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(2)))*lambda(2) &
76-
*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds)
77-
end do
78-
79-
do i = momxb + 1, momxe
80-
L(i) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(2)))*lambda(2) &
81-
*(dvel_ds(dir_idx(i - contxe)))
82-
end do
83-
84-
do i = E_idx, advxe - 1
85-
L(i) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(2)))*lambda(2) &
86-
*(dadv_ds(i - momxe))
87-
end do
88-
89-
L(advxe) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(3)))*lambda(3) &
90-
*(dpres_ds + rho*c*dvel_ds(dir_idx(1)))
91-
92-
if (chemistry) then
93-
do i = chemxb, chemxe
94-
L(i) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(2)))*lambda(2) &
95-
*(dYs_ds(i - chemxb + 1))
96-
end do
97-
end if
122+
lambda_factor = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(2)))
123+
call s_fill_density_L(L, lambda_factor, lambda(2), c, mf, dalpha_rho_ds, dpres_ds)
124+
call s_fill_velocity_L(L, lambda_factor, lambda(2), dvel_ds)
125+
call s_fill_advection_L(L, lambda_factor, lambda(2), dadv_ds)
126+
call s_fill_chemistry_L(L, lambda_factor, lambda(2), dYs_ds)
98127

128+
lambda_factor = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(3)))
129+
L(advxe) = lambda_factor*lambda(3)*(dpres_ds + rho*c*dvel_ds(dir_idx(1)))
99130
end subroutine s_compute_nonreflecting_subsonic_buffer_L
100-
!> The L variables for the nonreflecting subsonic inflow CBC
101-
!! see pg. 455, Thompson (1990). This nonreflecting subsonic
102-
!! CBC assumes an incoming flow and reduces the amplitude of
103-
!! any reflections caused by outgoing waves.
131+
132+
!> Nonreflecting subsonic inflow CBC (Thompson 1990, pg. 455)
104133
pure subroutine s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, dpres_ds, dvel_ds)
105134
#ifdef _CRAYFTN
106135
!DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_inflow_L
@@ -109,30 +138,15 @@ contains
109138
#endif
110139
real(wp), dimension(3), intent(in) :: lambda
111140
real(wp), dimension(sys_size), intent(inout) :: L
112-
real(wp), intent(in) :: rho, c
113-
real(wp), intent(in) :: dpres_ds
141+
real(wp), intent(in) :: rho, c, dpres_ds
114142
real(wp), dimension(num_dims), intent(in) :: dvel_ds
115143

116-
integer :: i
117-
118-
L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
119-
120-
do i = 2, advxe
121-
L(i) = 0._wp
122-
end do
123-
124-
if (chemistry) then
125-
do i = chemxb, chemxe
126-
L(i) = 0._wp
127-
end do
128-
end if
129-
144+
L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds)
145+
L(2:advxe) = 0._wp
146+
if (chemistry) L(chemxb:chemxe) = 0._wp
130147
end subroutine s_compute_nonreflecting_subsonic_inflow_L
131148

132-
!> The L variables for the nonreflecting subsonic outflow
133-
!! CBC see pg. 454 of Thompson (1990). This nonreflecting
134-
!! subsonic CBC presumes an outgoing flow and reduces the
135-
!! amplitude of any reflections caused by outgoing waves.
149+
!> Nonreflecting subsonic outflow CBC (Thompson 1990, pg. 454)
136150
pure subroutine s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds)
137151
#ifdef _CRAYFTN
138152
!DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_outflow_L
@@ -148,40 +162,15 @@ contains
148162
real(wp), dimension(num_fluids), intent(in) :: dadv_ds
149163
real(wp), dimension(num_species), intent(in) :: dYs_ds
150164

151-
integer :: i !> Generic loop iterator
152-
153-
L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
154-
155-
do i = 2, momxb
156-
L(i) = lambda(2)*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds)
157-
end do
158-
159-
do i = momxb + 1, momxe
160-
L(i) = lambda(2)*(dvel_ds(dir_idx(i - contxe)))
161-
end do
162-
163-
do i = E_idx, advxe - 1
164-
L(i) = lambda(2)*(dadv_ds(i - momxe))
165-
end do
166-
167-
! bubble index
165+
L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds)
166+
call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds)
167+
call s_fill_velocity_L(L, 1._wp, lambda(2), dvel_ds)
168+
call s_fill_advection_L(L, 1._wp, lambda(2), dadv_ds)
169+
call s_fill_chemistry_L(L, 1._wp, lambda(2), dYs_ds)
168170
L(advxe) = 0._wp
169-
170-
if (chemistry) then
171-
do i = chemxb, chemxe
172-
L(i) = lambda(2)*dYs_ds(i - chemxb + 1)
173-
end do
174-
end if
175-
176171
end subroutine s_compute_nonreflecting_subsonic_outflow_L
177172

178-
!> The L variables for the force-free subsonic outflow CBC,
179-
!! see pg. 454 of Thompson (1990). The force-free subsonic
180-
!! outflow sets to zero the sum of all of the forces which
181-
!! are acting on a fluid element for the normal coordinate
182-
!! direction to the boundary. As a result, a fluid element
183-
!! at the boundary is simply advected outward at the fluid
184-
!! velocity.
173+
!> Force-free subsonic outflow CBC (Thompson 1990, pg. 454)
185174
pure subroutine s_compute_force_free_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
186175
#ifdef _CRAYFTN
187176
!DIR$ INLINEALWAYS s_compute_force_free_subsonic_outflow_L
@@ -196,30 +185,14 @@ contains
196185
real(wp), dimension(num_dims), intent(in) :: dvel_ds
197186
real(wp), dimension(num_fluids), intent(in) :: dadv_ds
198187

199-
integer :: i !> Generic loop iterator
200-
201-
L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
202-
203-
do i = 2, momxb
204-
L(i) = lambda(2)*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds)
205-
end do
206-
207-
do i = momxb + 1, momxe
208-
L(i) = lambda(2)*(dvel_ds(dir_idx(i - contxe)))
209-
end do
210-
211-
do i = E_idx, advxe - 1
212-
L(i) = lambda(2)*(dadv_ds(i - momxe))
213-
end do
214-
188+
L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds)
189+
call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds)
190+
call s_fill_velocity_L(L, 1._wp, lambda(2), dvel_ds)
191+
call s_fill_advection_L(L, 1._wp, lambda(2), dadv_ds)
215192
L(advxe) = L(1) + 2._wp*rho*c*lambda(2)*dvel_ds(dir_idx(1))
216-
217193
end subroutine s_compute_force_free_subsonic_outflow_L
218194

219-
!> L variables for the constant pressure subsonic outflow
220-
!! CBC see pg. 455 Thompson (1990). The constant pressure
221-
!! subsonic outflow maintains a fixed pressure at the CBC
222-
!! boundary in absence of any transverse effects.
195+
!> Constant pressure subsonic outflow CBC (Thompson 1990, pg. 455)
223196
pure subroutine s_compute_constant_pressure_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
224197
#ifdef _CRAYFTN
225198
!DIR$ INLINEALWAYS s_compute_constant_pressure_subsonic_outflow_L
@@ -234,57 +207,26 @@ contains
234207
real(wp), dimension(num_dims), intent(in) :: dvel_ds
235208
real(wp), dimension(num_fluids), intent(in) :: dadv_ds
236209

237-
integer :: i !> Generic loop iterator
238-
239-
L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
240-
241-
do i = 2, momxb
242-
L(i) = lambda(2)*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds)
243-
end do
244-
245-
do i = momxb + 1, momxe
246-
L(i) = lambda(2)*(dvel_ds(dir_idx(i - contxe)))
247-
end do
248-
249-
do i = E_idx, advxe - 1
250-
L(i) = lambda(2)*(dadv_ds(i - momxe))
251-
end do
252-
210+
L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds)
211+
call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds)
212+
call s_fill_velocity_L(L, 1._wp, lambda(2), dvel_ds)
213+
call s_fill_advection_L(L, 1._wp, lambda(2), dadv_ds)
253214
L(advxe) = -L(1)
254-
255215
end subroutine s_compute_constant_pressure_subsonic_outflow_L
256216

257-
!> L variables for the supersonic inflow CBC, see pg. 453
258-
!! Thompson (1990). The supersonic inflow CBC is a steady
259-
!! state, or nearly a steady state, CBC in which only the
260-
!! transverse terms may generate a time dependence at the
261-
!! inflow boundary.
217+
!> Supersonic inflow CBC (Thompson 1990, pg. 453)
262218
pure subroutine s_compute_supersonic_inflow_L(L)
263219
#ifdef _CRAYFTN
264220
!DIR$ INLINEALWAYS s_compute_supersonic_inflow_L
265221
#else
266222
!$acc routine seq
267223
#endif
268224
real(wp), dimension(sys_size), intent(inout) :: L
269-
270-
integer :: i
271-
272-
do i = 1, advxe
273-
L(i) = 0._wp
274-
end do
275-
276-
if (chemistry) then
277-
do i = chemxb, chemxe
278-
L(i) = 0._wp
279-
end do
280-
end if
281-
225+
L(1:advxe) = 0._wp
226+
if (chemistry) L(chemxb:chemxe) = 0._wp
282227
end subroutine s_compute_supersonic_inflow_L
283228

284-
!> L variables for the supersonic outflow CBC, see pg. 453
285-
!! of Thompson (1990). For the supersonic outflow CBC, the
286-
!! flow evolution at the boundary is determined completely
287-
!! by the interior data.
229+
!> Supersonic outflow CBC (Thompson 1990, pg. 453)
288230
pure subroutine s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds)
289231
#ifdef _CRAYFTN
290232
!DIR$ INLINEALWAYS s_compute_supersonic_outflow_L
@@ -299,30 +241,12 @@ contains
299241
real(wp), dimension(num_dims), intent(in) :: dvel_ds
300242
real(wp), dimension(num_fluids), intent(in) :: dadv_ds
301243
real(wp), dimension(num_species), intent(in) :: dYs_ds
302-
integer :: i !< Generic loop iterator
303-
304-
L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
305-
306-
do i = 2, momxb
307-
L(i) = lambda(2)*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds)
308-
end do
309-
310-
do i = momxb + 1, momxe
311-
L(i) = lambda(2)*(dvel_ds(dir_idx(i - contxe)))
312-
end do
313-
314-
do i = E_idx, advxe - 1
315-
L(i) = lambda(2)*(dadv_ds(i - momxe))
316-
end do
317244

245+
L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds)
246+
call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds)
247+
call s_fill_velocity_L(L, 1._wp, lambda(2), dvel_ds)
248+
call s_fill_advection_L(L, 1._wp, lambda(2), dadv_ds)
249+
call s_fill_chemistry_L(L, 1._wp, lambda(2), dYs_ds)
318250
L(advxe) = lambda(3)*(dpres_ds + rho*c*dvel_ds(dir_idx(1)))
319-
320-
if (chemistry) then
321-
do i = chemxb, chemxe
322-
L(i) = lambda(2)*dYs_ds(i - chemxb + 1)
323-
end do
324-
end if
325-
326251
end subroutine s_compute_supersonic_outflow_L
327-
328252
end module m_compute_cbc

0 commit comments

Comments
0 (0)
Morty Proxy This is a proxified and sanitized view of the page, visit original site.