Bugfix for adaptive timestepping (#1666)

Co-authored-by: Brandon Blakeley <[email protected]> Co-authored-by: Ann Almgren <[email protected]>
erf-model · Jul 9, 2024 · fdf6c42 · fdf6c42
1 parent ded446b
commit fdf6c42
Show file tree

Hide file tree

Showing 3 changed files with 10 additions and 6 deletions.
diff --git a/Source/ERF.H b/Source/ERF.H
@@ -494,7 +494,7 @@ private:
     void setSpongeRefFromSounding (bool restarting);
 
     // a wrapper for estTimeStep()
-    void ComputeDt ();
+    void ComputeDt (int step = -1);
 
     // get plotfile name
     [[nodiscard]] std::string PlotFileName (int lev) const;

diff --git a/Source/ERF.cpp b/Source/ERF.cpp
@@ -340,7 +340,7 @@ ERF::Evolve ()
     {
         Print() << "\nCoarse STEP " << step+1 << " starts ..." << std::endl;
 
-        ComputeDt();
+        ComputeDt(step);
 
         // Make sure we have read enough of the boundary plane data to make it through this timestep
         if (input_bndry_planes)
@@ -1972,7 +1972,7 @@ ERF::Evolve_MB (int MBstep, int max_block_step)
 
         Print() << "\nCoarse STEP " << step+1 << " starts ..." << std::endl;
 
-        ComputeDt();
+        ComputeDt(step);
 
         // Make sure we have read enough of the boundary plane data to make it through this timestep
         if (input_bndry_planes)
@@ -2051,8 +2051,8 @@ ERF::writeNow(const Real cur_time, const Real dt_lev, const int nstep, const int
 
         const Real eps = std::numeric_limits<Real>::epsilon() * Real(10.0) * std::abs(cur_time);
 
-        int num_per_old = static_cast<int>(std::round((cur_time-eps-dt_lev) / plot_per));
-        int num_per_new = static_cast<int>(std::round((cur_time-eps       ) / plot_per));
+        int num_per_old = static_cast<int>(std::floor((cur_time-eps-dt_lev) / plot_per));
+        int num_per_new = static_cast<int>(std::floor((cur_time-eps       ) / plot_per));
 
         // Before using these, however, we must test for the case where we're
         // within machine epsilon of the next interval. In that case, increment

diff --git a/Source/TimeIntegration/ERF_ComputeTimestep.cpp b/Source/TimeIntegration/ERF_ComputeTimestep.cpp
@@ -8,7 +8,7 @@ using namespace amrex;
  *
  */
 void
-ERF::ComputeDt ()
+ERF::ComputeDt (int step)
 {
     Vector<Real> dt_tmp(finest_level+1);
 
@@ -25,6 +25,10 @@ ERF::ComputeDt ()
         dt_tmp[lev] = amrex::min(dt_tmp[lev], change_max*dt[lev]);
         n_factor *= nsubsteps[lev];
         dt_0 = amrex::min(dt_0, n_factor*dt_tmp[lev]);
+        if (step == 0){
+            dt_0 *= init_shrink;
+            Print() << "Timestep 0: shrink initial dt by " << init_shrink << std::endl;
+        }
     }
 
     // Limit dt's by the value of stop_time.