From a3f1611ab28dc9729f14cff464023f35daa6ef05 Mon Sep 17 00:00:00 2001 From: David Wiersema Date: Fri, 8 May 2026 00:03:41 -0700 Subject: [PATCH 1/2] Use double precision for simulation time, even when compiling in single precision. When compiling in single precision, the time step and current time were float32, which could lead to some unintended behavior when the number of time steps is large. For 6 hour Bomex simulations with fixed_dt=0.3, the actual dt is 0.29980469 until time 16384, after which dt=0.30078125. This causes the single precision run to reach 6 hrs in 20 less time steps relative to the double precision run. --- Source/ERF.H | 2 +- Source/ERF.cpp | 8 ++++---- Source/TimeIntegration/ERF_ComputeTimestep.cpp | 10 +++++----- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/Source/ERF.H b/Source/ERF.H index 08bc5eb4f..125e6282f 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -707,7 +707,7 @@ private: void setSpongeRefFromSounding (bool restarting); // a wrapper for estTimeStep() - void ComputeDt (int step = -1); + void ComputeDt (int step = -1, double cur_time_d = 0.0); // get plotfile name [[nodiscard]] std::string PlotFileName (int lev) const; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 6016ad726..2e1f121f6 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -609,8 +609,8 @@ ERF::Evolve () // // cur_time = t_new is elapsed time, not total time // stop_time is total time - // - Real cur_time = t_new[0]; + // Tracked in double to avoid float32 drift over many timesteps in single-precision builds. + double cur_time = static_cast(t_new[0]); // Take one coarse timestep by calling timeStep -- which recursively calls timeStep // for finer levels (with or without subcycling) @@ -622,7 +622,7 @@ ERF::Evolve () } Print() << "\nCoarse STEP " << step+1 << " starts ..." << std::endl; - ComputeDt(step); + ComputeDt(step, cur_time); // Make sure we have read enough of the boundary plane data to make it through this timestep if (input_bndry_planes) @@ -657,7 +657,7 @@ ERF::Evolve () int iteration = 1; timeStep(0, cur_time, iteration); - cur_time += dt[0]; + cur_time += static_cast(dt[0]); Print() << "Coarse STEP " << step+1 << " ends." << " TIME = " << cur_time << " DT = " << dt[0] << std::endl; diff --git a/Source/TimeIntegration/ERF_ComputeTimestep.cpp b/Source/TimeIntegration/ERF_ComputeTimestep.cpp index 9c3e1b06c..fd6342471 100644 --- a/Source/TimeIntegration/ERF_ComputeTimestep.cpp +++ b/Source/TimeIntegration/ERF_ComputeTimestep.cpp @@ -8,7 +8,7 @@ using namespace amrex; * */ void -ERF::ComputeDt (int step) +ERF::ComputeDt (int step, double cur_time_d) { Vector dt_tmp(finest_level+1); @@ -36,12 +36,12 @@ ERF::ComputeDt (int step) } // // Limit dt by the value of stop_time. - // Recall that stop_time is total time, but t_new is elapsed time, - // so we must add start_time to t_new + // Use double-precision cur_time_d (passed from Evolve) to avoid float32 drift + // causing premature dt shortening in single-precision builds. // const Real eps = Real(1.e-3)*dt_0; - if (t_new[0] + dt_0 > (stop_time - start_time) - eps) { - dt_0 = (stop_time - start_time) - t_new[0]; + if (cur_time_d + static_cast(dt_0) > static_cast(stop_time - start_time) - static_cast(eps)) { + dt_0 = static_cast(static_cast(stop_time - start_time) - cur_time_d); } dt[0] = dt_0; From bcfa503f2f7beec8dd25a854e75b20bec5b21e43 Mon Sep 17 00:00:00 2001 From: David Wiersema Date: Fri, 8 May 2026 12:20:50 -0700 Subject: [PATCH 2/2] For single precision, t_new[0] would drift several seconds from the true elapsed time due to representation of the time step as float32. This change assigns t_new[0] using the double cur_time after each time step, which avoids the accumulation of float32 quantization error. The correct timestamp is now provided to physics routines and the plt and chk Header files. There should be no change when compiled in double precision since t_new[0] is already a double, so the static_cast does nothing. --- Source/ERF.cpp | 2 ++ Source/TimeIntegration/ERF_ComputeTimestep.cpp | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 2e1f121f6..2d1565fe4 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -658,6 +658,8 @@ ERF::Evolve () timeStep(0, cur_time, iteration); cur_time += static_cast(dt[0]); + // Sync t_new[0] from accurate double to prevent float32 accumulation drift in SP builds. + t_new[0] = static_cast(cur_time); Print() << "Coarse STEP " << step+1 << " ends." << " TIME = " << cur_time << " DT = " << dt[0] << std::endl; diff --git a/Source/TimeIntegration/ERF_ComputeTimestep.cpp b/Source/TimeIntegration/ERF_ComputeTimestep.cpp index fd6342471..5cf40753a 100644 --- a/Source/TimeIntegration/ERF_ComputeTimestep.cpp +++ b/Source/TimeIntegration/ERF_ComputeTimestep.cpp @@ -36,8 +36,8 @@ ERF::ComputeDt (int step, double cur_time_d) } // // Limit dt by the value of stop_time. - // Use double-precision cur_time_d (passed from Evolve) to avoid float32 drift - // causing premature dt shortening in single-precision builds. + // Recall that stop_time is total time, but t_new is elapsed time, + // so we must add start_time to t_new // const Real eps = Real(1.e-3)*dt_0; if (cur_time_d + static_cast(dt_0) > static_cast(stop_time - start_time) - static_cast(eps)) {