ErrorControl.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <pup.h>
7 #include <pup_stl.h>
8 #include <utility>
9 
11 #include "Options/Options.hpp"
13 #include "Parallel/PupStlCpp17.hpp"
14 #include "Time/StepChoosers/StepChooser.hpp" // IWYU pragma: keep
15 #include "Time/Tags.hpp"
16 #include "Utilities/TMPL.hpp"
17 #include "Utilities/TypeTraits/IsIterable.hpp"
18 
19 /// \cond
20 class Time;
21 namespace Parallel {
22 template <typename Metavariables>
23 class GlobalCache;
24 } // namespace Parallel
25 /// \endcond
26 
27 namespace StepChoosers {
28 namespace Tags {
29 /// \brief The stepper error measure computed in the previous application of the
30 /// `StepChooser::ErrorControl` step chooser.
31 ///
32 /// \details This tag is templated on `EvolvedVariableTag`, as a separate error
33 /// measure should be stored for each evolved variables and corresponding
34 /// `TimeSteppers::History` in the system, because the stepper is separately
35 /// applied to each `TimeSteppers::History` object.
36 template <typename EvolvedVariableTag>
39 };
40 } // namespace Tags
41 
42 /*!
43  * \brief Suggests a step size based on a target absolute and/or relative error
44  * measure.
45  *
46  * \details The suggested step is calculated via a simple specialization of the
47  * scheme suggested in \cite Hairer1993. We first compute the aggregated error
48  * measure from the stepper error:
49  *
50  * \f[
51  * E = \max_i(|E_i| / sc_i),
52  * \f]
53  *
54  * where \f$E_i\f$ is the ODE error reported for each individual grid point,
55  * reported by the time stepper, and \f$sc_i\f$ is the step control measure
56  * determined by the tolerances:
57  *
58  * \f[
59  * sc_i = Atol_i + \max(|y_i|,|y_i + E_i|) Rtol_i,
60  * \f]
61  *
62  * and \f$y_i\f$ is the value of the function at the current step at grid point
63  * \f$i\f$.
64  *
65  * When no previous record of previous error is available, the step has size:
66  *
67  * \f[
68  * h_{\text{new}} = h \cdot \min\left(F_{\text{max}},
69  * \max\left(F_{\text{min}},
70  * \frac{F_{\text{safety}}}{E^{1/(q + 1)}}\right)\right),
71  * \f]
72  *
73  * where \f$h_{\text{new}}\f$ is the new suggested step size \f$h\f$ is the
74  * previous step size, \f$F_{\text{max}}\f$ is the maximum factor by which we
75  * allow the step to increase, \f$F_{\text{min}}\f$ is the minimum factor by
76  * which we allow the step to decrease. \f$F_{\text{safety}}\f$ is the safety
77  * factor on the computed error -- this forces the step size slightly lower
78  * than we would naively compute so that the result of the step will likely be
79  * within the target error. \f$q\f$ is the order of the stepper error
80  * calculation. Intuitively, we should change the step less drastically for a
81  * higher order stepper.
82  *
83  * After the first error calculation, the error \f$E\f$ is recorded in the \ref
84  * DataBoxGroup "DataBox" using tag
85  * `StepChoosers::Tags::PreviousStepError<EvolvedVariablesTag>`, and subsequent
86  * error calculations use a simple PI scheme suggested in \cite NumericalRecipes
87  * 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  * \frac{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 step.
97  *
98  * \note The template parameter `ErrorControlSelector` is used to disambiguate
99  * in the input-file options between `ErrorControl` step choosers that are
100  * based on different variables. This is needed if multiple systems are evolved
101  * in the same executable. The name used for the input file includes
102  * `ErrorControlSelector::name()` if it is provided.
103  */
104 template <typename EvolvedVariableTag,
105  typename ErrorControlSelector = NoSuchType>
106 class ErrorControl : public StepChooser<StepChooserUse::LtsStep> {
107  public:
108  using evolved_variable_type = typename EvolvedVariableTag::type;
109  using error_variable_type =
111  EvolvedVariableTag>::type;
112 
113  /// \cond
114  ErrorControl() = default;
115  explicit ErrorControl(CkMigrateMessage* /*unused*/) noexcept {}
116  using PUP::able::register_constructor;
118  /// \endcond
119 
120  static std::string name() noexcept {
121  if constexpr (std::is_same_v<ErrorControlSelector, NoSuchType>) {
122  return "ErrorControl";
123  } else {
124  return "ErrorControl(" + ErrorControlSelector::name() + ")";
125  }
126  }
127 
129  using type = double;
130  static constexpr Options::String help{"Target absolute tolerance"};
131  static type lower_bound() noexcept { return 0.0; }
132  };
133 
135  using type = double;
136  static constexpr Options::String help{"Target relative tolerance"};
137  static type lower_bound() noexcept { return 0.0; }
138  };
139 
140  struct MaxFactor {
141  using type = double;
142  static constexpr Options::String help{
143  "Maximum factor to increase the step by"};
144  static type lower_bound() noexcept { return 1.0; }
145  };
146 
147  struct MinFactor {
148  using type = double;
149  static constexpr Options::String help{
150  "Minimum factor to increase the step by"};
151  static type lower_bound() noexcept { return 0.0; }
152  static type upper_bound() noexcept { return 1.0; }
153  };
154 
155  struct SafetyFactor {
156  using type = double;
157  static constexpr Options::String help{
158  "Extra factor to apply to step estimate; can be used to decrease step "
159  "size to improve step acceptance rate."};
160  static type lower_bound() noexcept { return 0.0; }
161  };
162 
163  static constexpr Options::String help{
164  "Chooses a step based on a target relative and absolute error tolerance"};
165  using options = tmpl::list<AbsoluteTolerance, RelativeTolerance, MaxFactor,
166  MinFactor, SafetyFactor>;
167 
168  ErrorControl(const double absolute_tolerance, const double relative_tolerance,
169  const double max_factor, const double min_factor,
170  const double safety_factor) noexcept
171  : absolute_tolerance_{absolute_tolerance},
172  relative_tolerance_{relative_tolerance},
173  max_factor_{max_factor},
174  min_factor_{min_factor},
175  safety_factor_{safety_factor} {}
176 
177  using argument_tags =
178  tmpl::list<::Tags::HistoryEvolvedVariables<EvolvedVariableTag>,
181 
182  using return_tags = tmpl::list<Tags::PreviousStepError<EvolvedVariableTag>>;
183 
184  using simple_tags = tmpl::list<Tags::PreviousStepError<EvolvedVariableTag>>;
185 
186  template <typename Metavariables, typename History, typename TimeStepper>
187  std::pair<double, bool> operator()(
188  const gsl::not_null<std::optional<double>*> previous_step_error,
189  const History& history, const error_variable_type& error,
190  const bool& stepper_error_updated, const TimeStepper& stepper,
191  const double previous_step,
192  const Parallel::GlobalCache<Metavariables>& /*cache*/) const noexcept {
193  // request that the step size not be changed if there isn't a new error
194  // estimate
195  if (not stepper_error_updated) {
196  return std::make_pair(std::numeric_limits<double>::infinity(), true);
197  }
198  const double l_inf_error =
199  error_calc_impl(history.most_recent_value(), error);
200  double new_step;
201  if (not previous_step_error->has_value()) {
202  new_step = previous_step *
203  std::clamp(safety_factor_ *
204  pow(1.0 / std::max(l_inf_error, 1e-14),
205  1.0 / (stepper.error_estimate_order() + 1)),
206  min_factor_, max_factor_);
207  } else {
208  // From simple advice from Numerical Recipes 17.2.1 regarding a heuristic
209  // for PI step control.
210  const double alpha_factor = 0.7 / (stepper.error_estimate_order() + 1);
211  const double beta_factor = 0.4 / (stepper.error_estimate_order() + 1);
212  new_step =
213  previous_step *
214  std::clamp(
215  safety_factor_ *
216  pow(1.0 / std::max(l_inf_error, 1e-14), alpha_factor) *
217  pow(1.0 / std::max(previous_step_error->value(), 1e-14),
218  beta_factor),
219  min_factor_, max_factor_);
220  }
221  *previous_step_error = l_inf_error;
222  return std::make_pair(new_step, l_inf_error <= 1.0);
223  }
224 
225  void pup(PUP::er& p) noexcept override { // NOLINT
226  p | absolute_tolerance_;
227  p | relative_tolerance_;
228  p | min_factor_;
229  p | max_factor_;
230  p | safety_factor_;
231  }
232 
233  private:
234  template <typename EvolvedType, typename ErrorType>
235  double error_calc_impl(const EvolvedType& values,
236  const ErrorType& errors) const noexcept {
237  if constexpr (std::is_fundamental_v<std::remove_cv_t<EvolvedType>> or
238  tt::is_complex_of_fundamental_v<
240  return std::abs(errors) /
241  (absolute_tolerance_ +
242  relative_tolerance_ *
243  std::max(abs(values), abs(values + errors)));
244  } else if constexpr (tt::is_iterable_v<std::remove_cv_t<EvolvedType>>) {
245  double result = 0.0;
246  double recursive_call_result;
247  for (auto val_it = values.begin(), err_it = errors.begin();
248  val_it != values.end(); ++val_it, ++err_it) {
249  recursive_call_result = error_calc_impl(*val_it, *err_it);
250  if (recursive_call_result > result) {
251  result = recursive_call_result;
252  }
253  }
254  return result;
255  } else if constexpr (tt::is_a_v<Variables, std::remove_cv_t<EvolvedType>>) {
256  double result = 0.0;
257  double recursive_call_result;
258  tmpl::for_each<typename EvolvedType::tags_list>(
259  [this, &errors, &values, &recursive_call_result,
260  &result](auto tag_v) noexcept {
261  // clang erroneously warns that `this` is not used despite requiring
262  // it in the capture...
263  (void)this;
264  using tag = typename decltype(tag_v)::type;
265  recursive_call_result = error_calc_impl(
266  get<tag>(values), get<::Tags::StepperError<tag>>(errors));
267  if (recursive_call_result > result) {
268  result = recursive_call_result;
269  }
270  });
271  return result;
272  } else if constexpr (is_any_spin_weighted_v<
274  return error_calc_impl(values.data(), errors.data());
275  }
276  }
277 
278  double absolute_tolerance_ = std::numeric_limits<double>::signaling_NaN();
279  double relative_tolerance_ = std::numeric_limits<double>::signaling_NaN();
280  double max_factor_ = std::numeric_limits<double>::signaling_NaN();
281  double min_factor_ = std::numeric_limits<double>::signaling_NaN();
282  double safety_factor_ = std::numeric_limits<double>::signaling_NaN();
283 };
284 /// \cond
285 template <typename EvolvedVariableTag, typename ErrorControlSelector>
286 PUP::able::PUP_ID
287  ErrorControl<EvolvedVariableTag, ErrorControlSelector>::my_PUP_ID =
288  0; // NOLINT
289 /// \endcond
290 } // namespace StepChoosers
291 
292 namespace Tags {
293 namespace detail {
294 template <typename Tag>
295 using NoPrefix = Tag;
296 } // detail
297 
298 /// \ingroup TimeGroup
299 /// \brief Base tag for reporting whether the `ErrorControl` step chooser is in
300 /// use.
302 
303 /// \ingroup TimeGroup
304 /// \brief A tag that is true if the `ErrorControl` step chooser is one of the
305 /// option-created `Event`s.
306 template <template <typename> typename Prefix = detail::NoPrefix>
309  using type = bool;
310  template <typename Metavariables>
311  using option_tags = tmpl::list<Prefix<::OptionTags::StepChoosers>>;
312 
313  static constexpr bool pass_metavariables = true;
314  template <typename Metavariables>
315  static bool create_from_options(
316  const std::vector<
318  step_choosers) noexcept {
319  bool is_using_error_control = false;
320  // unwraps the factory-created classes to determine whether any of the
321  // step choosers are the ErrorControl step chooser.
322  tmpl::for_each<
323  tmpl::at<typename Metavariables::factory_creation::factory_classes,
325  [&is_using_error_control,
326  &step_choosers](auto step_chooser_type_v) noexcept {
327  if (is_using_error_control) {
328  return;
329  }
330  using step_chooser_type =
331  typename decltype(step_chooser_type_v)::type;
332  if constexpr (tt::is_a_v<::StepChoosers::ErrorControl,
333  step_chooser_type>) {
334  for (const auto& step_chooser : step_choosers) {
335  if (dynamic_cast<const step_chooser_type*>(&*step_chooser) !=
336  nullptr) {
337  is_using_error_control = true;
338  return;
339  }
340  }
341  }
342  });
343  return is_using_error_control;
344  }
345 };
346 
347 /// \ingroup TimeGroup
348 /// \brief Tag indicating that the `ErrorControl` step chooser will not be used.
349 ///
350 /// \details This can be useful when e.g. local time stepping is disabled and it
351 /// is known at compile-time that the `ErrorControl` step chooser cannot be
352 /// used.
354  : db::SimpleTag,
356  using type = bool;
357  using option_tags = tmpl::list<>;
358 
359  static constexpr bool pass_metavariables = false;
360  static bool create_from_options() noexcept { return false; }
361 };
362 } // namespace Tags
NoSuchType
Used to mark "no type" or "bad state" for metaprogramming.
Definition: NoSuchType.hpp:10
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:125
Tags::NeverUsingTimeSteppingErrorControl
Tag indicating that the ErrorControl step chooser will not be used.
Definition: ErrorControl.hpp:353
StepChooser< StepChooserUse::LtsStep >
Definition: StepChooser.hpp:143
std::string
DataBoxTag.hpp
CharmPupable.hpp
pow
constexpr decltype(auto) pow(const T &t) noexcept
Compute t^N where N is an integer (positive or negative)
Definition: ConstantExpressions.hpp:160
StepChoosers
Definition: ByBlock.hpp:33
utility
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
std::pair
StepChoosers::ErrorControl
Suggests a step size based on a target absolute and/or relative error measure.
Definition: ErrorControl.hpp:106
Options.hpp
std::vector
tt::is_a_v
constexpr bool is_a_v
Definition: IsA.hpp:50
StepChoosers::ErrorControl::SafetyFactor
Definition: ErrorControl.hpp:155
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
Tags::IsUsingTimeSteppingErrorControlBase
Base tag for reporting whether the ErrorControl step chooser is in use.
Definition: ErrorControl.hpp:301
WRAPPED_PUPable_decl_template
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
StepChoosers::ErrorControl::AbsoluteTolerance
Definition: ErrorControl.hpp:128
StepChoosers::ErrorControl::MinFactor
Definition: ErrorControl.hpp:147
TimeStepper
Definition: TimeStepper.hpp:47
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
db::BaseTag
Mark a (usually) empty struct as a base tag by inheriting from this.
Definition: Tag.hpp:69
StepChoosers::ErrorControl::MaxFactor
Definition: ErrorControl.hpp:140
Tags::StepperError
Tag for the stepper error measure.
Definition: Tags.hpp:118
Time
Definition: Time.hpp:29
Tags::TimeStepper
Tag for a TimeStepper of type StepperType.
Definition: Tags.hpp:214
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
std::remove_cv_t
std::optional< double >
Tags::StepperErrorUpdated
Tag indicating whether the stepper error has been updated on the current step.
Definition: Tags.hpp:129
StepChoosers::Tags::PreviousStepError
The stepper error measure computed in the previous application of the StepChooser::ErrorControl step ...
Definition: ErrorControl.hpp:37
StepChoosers::ErrorControl::RelativeTolerance
Definition: ErrorControl.hpp:134
tt::is_iterable_v
constexpr bool is_iterable_v
Definition: IsIterable.hpp:51
Tags::IsUsingTimeSteppingErrorControl
A tag that is true if the ErrorControl step chooser is one of the option-created Events.
Definition: ErrorControl.hpp:307
std::unique_ptr
StepChooser
StepChoosers suggest upper bounds on step sizes.
Definition: StepChooser.hpp:56
std::numeric_limits
Parallel
Functionality for parallelization.
Definition: ElementReceiveInterpPoints.hpp:13
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13