Update terminal output and add limiting factor to adaptive output

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

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)
 

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

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

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
+
+    !> 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(4), parameter :: dt_limiter_names = (/'ICFL', 'VCFL', 'CCFL', 'COLL'/)
 
 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(out)       :: 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

src/simulation/m_start_up.fpp

@@ -569,6 +569,8 @@ 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
             if (cfl_const_dt .and. t_step == 0) call s_compute_dt()
@@ -577,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
@@ -601,18 +610,21 @@ 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
                 eta_sec = wall_time_avg*real(t_step_stop - t_step, wp)
                 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