Skip to content

Commit 46d55d8

Browse files
sbryngelsonprathi-wind
sbryngelson
authored and
prathi-wind
committed
Make CBC not so terrible (MFlowCode#870)
1 parent 28d9700 commit 46d55d8

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,12 +1,10 @@
11
#:include 'directive_macros.fpp'
22
!>
33
!! @file m_compute_cbc.f90
4-
!! @brief Contains module m_compute_cbc
4+
!! @brief CBC computation module
55

66
module m_compute_cbc
7-
8-
use m_global_parameters !< Definitions of the global parameters
9-
7+
use m_global_parameters
108
implicit none
119

1210
private; public :: s_compute_slip_wall_L, &
@@ -19,11 +17,72 @@ module m_compute_cbc
1917
s_compute_supersonic_outflow_L
2018

2119
contains
20+
!> Base L1 calculation
21+
pure function f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) result(L1)
22+
!$acc routine seq
23+
real(wp), dimension(3), intent(in) :: lambda
24+
real(wp), intent(in) :: rho, c, dpres_ds
25+
real(wp), dimension(num_dims), intent(in) :: dvel_ds
26+
real(wp) :: L1
27+
L1 = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
28+
end function f_base_L1
29+
30+
!> Fill density L variables
31+
pure subroutine s_fill_density_L(L, lambda_factor, lambda2, c, mf, dalpha_rho_ds, dpres_ds)
32+
!$acc routine seq
33+
real(wp), dimension(sys_size), intent(inout) :: L
34+
real(wp), intent(in) :: lambda_factor, lambda2, c
35+
real(wp), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds
36+
real(wp), intent(in) :: dpres_ds
37+
integer :: i
38+
39+
do i = 2, momxb
40+
L(i) = lambda_factor*lambda2*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds)
41+
end do
42+
end subroutine s_fill_density_L
43+
44+
!> Fill velocity L variables
45+
pure subroutine s_fill_velocity_L(L, lambda_factor, lambda2, dvel_ds)
46+
!$acc routine seq
47+
real(wp), dimension(sys_size), intent(inout) :: L
48+
real(wp), intent(in) :: lambda_factor, lambda2
49+
real(wp), dimension(num_dims), intent(in) :: dvel_ds
50+
integer :: i
51+
52+
do i = momxb + 1, momxe
53+
L(i) = lambda_factor*lambda2*dvel_ds(dir_idx(i - contxe))
54+
end do
55+
end subroutine s_fill_velocity_L
56+
57+
!> Fill advection L variables
58+
pure subroutine s_fill_advection_L(L, lambda_factor, lambda2, dadv_ds)
59+
!$acc routine seq
60+
real(wp), dimension(sys_size), intent(inout) :: L
61+
real(wp), intent(in) :: lambda_factor, lambda2
62+
real(wp), dimension(num_fluids), intent(in) :: dadv_ds
63+
integer :: i
64+
65+
do i = E_idx, advxe - 1
66+
L(i) = lambda_factor*lambda2*dadv_ds(i - momxe)
67+
end do
68+
end subroutine s_fill_advection_L
69+
70+
!> Fill chemistry L variables
71+
pure subroutine s_fill_chemistry_L(L, lambda_factor, lambda2, dYs_ds)
72+
!$acc routine seq
73+
real(wp), dimension(sys_size), intent(inout) :: L
74+
real(wp), intent(in) :: lambda_factor, lambda2
75+
real(wp), dimension(num_species), intent(in) :: dYs_ds
76+
integer :: i
77+
78+
if (.not. chemistry) return
79+
80+
do i = chemxb, chemxe
81+
L(i) = lambda_factor*lambda2*dYs_ds(i - chemxb + 1)
82+
end do
83+
end subroutine s_fill_chemistry_L
2284

23-
!> The L variables for the slip wall CBC, see pg. 451 of
24-
!! Thompson (1990). At the slip wall (frictionless wall),
25-
!! the normal component of velocity is zero at all times,
26-
!! while the transverse velocities may be nonzero.
85+
!> Slip wall CBC (Thompson 1990, pg. 451)
2786
pure subroutine s_compute_slip_wall_L(lambda, L, rho, c, dpres_ds, dvel_ds)
2887
#ifdef _CRAYFTN
2988
!DIR$ INLINEALWAYS s_compute_slip_wall_L
@@ -32,26 +91,16 @@ contains
3291
#endif
3392
real(wp), dimension(3), intent(in) :: lambda
3493
real(wp), dimension(sys_size), intent(inout) :: L
35-
real(wp), intent(in) :: rho, c
36-
real(wp), intent(in) :: dpres_ds
94+
real(wp), intent(in) :: rho, c, dpres_ds
3795
real(wp), dimension(num_dims), intent(in) :: dvel_ds
38-
3996
integer :: i
4097

41-
L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
42-
43-
do i = 2, advxe
44-
L(i) = 0._wp
45-
end do
46-
98+
L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds)
99+
L(2:advxe - 1) = 0._wp
47100
L(advxe) = L(1)
48-
49101
end subroutine s_compute_slip_wall_L
50102

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

70-
integer :: i !< Generic loop iterator
120+
lambda_factor = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(1)))
121+
L(1) = lambda_factor*lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1)))
71122

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

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

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

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

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

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

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

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

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

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

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

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

0 commit comments

Comments
 (0)