Skip to content

Commit 1ab0f5e

Browse files
sbryngelsonsbryngelson
authored
Refac helpers for CFL conditions (#901)
1 parent 07d8f1f commit 1ab0f5e

File tree

2 files changed

+103
-99
lines changed

2 files changed

+103
-99
lines changed

‎.github/workflows/bench.yml

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
name: 'Benchmark'
22

3-
on:
3+
on:
4+
pull_request:
45
pull_request_review:
56
types: [submitted]
67
workflow_dispatch:
@@ -23,7 +24,10 @@ jobs:
2324

2425
self:
2526
name: Georgia Tech | Phoenix (NVHPC)
26-
if: github.repository == 'MFlowCode/MFC' && needs.file-changes.outputs.checkall == 'true' && ${{ github.event.review.state == 'approved' }}
27+
if: ${{ github.repository == 'MFlowCode/MFC' && needs.file-changes.outputs.checkall == 'true' && (
28+
(github.event_name == 'pull_request_review' && github.event.review.state == 'approved') ||
29+
(github.event_name == 'pull_request' && github.event.pull_request.user.login == 'sbryngelson')
30+
) }}
2731
needs: file-changes
2832
strategy:
2933
matrix:

‎src/simulation/m_sim_helpers.f90

Lines changed: 97 additions & 97 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,65 @@ module m_sim_helpers
1414

1515
contains
1616

17+
!> Computes the modified dtheta for Fourier filtering in azimuthal direction
18+
!! @param k y coordinate index
19+
!! @param l z coordinate index
20+
!! @return fltr_dtheta Modified dtheta value for cylindrical coordinates
21+
pure function f_compute_filtered_dtheta(k, l) result(fltr_dtheta)
22+
!$acc routine seq
23+
integer, intent(in) :: k, l
24+
real(wp) :: fltr_dtheta
25+
integer :: Nfq
26+
27+
if (grid_geometry == 3) then
28+
if (k == 0) then
29+
fltr_dtheta = 2._wp*pi*y_cb(0)/3._wp
30+
elseif (k <= fourier_rings) then
31+
Nfq = min(floor(2._wp*real(k, wp)*pi), (p + 1)/2 + 1)
32+
fltr_dtheta = 2._wp*pi*y_cb(k - 1)/real(Nfq, wp)
33+
else
34+
fltr_dtheta = y_cb(k - 1)*dz(l)
35+
end if
36+
else
37+
fltr_dtheta = 0._wp
38+
end if
39+
end function f_compute_filtered_dtheta
40+
41+
!> Computes inviscid CFL terms for multi-dimensional cases (2D/3D only)
42+
!! @param vel directional velocities
43+
!! @param c mixture speed of sound
44+
!! @param j x coordinate index
45+
!! @param k y coordinate index
46+
!! @param l z coordinate index
47+
!! @return cfl_terms computed CFL terms for 2D/3D cases
48+
pure function f_compute_multidim_cfl_terms(vel, c, j, k, l) result(cfl_terms)
49+
!$acc routine seq
50+
real(wp), dimension(num_vels), intent(in) :: vel
51+
real(wp), intent(in) :: c
52+
integer, intent(in) :: j, k, l
53+
real(wp) :: cfl_terms
54+
real(wp) :: fltr_dtheta
55+
56+
fltr_dtheta = f_compute_filtered_dtheta(k, l)
57+
58+
if (p > 0) then
59+
!3D
60+
if (grid_geometry == 3) then
61+
cfl_terms = min(dx(j)/(abs(vel(1)) + c), &
62+
dy(k)/(abs(vel(2)) + c), &
63+
fltr_dtheta/(abs(vel(3)) + c))
64+
else
65+
cfl_terms = min(dx(j)/(abs(vel(1)) + c), &
66+
dy(k)/(abs(vel(2)) + c), &
67+
dz(l)/(abs(vel(3)) + c))
68+
end if
69+
else
70+
!2D
71+
cfl_terms = min(dx(j)/(abs(vel(1)) + c), &
72+
dy(k)/(abs(vel(2)) + c))
73+
end if
74+
end function f_compute_multidim_cfl_terms
75+
1776
!> Computes enthalpy
1877
!! @param q_prim_vf cell centered primitive variables
1978
!! @param pres mixture pressure
@@ -77,7 +136,7 @@ pure subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, a
77136

78137
E = gamma*pres + pi_inf + 5.e-1_wp*rho*vel_sum + qv
79138

80-
! ENERGY ADJUSTMENTS FOR HYPERELASTIC ENERGY
139+
! Adjust energy for hyperelasticity
81140
if (hyperelasticity) then
82141
E = E + G*q_prim_vf(xiend + 1)%sf(j, k, l)
83142
end if
@@ -93,8 +152,8 @@ end subroutine s_compute_enthalpy
93152
!! @param j x index
94153
!! @param k y index
95154
!! @param l z index
96-
!! @param icfl_sf cell centered inviscid cfl number
97-
!! @param vcfl_sf (optional) cell centered viscous cfl number
155+
!! @param icfl_sf cell-centered inviscid cfl number
156+
!! @param vcfl_sf (optional) cell-centered viscous CFL number
98157
!! @param Rc_sf (optional) cell centered Rc
99158
pure subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl_sf, Rc_sf)
100159
!$acc routine seq
@@ -105,82 +164,48 @@ pure subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf,
105164
real(wp), dimension(2), intent(in) :: Re_l
106165
integer, intent(in) :: j, k, l
107166

108-
real(wp) :: fltr_dtheta !<
109-
!! Modified dtheta accounting for Fourier filtering in azimuthal direction.
110-
integer :: Nfq
167+
real(wp) :: fltr_dtheta
111168

112-
if (grid_geometry == 3) then
113-
if (k == 0) then
114-
fltr_dtheta = 2._wp*pi*y_cb(0)/3._wp
115-
elseif (k <= fourier_rings) then
116-
Nfq = min(floor(2._wp*real(k, wp)*pi), (p + 1)/2 + 1)
117-
fltr_dtheta = 2._wp*pi*y_cb(k - 1)/real(Nfq, wp)
118-
else
119-
fltr_dtheta = y_cb(k - 1)*dz(l)
120-
end if
169+
! Inviscid CFL calculation
170+
if (p > 0 .or. n > 0) then
171+
! 2D/3D
172+
icfl_sf(j, k, l) = dt/f_compute_multidim_cfl_terms(vel, c, j, k, l)
173+
else
174+
! 1D
175+
icfl_sf(j, k, l) = (dt/dx(j))*(abs(vel(1)) + c)
121176
end if
122177

123-
if (p > 0) then
124-
!3D
125-
if (grid_geometry == 3) then
126-
icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), &
127-
dy(k)/(abs(vel(2)) + c), &
128-
fltr_dtheta/(abs(vel(3)) + c))
129-
else
130-
icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), &
131-
dy(k)/(abs(vel(2)) + c), &
132-
dz(l)/(abs(vel(3)) + c))
133-
end if
134-
135-
if (viscous) then
136-
178+
! Viscous calculations
179+
if (viscous) then
180+
if (p > 0) then
181+
!3D
137182
if (grid_geometry == 3) then
183+
fltr_dtheta = f_compute_filtered_dtheta(k, l)
138184
vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) &
139185
/min(dx(j), dy(k), fltr_dtheta)**2._wp
140-
141186
Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), &
142187
dy(k)*(abs(vel(2)) + c), &
143188
fltr_dtheta*(abs(vel(3)) + c)) &
144189
/maxval(1._wp/Re_l)
145190
else
146191
vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) &
147192
/min(dx(j), dy(k), dz(l))**2._wp
148-
149193
Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), &
150194
dy(k)*(abs(vel(2)) + c), &
151195
dz(l)*(abs(vel(3)) + c)) &
152196
/maxval(1._wp/Re_l)
153197
end if
154-
155-
end if
156-
157-
elseif (n > 0) then
158-
!2D
159-
icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), &
160-
dy(k)/(abs(vel(2)) + c))
161-
162-
if (viscous) then
163-
198+
elseif (n > 0) then
199+
!2D
164200
vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k))**2._wp
165-
166201
Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), &
167202
dy(k)*(abs(vel(2)) + c)) &
168203
/maxval(1._wp/Re_l)
169-
170-
end if
171-
172-
else
173-
!1D
174-
icfl_sf(j, k, l) = (dt/dx(j))*(abs(vel(1)) + c)
175-
176-
if (viscous) then
177-
204+
else
205+
!1D
178206
vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/dx(j)**2._wp
179-
180207
Rc_sf(j, k, l) = dx(j)*(abs(vel(1)) + c)/maxval(1._wp/Re_l)
181-
182208
end if
183-
184209
end if
185210

186211
end subroutine s_compute_stability_from_dt
@@ -202,64 +227,39 @@ pure subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l)
202227
integer, intent(in) :: j, k, l
203228

204229
real(wp) :: icfl_dt, vcfl_dt
205-
real(wp) :: fltr_dtheta !<
206-
!! Modified dtheta accounting for Fourier filtering in azimuthal direction.
207-
208-
integer :: Nfq
230+
real(wp) :: fltr_dtheta
209231

210-
if (grid_geometry == 3) then
211-
if (k == 0) then
212-
fltr_dtheta = 2._wp*pi*y_cb(0)/3._wp
213-
elseif (k <= fourier_rings) then
214-
Nfq = min(floor(2._wp*real(k, wp)*pi), (p + 1)/2 + 1)
215-
fltr_dtheta = 2._wp*pi*y_cb(k - 1)/real(Nfq, wp)
216-
else
217-
fltr_dtheta = y_cb(k - 1)*dz(l)
218-
end if
232+
! Inviscid CFL calculation
233+
if (p > 0 .or. n > 0) then
234+
! 2D/3D cases
235+
icfl_dt = cfl_target*f_compute_multidim_cfl_terms(vel, c, j, k, l)
236+
else
237+
! 1D case
238+
icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c))
219239
end if
220240

221-
if (p > 0) then
222-
!3D
223-
if (grid_geometry == 3) then
224-
icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), &
225-
dy(k)/(abs(vel(2)) + c), &
226-
fltr_dtheta/(abs(vel(3)) + c))
227-
else
228-
icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), &
229-
dy(k)/(abs(vel(2)) + c), &
230-
dz(l)/(abs(vel(3)) + c))
231-
end if
232-
233-
if (viscous) then
241+
! Viscous calculations
242+
if (viscous) then
243+
if (p > 0) then
244+
!3D
234245
if (grid_geometry == 3) then
246+
fltr_dtheta = f_compute_filtered_dtheta(k, l)
235247
vcfl_dt = cfl_target*(min(dx(j), dy(k), fltr_dtheta)**2._wp) &
236248
/minval(1/(rho*Re_l))
237249
else
238250
vcfl_dt = cfl_target*(min(dx(j), dy(k), dz(l))**2._wp) &
239251
/minval(1/(rho*Re_l))
240252
end if
241-
end if
242-
243-
elseif (n > 0) then
244-
!2D
245-
icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), &
246-
dy(k)/(abs(vel(2)) + c))
247-
248-
if (viscous) then
253+
elseif (n > 0) then
254+
!2D
249255
vcfl_dt = cfl_target*(min(dx(j), dy(k))**2._wp)/maxval((1/Re_l)/rho)
250-
end if
251-
252-
else
253-
!1D
254-
icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c))
255-
256-
if (viscous) then
256+
else
257+
!1D
257258
vcfl_dt = cfl_target*(dx(j)**2._wp)/minval(1/(rho*Re_l))
258259
end if
259-
260260
end if
261261

262-
if (any(re_size > 0)) then
262+
if (any(Re_size > 0)) then
263263
max_dt(j, k, l) = min(icfl_dt, vcfl_dt)
264264
else
265265
max_dt(j, k, l) = icfl_dt

0 commit comments

Comments
 (0)