Line data    Source code

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