@@ -16,6 +16,65 @@ module m_sim_helpers
1616
1717contains
1818
19+ !> Computes the modified dtheta for Fourier filtering in azimuthal direction
20+ !! @param k y coordinate index
21+ !! @param l z coordinate index
22+ !! @return fltr_dtheta Modified dtheta value for cylindrical coordinates
23+ pure function f_compute_filtered_dtheta (k , l ) result(fltr_dtheta)
24+ !$acc routine seq
25+ integer , intent (in ) :: k, l
26+ real (wp) :: fltr_dtheta
27+ integer :: Nfq
28+
29+ if (grid_geometry == 3 ) then
30+ if (k == 0 ) then
31+ fltr_dtheta = 2._wp * pi* y_cb(0 )/ 3._wp
32+ elseif (k <= fourier_rings) then
33+ Nfq = min (floor(2._wp * real (k, wp)* pi), (p + 1 )/ 2 + 1 )
34+ fltr_dtheta = 2._wp * pi* y_cb(k - 1 )/ real (Nfq, wp)
35+ else
36+ fltr_dtheta = y_cb(k - 1 )* dz(l)
37+ end if
38+ else
39+ fltr_dtheta = 0._wp
40+ end if
41+ end function f_compute_filtered_dtheta
42+
43+ !> Computes inviscid CFL terms for multi- dimensional cases (2D / 3D only)
44+ !! @param vel directional velocities
45+ !! @param c mixture speed of sound
46+ !! @param j x coordinate index
47+ !! @param k y coordinate index
48+ !! @param l z coordinate index
49+ !! @return cfl_terms computed CFL terms for 2D / 3D cases
50+ pure function f_compute_multidim_cfl_terms (vel , c , j , k , l ) result(cfl_terms)
51+ !$acc routine seq
52+ real (wp), dimension (num_vels), intent (in ) :: vel
53+ real (wp), intent (in ) :: c
54+ integer , intent (in ) :: j, k, l
55+ real (wp) :: cfl_terms
56+ real (wp) :: fltr_dtheta
57+
58+ fltr_dtheta = f_compute_filtered_dtheta(k, l)
59+
60+ if (p > 0 ) then
61+ !3D
62+ if (grid_geometry == 3 ) then
63+ cfl_terms = min (dx(j)/ (abs (vel(1 )) + c), &
64+ dy(k)/ (abs (vel(2 )) + c), &
65+ fltr_dtheta/ (abs (vel(3 )) + c))
66+ else
67+ cfl_terms = min (dx(j)/ (abs (vel(1 )) + c), &
68+ dy(k)/ (abs (vel(2 )) + c), &
69+ dz(l)/ (abs (vel(3 )) + c))
70+ end if
71+ else
72+ !2D
73+ cfl_terms = min (dx(j)/ (abs (vel(1 )) + c), &
74+ dy(k)/ (abs (vel(2 )) + c))
75+ end if
76+ end function f_compute_multidim_cfl_terms
77+
1978 !> Computes enthalpy
2079 !! @param q_prim_vf cell centered primitive variables
2180 !! @param pres mixture pressure
@@ -76,7 +135,7 @@ contains
76135
77136 E = gamma* pres + pi_inf + 5.e-1_wp * rho* vel_sum + qv
78137
79- ! ENERGY ADJUSTMENTS FOR HYPERELASTIC ENERGY
138+ ! Adjust energy for hyperelasticity
80139 if (hyperelasticity) then
81140 E = E + G* q_prim_vf(xiend + 1 )%sf(j, k, l)
82141 end if
@@ -92,8 +151,8 @@ contains
92151 !! @param j x index
93152 !! @param k y index
94153 !! @param l z index
95- !! @param icfl_sf cell centered inviscid cfl number
96- !! @param vcfl_sf (optional) cell centered viscous cfl number
154+ !! @param icfl_sf cell- centered inviscid cfl number
155+ !! @param vcfl_sf (optional) cell- centered viscous CFL number
97156 !! @param Rc_sf (optional) cell centered Rc
98157 pure subroutine s_compute_stability_from_dt (vel , c , rho , Re_l , j , k , l , icfl_sf , vcfl_sf , Rc_sf )
99158 $:GPU_ROUTINE(parallelism= ' [seq]' )
@@ -104,82 +163,48 @@ contains
104163 real (wp), dimension (2 ), intent (in ) :: Re_l
105164 integer , intent (in ) :: j, k, l
106165
107- real (wp) :: fltr_dtheta !<
108- !! Modified dtheta accounting for Fourier filtering in azimuthal direction.
109- integer :: Nfq
166+ real (wp) :: fltr_dtheta
110167
111- if (grid_geometry == 3 ) then
112- if (k == 0 ) then
113- fltr_dtheta = 2._wp * pi* y_cb(0 )/ 3._wp
114- elseif (k <= fourier_rings) then
115- Nfq = min (floor(2._wp * real (k, wp)* pi), (p + 1 )/ 2 + 1 )
116- fltr_dtheta = 2._wp * pi* y_cb(k - 1 )/ real (Nfq, wp)
117- else
118- fltr_dtheta = y_cb(k - 1 )* dz(l)
119- end if
168+ ! Inviscid CFL calculation
169+ if (p > 0 .or. n > 0 ) then
170+ ! 2D / 3D
171+ icfl_sf(j, k, l) = dt/ f_compute_multidim_cfl_terms(vel, c, j, k, l)
172+ else
173+ ! 1D
174+ icfl_sf(j, k, l) = (dt/ dx(j))* (abs (vel(1 )) + c)
120175 end if
121176
122- if (p > 0 ) then
123- !3D
124- if (grid_geometry == 3 ) then
125- icfl_sf(j, k, l) = dt/ min (dx(j)/ (abs (vel(1 )) + c), &
126- dy(k)/ (abs (vel(2 )) + c), &
127- fltr_dtheta/ (abs (vel(3 )) + c))
128- else
129- icfl_sf(j, k, l) = dt/ min (dx(j)/ (abs (vel(1 )) + c), &
130- dy(k)/ (abs (vel(2 )) + c), &
131- dz(l)/ (abs (vel(3 )) + c))
132- end if
133-
134- if (viscous) then
135-
177+ ! Viscous calculations
178+ if (viscous) then
179+ if (p > 0 ) then
180+ !3D
136181 if (grid_geometry == 3 ) then
182+ fltr_dtheta = f_compute_filtered_dtheta(k, l)
137183 vcfl_sf(j, k, l) = maxval (dt/ Re_l/ rho) &
138184 / min (dx(j), dy(k), fltr_dtheta)** 2._wp
139-
140185 Rc_sf(j, k, l) = min (dx(j)* (abs (vel(1 )) + c), &
141186 dy(k)* (abs (vel(2 )) + c), &
142187 fltr_dtheta* (abs (vel(3 )) + c)) &
143188 / maxval (1._wp / Re_l)
144189 else
145190 vcfl_sf(j, k, l) = maxval (dt/ Re_l/ rho) &
146191 / min (dx(j), dy(k), dz(l))** 2._wp
147-
148192 Rc_sf(j, k, l) = min (dx(j)* (abs (vel(1 )) + c), &
149193 dy(k)* (abs (vel(2 )) + c), &
150194 dz(l)* (abs (vel(3 )) + c)) &
151195 / maxval (1._wp / Re_l)
152196 end if
153-
154- end if
155-
156- elseif (n > 0 ) then
157- !2D
158- icfl_sf(j, k, l) = dt/ min (dx(j)/ (abs (vel(1 )) + c), &
159- dy(k)/ (abs (vel(2 )) + c))
160-
161- if (viscous) then
162-
197+ elseif (n > 0 ) then
198+ !2D
163199 vcfl_sf(j, k, l) = maxval (dt/ Re_l/ rho)/ min (dx(j), dy(k))** 2._wp
164-
165200 Rc_sf(j, k, l) = min (dx(j)* (abs (vel(1 )) + c), &
166201 dy(k)* (abs (vel(2 )) + c)) &
167202 / maxval (1._wp / Re_l)
168-
169- end if
170-
171- else
172- !1D
173- icfl_sf(j, k, l) = (dt/ dx(j))* (abs (vel(1 )) + c)
174-
175- if (viscous) then
176-
203+ else
204+ !1D
177205 vcfl_sf(j, k, l) = maxval (dt/ Re_l/ rho)/ dx(j)** 2._wp
178-
179206 Rc_sf(j, k, l) = dx(j)* (abs (vel(1 )) + c)/ maxval (1._wp / Re_l)
180-
181207 end if
182-
183208 end if
184209
185210 end subroutine s_compute_stability_from_dt
@@ -201,64 +226,39 @@ contains
201226 integer , intent (in ) :: j, k, l
202227
203228 real (wp) :: icfl_dt, vcfl_dt
204- real (wp) :: fltr_dtheta !<
205- !! Modified dtheta accounting for Fourier filtering in azimuthal direction.
229+ real (wp) :: fltr_dtheta
206230
207- integer :: Nfq
208-
209- if (grid_geometry == 3 ) then
210- if (k == 0 ) then
211- fltr_dtheta = 2._wp * pi* y_cb(0 )/ 3._wp
212- elseif (k <= fourier_rings) then
213- Nfq = min (floor(2._wp * real (k, wp)* pi), (p + 1 )/ 2 + 1 )
214- fltr_dtheta = 2._wp * pi* y_cb(k - 1 )/ real (Nfq, wp)
215- else
216- fltr_dtheta = y_cb(k - 1 )* dz(l)
217- end if
231+ ! Inviscid CFL calculation
232+ if (p > 0 .or. n > 0 ) then
233+ ! 2D / 3D cases
234+ icfl_dt = cfl_target* f_compute_multidim_cfl_terms(vel, c, j, k, l)
235+ else
236+ ! 1D case
237+ icfl_dt = cfl_target* (dx(j)/ (abs (vel(1 )) + c))
218238 end if
219239
220- if (p > 0 ) then
221- !3D
222- if (grid_geometry == 3 ) then
223- icfl_dt = cfl_target* min (dx(j)/ (abs (vel(1 )) + c), &
224- dy(k)/ (abs (vel(2 )) + c), &
225- fltr_dtheta/ (abs (vel(3 )) + c))
226- else
227- icfl_dt = cfl_target* min (dx(j)/ (abs (vel(1 )) + c), &
228- dy(k)/ (abs (vel(2 )) + c), &
229- dz(l)/ (abs (vel(3 )) + c))
230- end if
231-
232- if (viscous) then
240+ ! Viscous calculations
241+ if (viscous) then
242+ if (p > 0 ) then
243+ !3D
233244 if (grid_geometry == 3 ) then
245+ fltr_dtheta = f_compute_filtered_dtheta(k, l)
234246 vcfl_dt = cfl_target* (min (dx(j), dy(k), fltr_dtheta)** 2._wp ) &
235247 / minval (1 / (rho* Re_l))
236248 else
237249 vcfl_dt = cfl_target* (min (dx(j), dy(k), dz(l))** 2._wp ) &
238250 / minval (1 / (rho* Re_l))
239251 end if
240- end if
241-
242- elseif (n > 0 ) then
243- !2D
244- icfl_dt = cfl_target* min (dx(j)/ (abs (vel(1 )) + c), &
245- dy(k)/ (abs (vel(2 )) + c))
246-
247- if (viscous) then
252+ elseif (n > 0 ) then
253+ !2D
248254 vcfl_dt = cfl_target* (min (dx(j), dy(k))** 2._wp )/ maxval ((1 / Re_l)/ rho)
249- end if
250-
251- else
252- !1D
253- icfl_dt = cfl_target* (dx(j)/ (abs (vel(1 )) + c))
254-
255- if (viscous) then
255+ else
256+ !1D
256257 vcfl_dt = cfl_target* (dx(j)** 2._wp )/ minval (1 / (rho* Re_l))
257258 end if
258-
259259 end if
260260
261- if (any (re_size > 0 )) then
261+ if (any (Re_size > 0 )) then
262262 max_dt(j, k, l) = min (icfl_dt, vcfl_dt)
263263 else
264264 max_dt(j, k, l) = icfl_dt
0 commit comments