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..2d1565fe4 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,9 @@ ERF::Evolve () int iteration = 1; timeStep(0, cur_time, iteration); - cur_time += dt[0]; + 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 9c3e1b06c..5cf40753a 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); @@ -40,8 +40,8 @@ ERF::ComputeDt (int step) // so we must add start_time to t_new // 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;