Skip to content

Commit a850436

Browse files
wilfonbawilfonba
committed
add ccfl code
1 parent 17ba8e0 commit a850436

3 files changed

Lines changed: 78 additions & 19 deletions

File tree

‎src/common/m_mpi_common.fpp‎

Lines changed: 12 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -300,7 +300,7 @@ contains
300300
!! single process, within its assigned section of the computational domain. Finally, note that the global extrema values are
301301
!! only bookkeept on the rank 0 processor.
302302
impure subroutine s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, vcfl_max_loc, Rc_min_loc, bubs_loc, icfl_max_glb, &
303-
& vcfl_max_glb, Rc_min_glb, bubs_glb)
303+
& vcfl_max_glb, Rc_min_glb, bubs_glb, ccfl_max_loc, ccfl_max_glb)
304304

305305
real(wp), intent(in) :: icfl_max_loc
306306
real(wp), intent(in) :: vcfl_max_loc
@@ -310,10 +310,13 @@ contains
310310
real(wp), intent(out) :: vcfl_max_glb
311311
real(wp), intent(out) :: Rc_min_glb
312312
integer, intent(out) :: bubs_glb
313+
real(wp), intent(in) :: ccfl_max_loc
314+
real(wp), intent(out) :: ccfl_max_glb
313315

314316
icfl_max_glb = icfl_max_loc
315317
vcfl_max_glb = vcfl_max_loc
316318
Rc_min_glb = Rc_min_loc
319+
ccfl_max_glb = ccfl_max_loc
317320

318321
#ifdef MFC_SIMULATION
319322
#ifdef MFC_MPI
@@ -328,6 +331,10 @@ contains
328331
call MPI_REDUCE(Rc_min_loc, Rc_min_glb, 1, mpi_p, MPI_MIN, 0, MPI_COMM_WORLD, ierr)
329332
end if
330333

334+
if (surface_tension) then
335+
call MPI_REDUCE(ccfl_max_loc, ccfl_max_glb, 1, mpi_p, MPI_MAX, 0, MPI_COMM_WORLD, ierr)
336+
end if
337+
331338
if (bubbles_lagrange) then
332339
call MPI_REDUCE(bubs_loc, bubs_glb, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ierr)
333340
end if
@@ -341,6 +348,10 @@ contains
341348
Rc_min_glb = Rc_min_loc
342349
end if
343350

351+
if (surface_tension) then
352+
ccfl_max_glb = ccfl_max_loc
353+
end if
354+
344355
if (bubbles_lagrange) bubs_glb = bubs_loc
345356
#endif
346357
#endif

‎src/simulation/m_data_output.fpp‎

Lines changed: 25 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -32,10 +32,11 @@ module m_data_output
3232
real(wp), public, allocatable, dimension(:,:) :: c_mass
3333
$:GPU_DECLARE(create='[c_mass]')
3434

35-
!> @name ICFL, VCFL, and Rc stability criteria extrema over all the time-steps
35+
!> @name ICFL, VCFL, CCFL, and Rc stability criteria extrema over all the time-steps
3636
!> @{
3737
real(wp) :: icfl_max !< ICFL criterion maximum
3838
real(wp) :: vcfl_max !< VCFL criterion maximum
39+
real(wp) :: ccfl_max !< CCFL criterion maximum
3940
real(wp) :: Rc_min !< Rc criterion maximum
4041
!> @}
4142

@@ -87,6 +88,10 @@ contains
8788

8889
write (3, '(13X,A9,13X,A10,13X,A10,13X,A10)', advance="no") trim('Time-step'), trim('dt'), trim('Time'), trim('ICFL Max')
8990

91+
if (surface_tension) then
92+
write (3, '(13X,A10)', advance="no") trim('CCFL Max')
93+
end if
94+
9095
if (viscous) then
9196
write (3, '(13X,A10,13X,A16)', advance="no") trim('VCFL Max'), trim('Rc Min')
9297
end if
@@ -182,15 +187,16 @@ contains
182187
real(wp) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids
183188
real(wp) :: ccfl_max_loc, ccfl_max_glb !< CCFL stability extrema on local and global grids
184189
real(wp) :: Rc_min_loc, Rc_min_glb !< Rc stability extrema on local and global grids
185-
real(wp) :: icfl, vcfl, Rc
190+
real(wp) :: icfl, vcfl, ccfl, Rc
186191
integer :: fl !< Fluid loop iterator
187192

188193
icfl_max_loc = 0._wp
189194
vcfl_max_loc = 0._wp
195+
ccfl_max_loc = 0._wp
190196
Rc_min_loc = huge(1.0_wp)
191197
! Computing Stability Criteria at Current Time-step
192198
$:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, icfl, &
193-
& vcfl, Rc, fl]', reduction='[[icfl_max_loc, vcfl_max_loc], [Rc_min_loc]]', reductionOp='[max, min]')
199+
& vcfl, Rc, ccfl, fl]', reduction='[[icfl_max_loc, vcfl_max_loc, ccfl_max_loc], [Rc_min_loc]]', reductionOp='[max, min]')
194200
do l = 0, p
195201
do k = 0, n
196202
do j = 0, m
@@ -210,14 +216,11 @@ contains
210216
Re(1) = 1._wp/max(Re(1), sgm_eps)
211217
end if
212218

213-
if (viscous) then
214-
call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl, vcfl, Rc)
215-
else
216-
call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl)
217-
end if
219+
call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl, vcfl, Rc, ccfl)
218220

219221
icfl_max_loc = max(icfl_max_loc, icfl)
220222
vcfl_max_loc = max(vcfl_max_loc, merge(vcfl, 0.0_wp, viscous))
223+
ccfl_max_loc = max(ccfl_max_loc, merge(ccfl, 0.0_wp, surface_tension))
221224
Rc_min_loc = min(Rc_min_loc, merge(Rc, huge(1.0_wp), viscous))
222225
end do
223226
end do
@@ -227,16 +230,21 @@ contains
227230

228231
if (num_procs > 1) then
229232
call s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, vcfl_max_loc, Rc_min_loc, n_el_bubs_loc, icfl_max_glb, &
230-
& vcfl_max_glb, Rc_min_glb, n_el_bubs_glb)
233+
& vcfl_max_glb, Rc_min_glb, n_el_bubs_glb, ccfl_max_loc, ccfl_max_glb)
231234
else
232235
icfl_max_glb = icfl_max_loc
233236
if (viscous) vcfl_max_glb = vcfl_max_loc
234237
if (viscous) Rc_min_glb = Rc_min_loc
238+
if (surface_tension) ccfl_max_glb = ccfl_max_loc
235239
if (bubbles_lagrange) n_el_bubs_glb = n_el_bubs_loc
236240
end if
237241

238242
if (icfl_max_glb > icfl_max) icfl_max = icfl_max_glb
239243

244+
if (surface_tension) then
245+
if (ccfl_max_glb > ccfl_max) ccfl_max = ccfl_max_glb
246+
end if
247+
240248
if (viscous) then
241249
if (vcfl_max_glb > vcfl_max) vcfl_max = vcfl_max_glb
242250
if (Rc_min_glb < Rc_min) Rc_min = Rc_min_glb
@@ -245,6 +253,10 @@ contains
245253
if (proc_rank == 0) then
246254
write (3, '(13X,I9,13X,F10.6,13X,F10.6,13X,F10.6)', advance="no") t_step, dt, mytime, icfl_max_glb
247255

256+
if (surface_tension) then
257+
write (3, '(13X,F10.6)', advance="no") ccfl_max_glb
258+
end if
259+
248260
if (viscous) then
249261
write (3, '(13X,F10.6,13X,ES16.6)', advance="no") vcfl_max_glb, Rc_min_glb
250262
end if
@@ -1642,6 +1654,7 @@ contains
16421654
write (3, '(A)') ''
16431655

16441656
write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max
1657+
if (surface_tension) write (3, '(A,F9.6)') 'CCFL Max: ', ccfl_max
16451658
if (viscous) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max
16461659
if (viscous) write (3, '(A,ES16.6)') 'Rc Min: ', Rc_min
16471660

@@ -1683,6 +1696,9 @@ contains
16831696

16841697
if (run_time_info) then
16851698
icfl_max = 0._wp
1699+
if (surface_tension) then
1700+
ccfl_max = 0._wp
1701+
end if
16861702
if (viscous) then
16871703
vcfl_max = 0._wp
16881704
Rc_min = 1.e12_wp

‎src/simulation/m_sim_helpers.fpp‎

Lines changed: 41 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -138,13 +138,13 @@ contains
138138
end subroutine s_compute_enthalpy
139139

140140
!> Computes stability criterion for a specified dt
141-
subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl, vcfl, Rc)
141+
subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl, vcfl, Rc, ccfl)
142142

143143
$:GPU_ROUTINE(parallelism='[seq]')
144144
real(wp), intent(in), dimension(num_vels) :: vel
145145
real(wp), intent(in) :: c, rho
146146
real(wp), intent(inout) :: icfl
147-
real(wp), intent(inout), optional :: vcfl, Rc
147+
real(wp), intent(inout) :: vcfl, Rc, ccfl
148148
real(wp), dimension(2), intent(in) :: Re_l
149149
integer, intent(in) :: j, k, l
150150
real(wp) :: fltr_dtheta
@@ -178,6 +178,24 @@ contains
178178
end if
179179
end if
180180

181+
! Capillary CFL calculation
182+
if (surface_tension) then
183+
if (p > 0) then
184+
#:if not MFC_CASE_OPTIMIZATION or num_dims > 2
185+
if (grid_geometry == 3) then
186+
fltr_dtheta = f_compute_filtered_dtheta(k, l)
187+
ccfl = dt*sqrt(2._wp*pi*sigma/(rho*min(dx(j), dy(k), fltr_dtheta)**3._wp))
188+
else
189+
ccfl = dt*sqrt(2._wp*pi*sigma/(rho*min(dx(j), dy(k), dz(l))**3._wp))
190+
end if
191+
#:endif
192+
else if (n > 0) then
193+
ccfl = dt*sqrt(2._wp*pi*sigma/(rho*min(dx(j), dy(k))**3._wp))
194+
else
195+
ccfl = dt*sqrt(2._wp*pi*sigma/(rho*dx(j)**3._wp))
196+
end if
197+
end if
198+
181199
end subroutine s_compute_stability_from_dt
182200

183201
!> Computes dt for a specified CFL number
@@ -189,14 +207,14 @@ contains
189207
real(wp), intent(inout) :: max_dt
190208
real(wp), dimension(2), intent(in) :: Re_l
191209
integer, intent(in) :: j, k, l
192-
real(wp) :: icfl_dt, vcfl_dt
210+
real(wp) :: vcfl_dt, ccfl_dt
193211
real(wp) :: fltr_dtheta
194212

195213
! Inviscid CFL calculation
196214
if (p > 0 .or. n > 0) then
197-
icfl_dt = cfl_target*f_compute_multidim_cfl_terms(vel, c, j, k, l)
215+
max_dt = cfl_target*f_compute_multidim_cfl_terms(vel, c, j, k, l)
198216
else
199-
icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c))
217+
max_dt = cfl_target*(dx(j)/(abs(vel(1)) + c))
200218
end if
201219

202220
! Viscous calculations
@@ -213,12 +231,26 @@ contains
213231
else
214232
vcfl_dt = cfl_target*(dx(j)**2._wp)/maxval(1/(rho*Re_l))
215233
end if
234+
max_dt = min(max_dt, vcfl_dt)
216235
end if
217236

218-
if (viscous) then
219-
max_dt = min(icfl_dt, vcfl_dt)
220-
else
221-
max_dt = icfl_dt
237+
! Capillary CFL calculations
238+
if (surface_tension) then
239+
if (p > 0) then
240+
#:if not MFC_CASE_OPTIMIZATION or num_dims > 2
241+
if (grid_geometry == 3) then
242+
fltr_dtheta = f_compute_filtered_dtheta(k, l)
243+
ccfl_dt = cfl_target*sqrt(rho*min(dx(j), dy(k), fltr_dtheta)**3._wp/(2._wp*pi*sigma))
244+
else
245+
ccfl_dt = cfl_target*sqrt(rho*min(dx(j), dy(k), dz(l))**3._wp/(2._wp*pi*sigma))
246+
end if
247+
#:endif
248+
else if (n > 0) then
249+
ccfl_dt = cfl_target*sqrt(rho*min(dx(j), dy(k))**3._wp/(2._wp*pi*sigma))
250+
else
251+
ccfl_dt = cfl_target*sqrt(rho*dx(j)**3._wp/(2._wp*pi*sigma))
252+
end if
253+
max_dt = min(max_dt, ccfl_dt)
222254
end if
223255

224256
end subroutine s_compute_dt_from_cfl

0 commit comments

Comments
 (0)