Skip to content

Commit b007e50

Browse files
sbryngelsonsbryngelson
committed
fixes
1 parent 47fa3d7 commit b007e50

File tree

1 file changed

+76
-54
lines changed

1 file changed

+76
-54
lines changed

‎src/simulation/m_sim_helpers.f90‎

Lines changed: 76 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -14,29 +14,16 @@ module m_sim_helpers
1414

1515
contains
1616

17-
!> Computes CFL-related terms (handles both Fourier filtering and inviscid CFL)
18-
!! This function consolidates the repeated CFL calculation patterns found in both
19-
!! s_compute_stability_from_dt and s_compute_dt_from_cfl subroutines.
20-
!! It handles Fourier filtering for cylindrical coordinates and computes inviscid
21-
!! CFL terms for 1D, 2D, and 3D cases.
22-
!! @param vel directional velocities
23-
!! @param c mixture speed of sound
24-
!! @param j x coordinate index
25-
!! @param k y coordinate index
17+
!> Computes the modified dtheta for Fourier filtering in azimuthal direction
18+
!! @param k y coordinate index
2619
!! @param l z coordinate index
27-
!! @param for_dt_calc logical flag: true for dt calculation, false for stability calculation
28-
!! @return cfl_terms computed CFL terms (with appropriate scaling applied)
29-
pure function f_compute_cfl_terms(vel, c, j, k, l, for_dt_calc) result(cfl_terms)
20+
!! @return fltr_dtheta Modified dtheta value for cylindrical coordinates
21+
pure function f_compute_filtered_dtheta(k, l) result(fltr_dtheta)
3022
!$acc routine seq
31-
real(wp), dimension(num_vels), intent(in) :: vel
32-
real(wp), intent(in) :: c
33-
integer, intent(in) :: j, k, l
34-
logical, intent(in) :: for_dt_calc
35-
real(wp) :: cfl_terms
23+
integer, intent(in) :: k, l
3624
real(wp) :: fltr_dtheta
3725
integer :: Nfq
3826

39-
! Compute filtered dtheta for cylindrical coordinates
4027
if (grid_geometry == 3) then
4128
if (k == 0) then
4229
fltr_dtheta = 2._wp*pi*y_cb(0)/3._wp
@@ -46,38 +33,45 @@ pure function f_compute_cfl_terms(vel, c, j, k, l, for_dt_calc) result(cfl_terms
4633
else
4734
fltr_dtheta = y_cb(k - 1)*dz(l)
4835
end if
36+
else
37+
fltr_dtheta = 0._wp
4938
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)
5057

51-
! Compute CFL terms based on dimensionality
5258
if (p > 0) then
5359
!3D
5460
if (grid_geometry == 3) then
5561
cfl_terms = min(dx(j)/(abs(vel(1)) + c), &
56-
dy(k)/(abs(vel(2)) + c), &
57-
fltr_dtheta/(abs(vel(3)) + c))
62+
dy(k)/(abs(vel(2)) + c), &
63+
fltr_dtheta/(abs(vel(3)) + c))
5864
else
5965
cfl_terms = min(dx(j)/(abs(vel(1)) + c), &
60-
dy(k)/(abs(vel(2)) + c), &
61-
dz(l)/(abs(vel(3)) + c))
66+
dy(k)/(abs(vel(2)) + c), &
67+
dz(l)/(abs(vel(3)) + c))
6268
end if
63-
elseif (n > 0) then
69+
else
6470
!2D
6571
cfl_terms = min(dx(j)/(abs(vel(1)) + c), &
66-
dy(k)/(abs(vel(2)) + c))
67-
else
68-
!1D - special handling for different calculation types
69-
if (for_dt_calc) then
70-
cfl_terms = dx(j)/(abs(vel(1)) + c)
71-
else
72-
cfl_terms = (1._wp/dx(j))*(abs(vel(1)) + c)
73-
end if
74-
end if
75-
76-
! Apply appropriate factor for stability vs dt calculation
77-
if (.not. for_dt_calc .and. (p > 0 .or. n > 0)) then
78-
cfl_terms = 1._wp/cfl_terms
72+
dy(k)/(abs(vel(2)) + c))
7973
end if
80-
end function f_compute_cfl_terms
74+
end function f_compute_multidim_cfl_terms
8175

8276
!> Computes enthalpy
8377
!! @param q_prim_vf cell centered primitive variables
@@ -170,26 +164,44 @@ pure subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf,
170164
real(wp), dimension(2), intent(in) :: Re_l
171165
integer, intent(in) :: j, k, l
172166

167+
real(wp) :: fltr_dtheta
168+
173169
! Inviscid CFL calculation
174-
icfl_sf(j, k, l) = dt*f_compute_cfl_terms(vel, c, j, k, l, .false.)
170+
if (p > 0 .or. n > 0) then
171+
! 2D/3D cases
172+
icfl_sf(j, k, l) = dt/f_compute_multidim_cfl_terms(vel, c, j, k, l)
173+
else
174+
! 1D case - exact original logic
175+
icfl_sf(j, k, l) = (dt/dx(j))*(abs(vel(1)) + c)
176+
end if
175177

176-
! Viscous calculations (simplified with common patterns)
178+
! Viscous calculations
177179
if (viscous) then
180+
fltr_dtheta = f_compute_filtered_dtheta(k, l)
181+
178182
if (p > 0) then
179183
!3D
180184
if (grid_geometry == 3) then
181-
vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k), &
182-
2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp)))**2._wp
183-
Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), dy(k)*(abs(vel(2)) + c), &
184-
2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp))*(abs(vel(3)) + c))/maxval(1._wp/Re_l)
185+
vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) &
186+
/min(dx(j), dy(k), fltr_dtheta)**2._wp
187+
Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), &
188+
dy(k)*(abs(vel(2)) + c), &
189+
fltr_dtheta*(abs(vel(3)) + c)) &
190+
/maxval(1._wp/Re_l)
185191
else
186-
vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k), dz(l))**2._wp
187-
Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), dy(k)*(abs(vel(2)) + c), dz(l)*(abs(vel(3)) + c))/maxval(1._wp/Re_l)
192+
vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) &
193+
/min(dx(j), dy(k), dz(l))**2._wp
194+
Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), &
195+
dy(k)*(abs(vel(2)) + c), &
196+
dz(l)*(abs(vel(3)) + c)) &
197+
/maxval(1._wp/Re_l)
188198
end if
189199
elseif (n > 0) then
190200
!2D
191201
vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k))**2._wp
192-
Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), dy(k)*(abs(vel(2)) + c))/maxval(1._wp/Re_l)
202+
Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), &
203+
dy(k)*(abs(vel(2)) + c)) &
204+
/maxval(1._wp/Re_l)
193205
else
194206
!1D
195207
vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/dx(j)**2._wp
@@ -216,26 +228,36 @@ pure subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l)
216228
integer, intent(in) :: j, k, l
217229

218230
real(wp) :: icfl_dt, vcfl_dt
231+
real(wp) :: fltr_dtheta
219232

220233
! Inviscid CFL calculation
221-
icfl_dt = cfl_target*f_compute_cfl_terms(vel, c, j, k, l, .true.)
234+
if (p > 0 .or. n > 0) then
235+
! 2D/3D cases
236+
icfl_dt = cfl_target*f_compute_multidim_cfl_terms(vel, c, j, k, l)
237+
else
238+
! 1D case - exact original logic
239+
icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c))
240+
end if
222241

223-
! Viscous calculations (simplified with common patterns)
242+
! Viscous calculations
224243
if (viscous) then
244+
fltr_dtheta = f_compute_filtered_dtheta(k, l)
245+
225246
if (p > 0) then
226247
!3D
227248
if (grid_geometry == 3) then
228-
vcfl_dt = cfl_target*min(dx(j), dy(k), &
229-
2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp)))**2._wp/minval(1/(rho*Re_l))
249+
vcfl_dt = cfl_target*(min(dx(j), dy(k), fltr_dtheta)**2._wp) &
250+
/minval(1/(rho*Re_l))
230251
else
231-
vcfl_dt = cfl_target*min(dx(j), dy(k), dz(l))**2._wp/minval(1/(rho*Re_l))
252+
vcfl_dt = cfl_target*(min(dx(j), dy(k), dz(l))**2._wp) &
253+
/minval(1/(rho*Re_l))
232254
end if
233255
elseif (n > 0) then
234256
!2D
235-
vcfl_dt = cfl_target*min(dx(j), dy(k))**2._wp/maxval((1/Re_l)/rho)
257+
vcfl_dt = cfl_target*(min(dx(j), dy(k))**2._wp)/maxval((1/Re_l)/rho)
236258
else
237259
!1D
238-
vcfl_dt = cfl_target*dx(j)**2._wp/minval(1/(rho*Re_l))
260+
vcfl_dt = cfl_target*(dx(j)**2._wp)/minval(1/(rho*Re_l))
239261
end if
240262
end if
241263

0 commit comments

Comments
 (0)