Skip to content

Commit 4ac189a

Browse files
sbryngelsonprathi-wind
sbryngelson
authored and
prathi-wind
committed
Make CBC not so terrible (MFlowCode#870)
1 parent 7cb2035 commit 4ac189a

File tree

1 file changed

+112
-188
lines changed

1 file changed

+112
-188
lines changed

‎src/simulation/m_compute_cbc.fpp‎

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