Upper bound to time step when adaptive (#845)

* upper bound to time step when adaptive

* Apply suggestions from code review

docs/source/run/parameters.rst

@@ -65,6 +65,10 @@ General parameters
 * ``hipace.dt`` (`float` or `string`) optional (default `0.`)
     Time step to advance the particle beam. For adaptive time step, use ``"adaptive"``.
 
+* ``hipace.dt_max`` (`float`) optional (default `inf`)
+    Only used if ``hipace.dt = adaptive``. Upper bound of the adaptive time step: if the computed adaptive time step is is larger than ``dt_max``, then ``dt_max`` is used instead.
+    Useful when the plasma profile starts with a very low density (e.g. in the presence of a realistic density ramp), to avoid unreasonably large time steps.
+
 * ``hipace.nt_per_betatron`` (`Real`) optional (default `40.`)
     Only used when using adaptive time step (see ``hipace.dt`` above).
     Number of time steps per betatron period (of the full blowout regime).

src/utils/AdaptiveTimeStep.H

@@ -22,6 +22,7 @@ private:
 
     /** Number of time steps per betatron period for the adaptive time step */
     amrex::Real m_nt_per_betatron = 40.;
+    amrex::Real m_dt_max = std::numeric_limits<amrex::Real>::infinity();
 
 public:
     /** Whether to use an adaptive time step */

src/utils/AdaptiveTimeStep.cpp

@@ -30,6 +30,7 @@ AdaptiveTimeStep::AdaptiveTimeStep (const int nbeams)
     if (str_dt == "adaptive"){
         m_do_adaptive_time_step = true;
         queryWithParser(ppa, "nt_per_betatron", m_nt_per_betatron);
+        queryWithParser(ppa, "dt_max", m_dt_max);
     }
     DeprecatedInput("hipace", "do_adaptive_time_step", "dt = adaptive");
 
@@ -141,7 +142,7 @@ AdaptiveTimeStep::Calculate (amrex::Real& dt, MultiBeam& beams, amrex::Real plas
             std::min(m_timestep_data[ibeam][WhichDouble::MinUz], amrex::get<3>(res));
     }
 
-    // only the last box or at initialiyation the adaptive time step is calculated
+    // only the last box or at initialization the adaptive time step is calculated
     // from the full beam information
     if (it == 0 || initial)
     {
@@ -184,6 +185,7 @@ AdaptiveTimeStep::Calculate (amrex::Real& dt, MultiBeam& beams, amrex::Real plas
         }
         /* set the new time step */
         dt = *std::min_element(new_dts.begin(), new_dts.end());
-
+        // Make sure the new time step is smaller than the upper bound
+        dt = std::min(dt, m_dt_max);
     }
 }