Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <algorithm>
       7             : #include <cmath>
       8             : #include <limits>
       9             : #include <pup.h>
      10             : #include <string>
      11             : #include <type_traits>
      12             : #include <typeindex>
      13             : #include <typeinfo>
      14             : #include <unordered_map>
      15             : 
      16             : #include "Options/String.hpp"
      17             : #include "Time/RequestsStepperErrorTolerances.hpp"
      18             : #include "Time/StepChoosers/StepChooser.hpp"
      19             : #include "Time/StepperErrorTolerances.hpp"
      20             : #include "Time/Tags/StepperErrors.hpp"
      21             : #include "Time/TimeStepRequest.hpp"
      22             : #include "Utilities/Serialization/CharmPupable.hpp"
      23             : #include "Utilities/TMPL.hpp"
      24             : 
      25             : /// \cond
      26             : struct NoSuchType;
      27             : /// \endcond
      28             : 
      29             : namespace StepChoosers {
      30             : /*!
      31             :  * \brief Sets a goal based on time-stepper truncation error.
      32             :  *
      33             :  * \details The suggested step is calculated via a simple specialization of the
      34             :  * scheme suggested in \cite Hairer1993. We first compute the aggregated error
      35             :  * measure from the stepper error:
      36             :  *
      37             :  * \f[
      38             :  * E = \max_i(|E_i| / sc_i),
      39             :  * \f]
      40             :  *
      41             :  * where \f$E_i\f$ is the ODE error reported for each individual grid point,
      42             :  * reported by the time stepper, and \f$sc_i\f$ is the step control measure
      43             :  * determined by the tolerances:
      44             :  *
      45             :  * \f[
      46             :  * sc_i = Atol_i + \max(|y_i|,|y_i + E_i|) Rtol_i,
      47             :  * \f]
      48             :  *
      49             :  * and \f$y_i\f$ is the value of the function at the previous step at
      50             :  * grid point \f$i\f$.  (The estimate is more commonly done comparing
      51             :  * with the current step, but using the previous step avoids a memory
      52             :  * allocation in the TimeStepper and should not have a major effect on
      53             :  * the result.)
      54             :  *
      55             :  * When choosing a step size for LTS or when no record of previous
      56             :  * error is available, the step has size:
      57             :  *
      58             :  * \f[
      59             :  * h_{\text{new}} = h \cdot \min\left(F_{\text{max}},
      60             :  * \max\left(F_{\text{min}},
      61             :  * \frac{F_{\text{safety}}}{E^{1/(q + 1)}}\right)\right),
      62             :  * \f]
      63             :  *
      64             :  * where \f$h_{\text{new}}\f$ is the new suggested step size \f$h\f$ is the
      65             :  * previous step size, \f$F_{\text{max}}\f$ is the maximum factor by which we
      66             :  * allow the step to increase, \f$F_{\text{min}}\f$ is the minimum factor by
      67             :  * which we allow the step to decrease. \f$F_{\text{safety}}\f$ is the safety
      68             :  * factor on the computed error -- this forces the step size slightly lower
      69             :  * than we would naively compute so that the result of the step will likely be
      70             :  * within the target error. \f$q\f$ is the order of the stepper error
      71             :  * calculation. Intuitively, we should change the step less drastically for a
      72             :  * higher order stepper.
      73             :  *
      74             :  * When controlling slab size, after the first error calculation, the
      75             :  * error \f$E\f$ is recorded in the \ref DataBoxGroup "DataBox", and
      76             :  * subsequent error calculations use a simple PI scheme suggested in
      77             :  * \cite NumericalRecipes section 17.2.1:
      78             :  *
      79             :  * \f[
      80             :  * h_{\text{new}} = h \cdot \min\left(F_{\text{max}},
      81             :  * \max\left(F_{\text{min}},
      82             :  * F_{\text{safety}} E^{-0.7 / (q + 1)}
      83             :  * E_{\text{prev}}^{0.4 / (q + 1)}\right)\right),
      84             :  * \f]
      85             :  *
      86             :  * where \f$E_{\text{prev}}\f$ is the error computed in the previous
      87             :  * step.  This method is never used for choosing an LTS step because
      88             :  * the restriction of step size changes to factors of two was found to
      89             :  * interfere with the more gradual increase chosen by the PI
      90             :  * controller.
      91             :  *
      92             :  * \note The template parameter `ErrorControlSelector` is used to disambiguate
      93             :  * in the input-file options between `ErrorControl` step choosers that are
      94             :  * based on different variables. This is needed if multiple systems are evolved
      95             :  * in the same executable. The name used for the input file includes
      96             :  * `ErrorControlSelector::name()` if it is provided.
      97             :  */
      98             : template <typename StepChooserUse, typename EvolvedVariableTag,
      99             :           typename ErrorControlSelector = NoSuchType>
     100           1 : class ErrorControl : public StepChooser<StepChooserUse>,
     101             :                      public RequestsStepperErrorTolerances {
     102             :  public:
     103             :   /// \cond
     104             :   ErrorControl() = default;
     105             :   explicit ErrorControl(CkMigrateMessage* /*unused*/) {}
     106             :   using PUP::able::register_constructor;
     107             :   WRAPPED_PUPable_decl_template(ErrorControl);  // NOLINT
     108             :   /// \endcond
     109             : 
     110           0 :   static std::string name() {
     111             :     if constexpr (std::is_same_v<ErrorControlSelector, NoSuchType>) {
     112             :       return "ErrorControl";
     113             :     } else {
     114             :       return "ErrorControl(" + ErrorControlSelector::name() + ")";
     115             :     }
     116             :   }
     117             : 
     118           0 :   struct AbsoluteTolerance {
     119           0 :     using type = double;
     120           0 :     static constexpr Options::String help{"Target absolute tolerance"};
     121           0 :     static type lower_bound() { return 0.0; }
     122             :   };
     123             : 
     124           0 :   struct RelativeTolerance {
     125           0 :     using type = double;
     126           0 :     static constexpr Options::String help{"Target relative tolerance"};
     127           0 :     static type lower_bound() { return 0.0; }
     128             :   };
     129             : 
     130           0 :   struct MaxFactor {
     131           0 :     using type = double;
     132           0 :     static constexpr Options::String help{
     133             :         "Maximum factor to increase the step by"};
     134           0 :     static type lower_bound() { return 1.0; }
     135             :   };
     136             : 
     137           0 :   struct MinFactor {
     138           0 :     using type = double;
     139           0 :     static constexpr Options::String help{
     140             :         "Minimum factor to increase the step by"};
     141           0 :     static type lower_bound() { return 0.0; }
     142           0 :     static type upper_bound() { return 1.0; }
     143             :   };
     144             : 
     145           0 :   struct SafetyFactor {
     146           0 :     using type = double;
     147           0 :     static constexpr Options::String help{
     148             :         "Extra factor to apply to step estimate; can be used to decrease step "
     149             :         "size to improve step acceptance rate."};
     150           0 :     static type lower_bound() { return 0.0; }
     151             :   };
     152             : 
     153           0 :   static constexpr Options::String help{
     154             :       "Sets a goal based on time-stepper truncation error."};
     155           0 :   using options = tmpl::list<AbsoluteTolerance, RelativeTolerance, MaxFactor,
     156             :                              MinFactor, SafetyFactor>;
     157             : 
     158           0 :   ErrorControl(const double absolute_tolerance, const double relative_tolerance,
     159             :                const double max_factor, const double min_factor,
     160             :                const double safety_factor)
     161             :       : absolute_tolerance_{absolute_tolerance},
     162             :         relative_tolerance_{relative_tolerance},
     163             :         max_factor_{max_factor},
     164             :         min_factor_{min_factor},
     165             :         safety_factor_{safety_factor} {}
     166             : 
     167           0 :   using argument_tags = tmpl::list<::Tags::StepperErrors<EvolvedVariableTag>>;
     168             : 
     169           0 :   TimeStepRequest operator()(
     170             :       const typename ::Tags::StepperErrors<EvolvedVariableTag>::type& errors,
     171             :       const double /*previous_step*/) const {
     172             :     // Do not request that the step size be changed if there isn't a new error
     173             :     // estimate
     174             :     if (not errors[1].has_value()) {
     175             :       return {};
     176             :     }
     177             :     double new_step;
     178             :     if (std::is_same_v<StepChooserUse, ::StepChooserUse::LtsStep> or
     179             :         not errors[0].has_value() or errors[0]->order != errors[1]->order) {
     180             :       new_step =
     181             :           errors[1]->step_size.value() *
     182             :           std::clamp(safety_factor_ *
     183             :                          pow(1.0 / std::max(errors[1]->step_error(), 1e-14),
     184             :                              1.0 / (errors[1]->order + 1)),
     185             :                      min_factor_, max_factor_);
     186             :     } else {
     187             :       // From simple advice from Numerical Recipes 17.2.1 regarding a heuristic
     188             :       // for PI step control.
     189             :       const double alpha_factor = 0.7 / (errors[1]->order + 1);
     190             :       const double beta_factor = 0.4 / (errors[0]->order + 1);
     191             :       new_step =
     192             :           errors[1]->step_size.value() *
     193             :           std::clamp(
     194             :               safety_factor_ *
     195             :                   pow(1.0 / std::max(errors[1]->step_error(), 1e-14),
     196             :                       alpha_factor) *
     197             :                   pow(std::max(errors[0]->step_error(), 1e-14), beta_factor),
     198             :               min_factor_, max_factor_);
     199             :     }
     200             :     return ::TimeStepRequest{.size_goal = new_step};
     201             :   }
     202             : 
     203           1 :   bool uses_local_data() const override { return true; }
     204           1 :   bool can_be_delayed() const override { return true; }
     205             : 
     206           1 :   std::unordered_map<std::type_index, StepperErrorTolerances> tolerances()
     207             :       const override {
     208             :     return {{typeid(EvolvedVariableTag),
     209             :              {.estimates = StepperErrorTolerances::Estimates::StepperOrder,
     210             :               .absolute = absolute_tolerance_,
     211             :               .relative = relative_tolerance_}}};
     212             :   }
     213             : 
     214           0 :   void pup(PUP::er& p) override {  // NOLINT
     215             :     StepChooser<StepChooserUse>::pup(p);
     216             :     p | absolute_tolerance_;
     217             :     p | relative_tolerance_;
     218             :     p | min_factor_;
     219             :     p | max_factor_;
     220             :     p | safety_factor_;
     221             :   }
     222             : 
     223             :  private:
     224           0 :   double absolute_tolerance_ = std::numeric_limits<double>::signaling_NaN();
     225           0 :   double relative_tolerance_ = std::numeric_limits<double>::signaling_NaN();
     226           0 :   double max_factor_ = std::numeric_limits<double>::signaling_NaN();
     227           0 :   double min_factor_ = std::numeric_limits<double>::signaling_NaN();
     228           0 :   double safety_factor_ = std::numeric_limits<double>::signaling_NaN();
     229             : };
     230             : /// \cond
     231             : template <typename StepChooserUse, typename EvolvedVariableTag,
     232             :           typename ErrorControlSelector>
     233             : PUP::able::PUP_ID ErrorControl<StepChooserUse, EvolvedVariableTag,
     234             :                                ErrorControlSelector>::my_PUP_ID = 0;  // NOLINT
     235             : /// \endcond
     236             : }  // namespace StepChoosers