From d27767ba685b781dee115e0576dc465a48a79ea9 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Mon, 24 Aug 2026 11:07:54 -0400 Subject: [PATCH 1/3] update terminal output and add time-step limiting factor to adaptive cases --- src/common/m_mpi_common.fpp | 16 ++++++++++++++++ src/simulation/m_sim_helpers.fpp | 26 +++++++++++++++++--------- src/simulation/m_start_up.fpp | 12 ++++++++---- src/simulation/m_time_steppers.fpp | 22 +++++++++++++++------- 4 files changed, 56 insertions(+), 20 deletions(-) diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 061d8f1d6a..5530539b5a 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -441,6 +441,22 @@ contains end subroutine s_mpi_allreduce_min + !> Reduce a local real vector to its elementwise global minimum across all MPI ranks. + impure subroutine s_mpi_allreduce_min_vec(var_loc, var_glb) + + real(wp), dimension(:), intent(in) :: var_loc + real(wp), dimension(:), intent(out) :: var_glb + +#ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors + + call MPI_ALLREDUCE(var_loc, var_glb, size(var_loc), mpi_p, MPI_MIN, MPI_COMM_WORLD, ierr) +#else + var_glb = var_loc +#endif + + end subroutine s_mpi_allreduce_min_vec + !> Reduce a local real value to its global maximum across all MPI ranks. impure subroutine s_mpi_allreduce_max(var_loc, var_glb) diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index e50b4eed2b..712f2c86b3 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -14,7 +14,11 @@ module m_sim_helpers implicit none - private; public :: s_compute_enthalpy, s_compute_stability_from_dt, s_compute_dt_from_cfl + private; public :: s_compute_enthalpy, s_compute_stability_from_dt, s_compute_dt_from_cfl, dt_limiter, dt_limiter_names + + !> Stability criterion currently limiting the adaptive time step (ICFL, VCFL, or CCFL) + character(len=4) :: dt_limiter = 'none' + character(len=4), dimension(3), parameter :: dt_limiter_names = (/'ICFL', 'VCFL', 'CCFL'/) contains @@ -178,18 +182,22 @@ contains end subroutine s_compute_stability_from_dt - !> Computes dt for a specified CFL number + !> Computes the candidate dts for a specified CFL number: max_dt(1) from the inviscid, max_dt(2) the viscous, and max_dt(3) the + !! capillary criterion (huge where the criterion is inactive) subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) $:GPU_ROUTINE(parallelism='[seq]') real(wp), dimension(num_vels), intent(in) :: vel real(wp), intent(in) :: c, rho - real(wp), intent(inout) :: max_dt + real(wp), dimension(3), intent(inout) :: max_dt real(wp), dimension(2), intent(in) :: Re_l integer, intent(in) :: j, k, l real(wp) :: vcfl_dt, ccfl_dt real(wp) :: fltr_dtheta + max_dt(2) = huge(1._wp) + max_dt(3) = huge(1._wp) + ! Inviscid CFL calculation ! The multi-dimensional CFL terms are written out here rather than ! obtained from a shared helper procedure: NVHPC 25.5's fort2 segfaults @@ -199,15 +207,15 @@ contains #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 if (grid_geometry == 3) then fltr_dtheta = f_compute_filtered_dtheta(k, l) - max_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), dy(k)/(abs(vel(2)) + c), fltr_dtheta/(abs(vel(3)) + c)) + max_dt(1) = cfl_target*min(dx(j)/(abs(vel(1)) + c), dy(k)/(abs(vel(2)) + c), fltr_dtheta/(abs(vel(3)) + c)) else - max_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), dy(k)/(abs(vel(2)) + c), dz(l)/(abs(vel(3)) + c)) + max_dt(1) = cfl_target*min(dx(j)/(abs(vel(1)) + c), dy(k)/(abs(vel(2)) + c), dz(l)/(abs(vel(3)) + c)) end if #:endif else if (n > 0) then - max_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), dy(k)/(abs(vel(2)) + c)) + max_dt(1) = cfl_target*min(dx(j)/(abs(vel(1)) + c), dy(k)/(abs(vel(2)) + c)) else - max_dt = cfl_target*(dx(j)/(abs(vel(1)) + c)) + max_dt(1) = cfl_target*(dx(j)/(abs(vel(1)) + c)) end if ! Viscous calculations @@ -224,7 +232,7 @@ contains else vcfl_dt = cfl_target*(dx(j)**2._wp)/maxval(1/(rho*Re_l)) end if - max_dt = min(max_dt, vcfl_dt) + max_dt(2) = vcfl_dt end if ! Capillary CFL calculations @@ -243,7 +251,7 @@ contains else ccfl_dt = cfl_target*sqrt(rho*dx(j)**3._wp/(2._wp*pi*sigma)) end if - max_dt = min(max_dt, ccfl_dt) + max_dt(3) = ccfl_dt end if end subroutine s_compute_dt_from_cfl diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 8053beb9ce..411a1114b9 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -569,6 +569,7 @@ contains real(wp), intent(inout) :: time_avg integer :: i, eta_hh, eta_mm, eta_ss real(wp) :: eta_sec + character(len=8) :: lim_str !< Time-step limiter tag, e.g. ' (ICFL)' if (cfl_dt) then if (cfl_const_dt .and. t_step == 0) call s_compute_dt() @@ -601,8 +602,11 @@ contains eta_hh = int(eta_sec)/3600 eta_mm = mod(int(eta_sec), 3600)/60 eta_ss = mod(int(eta_sec), 60) - print '(" [", I3, "%] Time ", ES16.6, " dt = ", ES16.6, " @ Time Step = ", I8, " Time Avg = ", ES16.6, " Time/step = ", ES12.6, " ETA (HH:MM:SS) = ", I0, ":", I2.2, ":", I2.2)', & - & int(ceiling(100._wp*(mytime/t_stop))), mytime, dt, t_step, wall_time_avg, wall_time, eta_hh, eta_mm, eta_ss + lim_str = '' + if (cfl_adap_dt) lim_str = ' (' // dt_limiter // ')' + print '(" [", I3, "%] t = ", ES11.4, " dt = ", ES11.4, A, " @ step ", I0, " t/step ", ES9.2, "s (avg ", ES9.2, "s) ETA ", I0, ":", I2.2, ":", I2.2)', & + & int(ceiling(100._wp*(mytime/t_stop))), mytime, dt, trim(lim_str), t_step, wall_time, wall_time_avg, eta_hh, & + & eta_mm, eta_ss end if else if (proc_rank == 0 .and. mod(t_step - t_step_start, t_step_print) == 0) then @@ -610,9 +614,9 @@ contains eta_hh = int(eta_sec)/3600 eta_mm = mod(int(eta_sec), 3600)/60 eta_ss = mod(int(eta_sec), 60) - print '(" [", I3, "%] Time step ", I8, " of ", I0, " @ t_step = ", I8, " Time Avg = ", ES12.6, " Time/step= ", ES12.6, " ETA (HH:MM:SS) = ", I0, ":", I2.2, ":", I2.2)', & + print '(" [", I3, "%] step ", I0, " of ", I0, " (t_step ", I0, ") t/step ", ES9.2, "s (avg ", ES9.2, "s) ETA ", I0, ":", I2.2, ":", I2.2)', & & int(ceiling(100._wp*(real(t_step - t_step_start)/(t_step_stop - t_step_start + 1)))), & - & t_step - t_step_start + 1, t_step_stop - t_step_start + 1, t_step, wall_time_avg, wall_time, eta_hh, & + & t_step - t_step_start + 1, t_step_stop - t_step_start + 1, t_step, wall_time, wall_time_avg, eta_hh, & & eta_mm, eta_ss end if end if diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index daf1c2a734..4cd25c5a3d 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -655,8 +655,9 @@ contains real(wp) :: c !< Cell-avg. sound speed real(wp) :: H !< Cell-avg. enthalpy real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers - real(wp) :: max_dt - real(wp) :: dt_local + real(wp), dimension(3) :: max_dt !< Cell dt candidates (inviscid, viscous, capillary) + real(wp) :: icfl_dt_local, vcfl_dt_local, ccfl_dt_local + real(wp), dimension(3) :: dt_cfl_glb integer :: j, k, l !< Generic loop iterators integer :: fl !< Fluid loop iterator @@ -664,9 +665,11 @@ contains call s_convert_conservative_to_primitive_variables(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, idwint) end if - dt_local = huge(1.0_wp) + icfl_dt_local = huge(1.0_wp) + vcfl_dt_local = huge(1.0_wp) + ccfl_dt_local = huge(1.0_wp) $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, fl, max_dt]', & - & reduction='[[dt_local]]', reductionOp='[min]') + & reduction='[[icfl_dt_local, vcfl_dt_local, ccfl_dt_local]]', reductionOp='[min]') do l = 0, p do k = 0, n do j = 0, m @@ -693,18 +696,23 @@ contains call s_compute_dt_from_cfl(vel, c, max_dt, rho, Re, j, k, l) - dt_local = min(dt_local, max_dt) + icfl_dt_local = min(icfl_dt_local, max_dt(1)) + vcfl_dt_local = min(vcfl_dt_local, max_dt(2)) + ccfl_dt_local = min(ccfl_dt_local, max_dt(3)) end do end do end do $:END_GPU_PARALLEL_LOOP() if (num_procs == 1) then - dt = dt_local + dt_cfl_glb = (/icfl_dt_local, vcfl_dt_local, ccfl_dt_local/) else - call s_mpi_allreduce_min(dt_local, dt) + call s_mpi_allreduce_min_vec((/icfl_dt_local, vcfl_dt_local, ccfl_dt_local/), dt_cfl_glb) end if + dt = minval(dt_cfl_glb) + dt_limiter = dt_limiter_names(minloc(dt_cfl_glb, dim=1)) + $:GPU_UPDATE(device='[dt]') end subroutine s_compute_dt From f3b35a5c958f5ca3a47f7d36aa4a977b6e5155a6 Mon Sep 17 00:00:00 2001 From: Ben Wilfong Date: Tue, 25 Aug 2026 12:26:05 -0400 Subject: [PATCH 2/3] add collision based CFL limit and ramp factor --- docs/documentation/case.md | 6 +++ src/simulation/m_collisions.fpp | 26 ++++++++++--- src/simulation/m_global_parameters.fpp | 2 + src/simulation/m_sim_helpers.fpp | 4 +- src/simulation/m_start_up.fpp | 10 ++++- src/simulation/m_time_steppers.fpp | 51 +++++++++++++++++--------- toolchain/mfc/params/definitions.py | 17 ++++++++- toolchain/mfc/params/descriptions.py | 2 + 8 files changed, 91 insertions(+), 27 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 83f5a50fd5..42b47cc194 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -362,6 +362,7 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi | `coefficient_of_restitution` | Real | A number 0 to 1 describing how elastic IB collisions are | | `collision_model` | Integer | Integer to select the collision model being used for IB collisions. | | `collision_time` | Real | Amount of simulation time used to resolve collisions | +| `collision_temporal_resolution` | Integer | Minimum number of adaptive time steps used to resolve each collision | | `ib_coefficient_of_friction` | Real | Coefficient of friction used in IB collisions | These parameters should be prepended with `patch_ib(j)%` where $j$ is the patch index. @@ -412,6 +413,8 @@ Additional details on this specification can be found in [NACA airfoil](https:// - `collision_time` is approximately the amount of simulation time used to resolve collisions. This is handled by modifying the spring constant used to apply collision forces. +- `collision_temporal_resolution` restricts the adaptive time step (`cfl_adap_dt`) to at most `collision_time / collision_temporal_resolution` while any collision is occurring, so that each collision is resolved with at least that many time steps. Pairing it with `ramp_ratio` limits how quickly the time step grows back once the collision ends. + - `ib_coefficient_of_friction` is the coefficient of friction used in IB collisions. - `ib_neighborhood_radius` controls the size of the neighborhood size. A value of $r$ indicates that any given rank is aware of IBs up to $r$ ranks away. This value defaults to 0, which leaves the radius unset so that it is selected automatically. This parameter is required to strong-scale a case when IBs eventually grow to be larger than one full processor domain wide. @@ -536,6 +539,7 @@ See @ref equations "Equations" for the mathematical models these parameters cont | `cfl_const_dt` | Logical | CFL based non-adaptive time-stepping | | `cfl_dt` | Logical | Enable CFL-based time stepping | | `cfl_target` | Real | Specified CFL value | +| `ramp_ratio` | Real | Maximum factor by which the adaptive time step may grow per time step | | `n_start` | Integer | Save file from which to start simulation | | `t_save` | Real | Time duration between data output | | `t_stop` | Real | Simulation stop time | @@ -698,6 +702,8 @@ To restart the simulation from $k$-th time step, set `t_step_start = k`; see @re - `cfl_target` specifies the target CFL value +- `ramp_ratio` limits how much the adaptive time step can grow from one time step to the next: `dt` is capped at `ramp_ratio` times the previous `dt`. Must be at least 1. When unset, the time step growth is unlimited. + - `n_start` specifies the save file to start at - `t_save` specifies the time interval between data output during the simulation diff --git a/src/simulation/m_collisions.fpp b/src/simulation/m_collisions.fpp index fadaa7ece0..329d43488f 100644 --- a/src/simulation/m_collisions.fpp +++ b/src/simulation/m_collisions.fpp @@ -21,7 +21,7 @@ module m_collisions implicit none private; public :: s_apply_collision_forces, s_initialize_collisions_module, s_finalize_collisions_module, & - & f_local_rank_owns_location, f_neighborhood_ranks_own_location, ib_gbl_idx_lookup + & f_local_rank_owns_location, f_neighborhood_ranks_own_location, ib_gbl_idx_lookup, collisions_active ! overlap distances for computing collisions integer, allocatable, dimension(:,:) :: collision_lookup real(wp), allocatable, dimension(:,:) :: wall_overlap_distances @@ -32,6 +32,9 @@ module m_collisions integer, dimension(:), allocatable :: ib_gbl_idx_lookup $:GPU_DECLARE(create='[ib_gbl_idx_lookup]') + !> true when any IB-IB or IB-wall contact was detected on this rank since the last adaptive-dt computation + logical :: collisions_active + contains subroutine s_initialize_collisions_module() @@ -47,6 +50,7 @@ contains @:ALLOCATE(wall_overlap_distances(num_local_ibs_max*27, 6)) wall_overlap_distances = 0 + collisions_active = .false. $:GPU_UPDATE(device='[wall_overlap_distances]') $:GPU_UPDATE(device='[ib_coefficient_of_friction]') @@ -59,15 +63,19 @@ contains type(integer_field), intent(in) :: ib_markers real(wp), dimension(num_ibs, 3), intent(inout) :: forces, torques integer :: num_considered_collisions + logical :: any_wall_collision ! return if no collisions if (collision_model == 0) return ! get is distance used in the force calculation with each IB and each wall - call s_detect_wall_collisions() + call s_detect_wall_collisions(any_wall_collision) call s_detect_ib_collisions(ghost_points, ib_markers, num_gps, num_considered_collisions) + ! accumulate across RK stages; consumed (and reset) by s_compute_dt once per time step + collisions_active = collisions_active .or. any_wall_collision .or. (num_considered_collisions > 0) + select case (collision_model) case (1) ! soft sphere model call s_apply_wall_collision_forces_soft_sphere(forces, torques) @@ -338,14 +346,16 @@ contains end subroutine s_detect_ib_collisions !> @brief uses boundary conditions and particle locations to check for wall conditions - subroutine s_detect_wall_collisions() + subroutine s_detect_wall_collisions(any_wall_collision) - integer :: gp_idx, i, j, k, patch_id - real(wp) :: edge_location, overlap_distance + logical, intent(out) :: any_wall_collision + integer :: gp_idx, i, j, k, patch_id + real(wp) :: edge_location, overlap_distance, max_overlap ! iterate over all ghost points to detect the one that is most-overlapping in each direction - $:GPU_PARALLEL_LOOP(private='[patch_id, edge_location, overlap_distance]') + max_overlap = 0._wp + $:GPU_PARALLEL_LOOP(private='[patch_id, edge_location, overlap_distance]', reduction='[[max_overlap]]', reductionOp='[max]') do patch_id = 1, num_ibs #:for X, DIR, IDX in [('x', 1, 1), ('y', 2, 3), ('z', 3, 5)] ! check if the boundaries are either of the two conditions we should compute collisions with @@ -360,6 +370,7 @@ contains overlap_distance = 0._wp end if wall_overlap_distances(patch_id, ${IDX}$) = overlap_distance + max_overlap = max(max_overlap, overlap_distance) end if if (ib_bc_${X}$%end == BC_SLIP_WALL .or. ib_bc_${X}$%end == BC_NO_SLIP_WALL) then @@ -370,11 +381,14 @@ contains overlap_distance = 0._wp end if wall_overlap_distances(patch_id, ${IDX}$ + 1) = overlap_distance + max_overlap = max(max_overlap, overlap_distance) end if #:endfor end do $:END_GPU_PARALLEL_LOOP() + any_wall_collision = max_overlap > 0._wp + end subroutine s_detect_wall_collisions !> @brief function checks if this local MPI processor owns this specific collision diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 8bfc807d71..fc4e4e30e8 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -358,6 +358,7 @@ contains dt = dflt_real cfl_dt = .false. cfl_target = dflt_real + ramp_ratio = dflt_real t_step_stop = dflt_int t_step_save = dflt_int @@ -488,6 +489,7 @@ contains ! Immersed Boundaries (sim-specific extras) ib_neighborhood_radius = 0 collision_model = 0 + collision_temporal_resolution = 0 coefficient_of_restitution = dflt_real collision_time = dflt_real ib_coefficient_of_friction = dflt_real diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index 712f2c86b3..4210ff50cf 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -16,9 +16,9 @@ module m_sim_helpers private; public :: s_compute_enthalpy, s_compute_stability_from_dt, s_compute_dt_from_cfl, dt_limiter, dt_limiter_names - !> Stability criterion currently limiting the adaptive time step (ICFL, VCFL, or CCFL) + !> Criterion currently limiting the adaptive time step (ICFL, VCFL, CCFL, the collision cap, or the ramp limiter) character(len=4) :: dt_limiter = 'none' - character(len=4), dimension(3), parameter :: dt_limiter_names = (/'ICFL', 'VCFL', 'CCFL'/) + character(len=4), dimension(4), parameter :: dt_limiter_names = (/'ICFL', 'VCFL', 'CCFL', 'COLL'/) contains diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 0bb5f9dee8..73494daff4 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -569,6 +569,7 @@ contains real(wp), intent(inout) :: time_avg integer :: i, eta_hh, eta_mm, eta_ss real(wp) :: eta_sec + real(wp) :: dt_floor character(len=8) :: lim_str !< Time-step limiter tag, e.g. ' (ICFL)' if (cfl_dt) then @@ -578,7 +579,14 @@ contains if (t_step == 0) dt_init = dt - if (dt < 1.e-3_wp*dt_init .and. cfl_adap_dt .and. proc_rank == 0) then + ! the collision restriction deliberately drops dt to collision_time/collision_temporal_resolution, so lower the + ! runaway-dt abort threshold below that cap when it is enabled + dt_floor = 1.e-3_wp*dt_init + if (collision_model > 0 .and. collision_temporal_resolution > 0) then + dt_floor = min(dt_floor, 1.e-3_wp*collision_time/real(collision_temporal_resolution, wp)) + end if + + if (dt < dt_floor .and. cfl_adap_dt .and. proc_rank == 0) then print *, "Delta t = ", dt call s_mpi_abort("Delta t has become too small") end if diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 4cd25c5a3d..0d031b8e28 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -17,6 +17,7 @@ module m_time_steppers use m_bubbles_EE use m_bubbles_EL use m_ibm + use m_collisions, only: collisions_active use m_mpi_proxy use m_boundary_common use m_helper @@ -647,27 +648,30 @@ contains real(wp), dimension(num_vels) :: vel !< Cell-avg. velocity real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction #:endif - real(wp) :: vel_sum !< Cell-avg. velocity sum - real(wp) :: pres !< Cell-avg. pressure - real(wp) :: gamma !< Cell-avg. sp. heat ratio - real(wp) :: pi_inf !< Cell-avg. liquid stiffness function - real(wp) :: qv !< Cell-avg. fluid reference energy - real(wp) :: c !< Cell-avg. sound speed - real(wp) :: H !< Cell-avg. enthalpy - real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers - real(wp), dimension(3) :: max_dt !< Cell dt candidates (inviscid, viscous, capillary) - real(wp) :: icfl_dt_local, vcfl_dt_local, ccfl_dt_local - real(wp), dimension(3) :: dt_cfl_glb - integer :: j, k, l !< Generic loop iterators - integer :: fl !< Fluid loop iterator + real(wp) :: vel_sum !< Cell-avg. velocity sum + real(wp) :: pres !< Cell-avg. pressure + real(wp) :: gamma !< Cell-avg. sp. heat ratio + real(wp) :: pi_inf !< Cell-avg. liquid stiffness function + real(wp) :: qv !< Cell-avg. fluid reference energy + real(wp) :: c !< Cell-avg. sound speed + real(wp) :: H !< Cell-avg. enthalpy + real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers + real(wp), dimension(3) :: max_dt !< Cell dt candidates (inviscid, viscous, capillary) + real(wp) :: icfl_dt_local, vcfl_dt_local, ccfl_dt_local, coll_dt_local + real(wp), dimension(4) :: dt_candidates_glb !< Global dt candidates (ICFL, VCFL, CCFL, collision cap) + real(wp) :: dt_prev + integer :: j, k, l !< Generic loop iterators + integer :: fl !< Fluid loop iterator if (.not. igr) then call s_convert_conservative_to_primitive_variables(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, idwint) end if + dt_prev = dt icfl_dt_local = huge(1.0_wp) vcfl_dt_local = huge(1.0_wp) ccfl_dt_local = huge(1.0_wp) + coll_dt_local = huge(1.0_wp) $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, fl, max_dt]', & & reduction='[[icfl_dt_local, vcfl_dt_local, ccfl_dt_local]]', reductionOp='[min]') do l = 0, p @@ -704,14 +708,27 @@ contains end do $:END_GPU_PARALLEL_LOOP() + ! restrict the time step so an ongoing collision spans at least collision_temporal_resolution time steps; the collision + ! flag is rank-local, so the cap enters as a candidate before the global elementwise min propagates it to all ranks + if (collision_model > 0 .and. collision_temporal_resolution > 0) then + if (collisions_active) coll_dt_local = collision_time/real(collision_temporal_resolution, wp) + collisions_active = .false. + end if + if (num_procs == 1) then - dt_cfl_glb = (/icfl_dt_local, vcfl_dt_local, ccfl_dt_local/) + dt_candidates_glb = (/icfl_dt_local, vcfl_dt_local, ccfl_dt_local, coll_dt_local/) else - call s_mpi_allreduce_min_vec((/icfl_dt_local, vcfl_dt_local, ccfl_dt_local/), dt_cfl_glb) + call s_mpi_allreduce_min_vec((/icfl_dt_local, vcfl_dt_local, ccfl_dt_local, coll_dt_local/), dt_candidates_glb) end if - dt = minval(dt_cfl_glb) - dt_limiter = dt_limiter_names(minloc(dt_cfl_glb, dim=1)) + dt = minval(dt_candidates_glb) + dt_limiter = dt_limiter_names(minloc(dt_candidates_glb, dim=1)) + + ! limit how much the time step can grow relative to the previous step + if (ramp_ratio > 0._wp .and. dt_prev > 0._wp .and. ramp_ratio*dt_prev < dt) then + dt = ramp_ratio*dt_prev + dt_limiter = 'RAMP' + end if $:GPU_UPDATE(device='[dt]') diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 5f218f5d05..7bcedaf625 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -365,6 +365,8 @@ def get_value_label(param_name: str, value: int) -> str: "t_step_save": {"min": 1}, "t_step_print": {"min": 1}, "cfl_target": {"min": 0}, + "collision_temporal_resolution": {"min": 1}, + "ramp_ratio": {"min": 1}, # WENO "weno_eps": {"min": 0}, # MUSCL @@ -433,6 +435,16 @@ def get_value_label(param_name: str, value: int) -> str: "requires": ["ib", "coefficient_of_restitution", "collision_time"], } }, + "collision_temporal_resolution": { + "when_set": { + "requires": ["collision_model", "cfl_adap_dt"], + } + }, + "ramp_ratio": { + "when_set": { + "requires": ["cfl_adap_dt"], + } + }, "acoustic_source": { "when_true": { "requires": ["num_source"], @@ -591,7 +603,7 @@ def _load(): _r(n, INT, {"time"}) _r("dt", REAL, {"time"}, math=r"\f$\Delta t\f$") _r("cfl_target", REAL, {"time"}, math=r"\f$\mathrm{CFL}\f$") - for n in ["adap_dt_tol", "t_stop", "t_save"]: + for n in ["adap_dt_tol", "t_stop", "t_save", "ramp_ratio"]: _r(n, REAL, {"time"}) for n in ["cfl_adap_dt", "cfl_const_dt", "cfl_dt", "adap_dt"]: _r(n, LOG, {"time"}) @@ -664,6 +676,7 @@ def _load(): _r("ib_neighborhood_radius", INT, {"ib"}) _r("ib", LOG, {"ib"}) _r("collision_model", INT, {"ib"}) + _r("collision_temporal_resolution", INT, {"ib"}) _r("coefficient_of_restitution", REAL, {"ib"}) _r("collision_time", REAL, {"ib"}) _r("ib_coefficient_of_friction", REAL, {"ib"}) @@ -1386,6 +1399,8 @@ def _decl(targets: set, *names: str) -> None: "turb_pos", "synth_L", "collision_model", + "collision_temporal_resolution", + "ramp_ratio", "coefficient_of_restitution", "collision_time", "ib_coefficient_of_friction", diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index 589d00b553..be48a49426 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -279,7 +279,9 @@ # IB collision parameters "coefficient_of_restitution": "Real number describing the elasticity of collisions from 0 (perfectly ineleastic) to 1 (perfectly elastic)", "collision_model": "Integer selecting the collision model being used. 0 for no collision. 1 for soft-sphere collisions", + "collision_temporal_resolution": "Minimum number of adaptive time steps used to resolve collision_time while collisions are occurring", "collision_time": "Amount of simulation time each collision will take to resolve", + "ramp_ratio": "Maximum factor by which the adaptive time step may grow from one time step to the next", "ib_coefficient_of_friction": "coefficient of friction used in IB collisions", } From aea8cda950d04cf744c261a1b6797a41557efb14 Mon Sep 17 00:00:00 2001 From: Ben Wilfong Date: Tue, 25 Aug 2026 12:36:24 -0400 Subject: [PATCH 3/3] review comments --- src/simulation/m_sim_helpers.fpp | 2 +- src/simulation/m_time_steppers.fpp | 10 ++++++++-- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index 4210ff50cf..ad4be13697 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -189,7 +189,7 @@ contains $:GPU_ROUTINE(parallelism='[seq]') real(wp), dimension(num_vels), intent(in) :: vel real(wp), intent(in) :: c, rho - real(wp), dimension(3), intent(inout) :: max_dt + real(wp), dimension(3), intent(out) :: max_dt real(wp), dimension(2), intent(in) :: Re_l integer, intent(in) :: j, k, l real(wp) :: vcfl_dt, ccfl_dt diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 0d031b8e28..ee4748ec92 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -658,6 +658,7 @@ contains real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers real(wp), dimension(3) :: max_dt !< Cell dt candidates (inviscid, viscous, capillary) real(wp) :: icfl_dt_local, vcfl_dt_local, ccfl_dt_local, coll_dt_local + real(wp), dimension(4) :: dt_candidates_loc !< Rank-local dt candidates (ICFL, VCFL, CCFL, collision cap) real(wp), dimension(4) :: dt_candidates_glb !< Global dt candidates (ICFL, VCFL, CCFL, collision cap) real(wp) :: dt_prev integer :: j, k, l !< Generic loop iterators @@ -715,10 +716,15 @@ contains collisions_active = .false. end if + dt_candidates_loc(1) = icfl_dt_local + dt_candidates_loc(2) = vcfl_dt_local + dt_candidates_loc(3) = ccfl_dt_local + dt_candidates_loc(4) = coll_dt_local + if (num_procs == 1) then - dt_candidates_glb = (/icfl_dt_local, vcfl_dt_local, ccfl_dt_local, coll_dt_local/) + dt_candidates_glb = dt_candidates_loc else - call s_mpi_allreduce_min_vec((/icfl_dt_local, vcfl_dt_local, ccfl_dt_local, coll_dt_local/), dt_candidates_glb) + call s_mpi_allreduce_min_vec(dt_candidates_loc, dt_candidates_glb) end if dt = minval(dt_candidates_glb)