From d282849e49992f38351b0cefe10109615377f36b Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 7 Oct 2026 01:39:58 -0400 Subject: [PATCH] Fixed-dt runs: keep the clock on t_step*dt so restarts continue bitwise With a fixed dt, a restart starts its clock at t_step*dt (p_main) and evaluates prescribed IB kinematics there (s_ibm_setup), but an uninterrupted run accumulates mytime = mytime + dt, which drifts (1290 ulps by step 27000 at dt = 6.17e-4). Stage times were likewise mytime + dt. And the last step of every run trimmed dt to finaltime - mytime, so a chunk's final step used a dt off by that drift. A restarted run therefore never continued the uninterrupted one exactly; near-tie immersed-boundary cells amplify such round-off to visible local differences. For fixed dt, advance mytime as (t_step + 1)*dt, form the IB stage times from t_step the same way, and drop the fixed-dt trim (the step count already lands on t_step_stop). cfl_dt runs are unchanged. Co-Authored-By: Claude --- src/simulation/m_global_parameters.fpp | 3 +-- src/simulation/m_start_up.fpp | 17 ++++++++++------- src/simulation/m_time_steppers.fpp | 16 +++++++++++----- src/simulation/p_main.fpp | 1 - 4 files changed, 22 insertions(+), 15 deletions(-) diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 0637b76923..1a3aecde83 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -286,8 +286,7 @@ module m_global_parameters !> @{ !> @} - real(wp) :: mytime !< Current simulation time - real(wp) :: finaltime !< Final simulation time + real(wp) :: mytime !< Current simulation time type(pres_field), allocatable, dimension(:) :: pb_ts type(pres_field), allocatable, dimension(:) :: mv_ts diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index df311bca2b..2a7181595d 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -601,16 +601,13 @@ contains end if end if + ! Land exactly on t_stop. A fixed dt already lands on t_step_stop by step count; trimming it to t_step_stop*dt - mytime + ! would only change the run's last dt by round-off, which a run continuing past that step does not see if (cfl_dt) then if ((mytime + dt) >= t_stop) then dt = t_stop - mytime $:GPU_UPDATE(device='[dt]') end if - else - if ((mytime + dt) >= finaltime) then - dt = finaltime - mytime - $:GPU_UPDATE(device='[dt]') - end if end if if (cfl_dt) then @@ -652,8 +649,14 @@ contains call s_tvd_rk(t_step, time_avg, time_stepper) end if - ! Advance time after RK so source terms see current-step time - mytime = mytime + dt + ! Advance time after RK so source terms see current-step time. With a fixed dt, use the same t_step*dt a restart + ! starts from (p_main): a running sum drifts from it (1290 ulps by step 27000), so a restarted run would see the + ! prescribed IB kinematics, inflow ramps and forcing at slightly different times than the run it continues + if (cfl_dt) then + mytime = mytime + dt + else + mytime = (t_step + 1)*dt + end if if (relax) call s_infinite_relaxation_k(q_cons_ts(1)%vf) diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 5fddc3f132..025947bfc1 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -564,7 +564,7 @@ contains if (ib) then ! check if any IBMS are moving, and if so, update the markers, ghost points, levelsets, and levelset norms if (moving_immersed_boundary_flag) then - call s_propagate_immersed_boundaries(s) + call s_propagate_immersed_boundaries(s, t_step) end if ! update the ghost fluid properties point values based on IB state @@ -831,9 +831,9 @@ contains end subroutine s_apply_synthetic_turbulence_force !> Update immersed boundary positions and velocities at the current Runge-Kutta stage - subroutine s_propagate_immersed_boundaries(s) + subroutine s_propagate_immersed_boundaries(s, t_step) - integer, intent(in) :: s + integer, intent(in) :: s, t_step integer :: i integer :: gbl_id ! used for analytic ib patch motion real(wp) :: t_stage ! time of the state produced by RK stage s (used by prescribed kinematics) @@ -842,8 +842,14 @@ contains if (moving_immersed_boundary_flag) call s_compute_ib_forces(q_prim_vf, fluid_pp) - t_stage = mytime + dt - if (time_stepper == time_stepper_rk3 .and. s == 2) t_stage = mytime + 0.5_wp*dt + if (cfl_dt) then + t_stage = mytime + dt + if (time_stepper == time_stepper_rk3 .and. s == 2) t_stage = mytime + 0.5_wp*dt + else + ! The same t_step*dt form a restart evaluates the kinematics at, so a restart sees bitwise the same body + t_stage = (t_step + 1)*dt + if (time_stepper == time_stepper_rk3 .and. s == 2) t_stage = (t_step + 0.5_wp)*dt + end if $:GPU_PARALLEL_LOOP(private='[i, gbl_id]', copyin='[s, t_stage]') do i = 1, num_ibs diff --git a/src/simulation/p_main.fpp b/src/simulation/p_main.fpp index 0111543fee..2792a69911 100644 --- a/src/simulation/p_main.fpp +++ b/src/simulation/p_main.fpp @@ -54,7 +54,6 @@ program p_main else mytime = t_step*dt end if - finaltime = t_step_stop*dt end if call nvtxEndRange ! INIT