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/Registration.hpp"
17 #include "Utilities/TMPL.hpp"
18 #include "Utilities/TypeTraits/IsIterable.hpp"
19 
20 /// \cond
21 class Time;
22 namespace Parallel {
23 template <typename Metavariables>
24 class GlobalCache;
25 } // namespace Parallel
26 /// \endcond
27 
28 namespace StepChoosers {
29 template <typename EvolvedType, typename ErrorControlSelector,
30  typename StepChooserRegistrars>
32 
33 namespace Registrars {
34 template <typename EvolvedType, typename ErrorControlSelector = NoSuchType>
36  EvolvedType, ErrorControlSelector>;
37 } // namespace Registrars
38 
39 namespace Tags {
42 };
43 } // namespace Tags
44 
45 /*!
46  * \brief Suggests a step size based on a target absolute and/or relative error
47  * measure.
48  *
49  * \details The suggested step is calculated via a simple specialization of the
50  * scheme suggested in \cite Hairer1993. We first compute the aggregated error
51  * measure from the stepper error:
52  *
53  * \f[
54  * E = \max_i(|E_i| / sc_i),
55  * \f]
56  *
57  * where \f$E_i\f$ is the ODE error reported for each individual grid point,
58  * reported by the time stepper, and \f$sc_i\f$ is the step control measure
59  * determined by the tolerances:
60  *
61  * \f[
62  * sc_i = Atol_i + \max(|y_i|,|y_i + E_i|) Rtol_i,
63  * \f]
64  *
65  * and \f$y_i\f$ is the value of the function at the current step at grid point
66  * \f$i\f$.
67  *
68  * When no previous record of previous error is available, the step has size:
69  *
70  * \f[
71  * h_{\text{new}} = h \cdot \min\left(F_{\text{max}},
72  * \max\left(F_{\text{min}},
73  * \frac{F_{\text{safety}}}{E^{1/(q + 1)}}\right)\right),
74  * \f]
75  *
76  * where \f$h_{\text{new}}\f$ is the new suggested step size \f$h\f$ is the
77  * previous step size, \f$F_{\text{max}}\f$ is the maximum factor by which we
78  * allow the step to increase, \f$F_{\text{min}}\f$ is the minimum factor by
79  * which we allow the step to decrease. \f$F_{\text{safety}}\f$ is the safety
80  * factor on the computed error -- this forces the step size slightly lower
81  * than we would naively compute so that the result of the step will likely be
82  * within the target error. \f$q\f$ is the order of the stepper error
83  * calculation. Intuitively, we should change the step less drastically for a
84  * higher order stepper.
85  *
86  * After the first error calculation, the error \f$E\f$ is recorded, and
87  * subsequent error calculations use a simple PI scheme suggested in
88  * \cite NumericalRecipes section 17.2.1:
89  *
90  * \f[
91  * h_{\text{new}} = h \cdot \min\left(F_{\text{max}},
92  * \max\left(F_{\text{min}},
93  * \frac{F_{\text{safety}}}{E^{0.7 / (q + 1)}
94  * E_{\text{prev}}^{0.4 / (q + 1)}}\right)\right),
95  * \f]
96  *
97  * where \f$E_{\text{prev}}\f$ is the error computed in the previous step.
98  *
99  * \note The template parameter `ErrorControlSelector` is used to disambiguate
100  * in the input-file options between `ErrorControl` step choosers that are
101  * based on different variables. This is needed if multiple systems are evolved
102  * in the same executable. The name used for the input file includes
103  * `ErrorControlSelector::name()` if it is provided.
104  */
105 template <
106  typename EvolvedVariableTag, typename ErrorControlSelector = NoSuchType,
107  typename StepChooserRegistrars = tmpl::list<
109 class ErrorControl : public StepChooser<StepChooserRegistrars> {
110  public:
111  using evolved_variable_type = typename EvolvedVariableTag::type;
112  using error_variable_type =
114  EvolvedVariableTag>::type;
115 
116  /// \cond
117  ErrorControl() = default;
118  explicit ErrorControl(CkMigrateMessage* /*unused*/) noexcept {}
119  using PUP::able::register_constructor;
120  WRAPPED_PUPable_decl_template(ErrorControl); // NOLINT
121  /// \endcond
122 
123  static std::string name() noexcept {
124  if constexpr (std::is_same_v<ErrorControlSelector, NoSuchType>) {
125  return "ErrorControl";
126  } else {
127  return "ErrorControl(" + ErrorControlSelector::name() + ")";
128  }
129  }
130 
132  using type = double;
133  static constexpr Options::String help{"Target absolute tolerance"};
134  static type lower_bound() noexcept { return 0.0; }
135  };
136 
138  using type = double;
139  static constexpr Options::String help{"Target relative tolerance"};
140  static type lower_bound() noexcept { return 0.0; }
141  };
142 
143  struct MaxFactor {
144  using type = double;
145  static constexpr Options::String help{
146  "Maximum factor to increase the step by"};
147  static type lower_bound() noexcept { return 1.0; }
148  };
149 
150  struct MinFactor {
151  using type = double;
152  static constexpr Options::String help{
153  "Minimum factor to increase the step by"};
154  static type lower_bound() noexcept { return 0.0; }
155  static type upper_bound() noexcept { return 1.0; }
156  };
157 
158  struct SafetyFactor {
159  using type = double;
160  static constexpr Options::String help{
161  "Extra factor to apply to step estimate; can be used to decrease step "
162  "size to improve step acceptance rate."};
163  static type lower_bound() noexcept { return 0.0; }
164  };
165 
166  static constexpr Options::String help{
167  "Chooses a step based on a target relative and absolute error tolerance"};
168  using options = tmpl::list<AbsoluteTolerance, RelativeTolerance, MaxFactor,
169  MinFactor, SafetyFactor>;
170 
171  ErrorControl(const double absolute_tolerance, const double relative_tolerance,
172  const double max_factor, const double min_factor,
173  const double safety_factor) noexcept
174  : absolute_tolerance_{absolute_tolerance},
175  relative_tolerance_{relative_tolerance},
176  max_factor_{max_factor},
177  min_factor_{min_factor},
178  safety_factor_{safety_factor} {}
179 
180  static constexpr UsableFor usable_for = UsableFor::OnlyLtsStepChoice;
181 
182  using argument_tags =
183  tmpl::list<::Tags::HistoryEvolvedVariables<EvolvedVariableTag>,
186 
187  using return_tags = tmpl::list<Tags::PreviousStepError>;
188 
189  using simple_tags = tmpl::list<Tags::PreviousStepError>;
190 
191  template <typename Metavariables, typename History, typename TimeStepper>
192  std::pair<double, bool> operator()(
193  const gsl::not_null<std::optional<double>*> previous_step_error,
194  const History& history, const error_variable_type& error,
195  const bool& stepper_error_updated, const TimeStepper& stepper,
196  const double previous_step,
197  const Parallel::GlobalCache<Metavariables>& /*cache*/) const noexcept {
198  // request that the step size not be changed if there isn't a new error
199  // estimate
200  if (not stepper_error_updated) {
201  return std::make_pair(std::numeric_limits<double>::infinity(), true);
202  }
203  const double l_inf_error =
204  error_calc_impl((history.end() - 1).value(), error);
205  double new_step;
206  if (not previous_step_error->has_value()) {
207  new_step = previous_step *
208  std::clamp(safety_factor_ *
209  pow(1.0 / std::max(l_inf_error, 1e-14),
210  1.0 / (stepper.error_estimate_order() + 1)),
211  min_factor_, max_factor_);
212  } else {
213  // From simple advice from Numerical Recipes 17.2.1 regarding a heuristic
214  // for PI step control.
215  const double alpha_factor = 0.7 / (stepper.error_estimate_order() + 1);
216  const double beta_factor = 0.4 / (stepper.error_estimate_order() + 1);
217  new_step =
218  previous_step *
219  std::clamp(
220  safety_factor_ *
221  pow(1.0 / std::max(l_inf_error, 1e-14), alpha_factor) *
222  pow(1.0 / std::max(previous_step_error->value(), 1e-14),
223  beta_factor),
224  min_factor_, max_factor_);
225  }
226  *previous_step_error = l_inf_error;
227  return std::make_pair(new_step, l_inf_error <= 1.0);
228  }
229 
230  void pup(PUP::er& p) noexcept override { // NOLINT
231  p | absolute_tolerance_;
232  p | relative_tolerance_;
233  p | min_factor_;
234  p | max_factor_;
235  p | safety_factor_;
236  }
237 
238  private:
239  template <typename EvolvedType, typename ErrorType>
240  double error_calc_impl(const EvolvedType& values,
241  const ErrorType& errors) const noexcept {
242  if constexpr (std::is_fundamental_v<std::remove_cv_t<EvolvedType>> or
243  tt::is_complex_of_fundamental_v<
245  return std::abs(errors) /
246  (absolute_tolerance_ +
247  relative_tolerance_ *
248  std::max(abs(values), abs(values + errors)));
249  } else if constexpr (tt::is_iterable_v<std::remove_cv_t<EvolvedType>>) {
250  double result = 0.0;
251  double recursive_call_result;
252  for (auto val_it = values.begin(), err_it = errors.begin();
253  val_it != values.end(); ++val_it, ++err_it) {
254  recursive_call_result = error_calc_impl(*val_it, *err_it);
255  if (recursive_call_result > result) {
256  result = recursive_call_result;
257  }
258  }
259  return result;
260  } else if constexpr (tt::is_a_v<Variables, std::remove_cv_t<EvolvedType>>) {
261  double result = 0.0;
262  double recursive_call_result;
263  tmpl::for_each<typename EvolvedType::tags_list>(
264  [this, &errors, &values, &recursive_call_result,
265  &result](auto tag_v) noexcept {
266  // clang erroneously warns that `this` is not used despite requiring
267  // it in the capture...
268  (void)this;
269  using tag = typename decltype(tag_v)::type;
270  recursive_call_result = error_calc_impl(
271  get<tag>(values), get<::Tags::StepperError<tag>>(errors));
272  if (recursive_call_result > result) {
273  result = recursive_call_result;
274  }
275  });
276  return result;
277  }
278  }
279 
280  double absolute_tolerance_ = std::numeric_limits<double>::signaling_NaN();
281  double relative_tolerance_ = std::numeric_limits<double>::signaling_NaN();
282  double max_factor_ = std::numeric_limits<double>::signaling_NaN();
283  double min_factor_ = std::numeric_limits<double>::signaling_NaN();
284  double safety_factor_ = std::numeric_limits<double>::signaling_NaN();
285 };
286 /// \cond
287 template <typename EvolvedVariableTag, typename ErrorControlSelector,
288  typename StepChooserRegistrars>
289 PUP::able::PUP_ID ErrorControl<EvolvedVariableTag, ErrorControlSelector,
290  StepChooserRegistrars>::my_PUP_ID = 0; // NOLINT
291 /// \endcond
292 } // namespace StepChoosers
293 
294 namespace Tags {
295 /// Determines at runtime whether the `ErrorControl` step chooser has been
296 /// selected.
297 template <typename StepChooserRegistrars>
299  db::SimpleTag {
300  using type = bool;
301  template <typename Metavariables>
302  using option_tags = tmpl::flatten<tmpl::list<
303  tmpl::conditional_t<(tmpl::size<StepChooserRegistrars>::value > 0),
305  tmpl::list<>>>>;
306 
307  static constexpr bool pass_metavariables = true;
308  template <typename Metavariables>
309  static bool create_from_options(
311  step_choosers) noexcept {
312  bool is_using_error_control = false;
313  // unwraps the factory-created classes to determine whether any of the
314  // step choosers are the ErrorControl step chooser.
315  tmpl::for_each<
316  typename StepChooser<StepChooserRegistrars>::creatable_classes>(
317  [&is_using_error_control,
318  &step_choosers](auto step_chooser_type_v) noexcept {
319  if (is_using_error_control) {
320  return;
321  }
322  using step_chooser_type =
323  typename decltype(step_chooser_type_v)::type;
324  if constexpr (tt::is_a_v<::StepChoosers::ErrorControl,
325  step_chooser_type>) {
326  for (const auto& step_chooser : step_choosers) {
327  if (dynamic_cast<const step_chooser_type*>(&*step_chooser) !=
328  nullptr) {
329  is_using_error_control = true;
330  return;
331  }
332  }
333  }
334  });
335  return is_using_error_control;
336  }
337 };
338 
339 template <>
342  using type = bool;
343  using option_tags = tmpl::list<>;
344  static constexpr bool pass_metavariables = false;
345  static bool create_from_options() noexcept { return false; }
346 };
347 } // namespace Tags
NoSuchType
Used to mark "no type" or "bad state" for metaprogramming.
Definition: NoSuchType.hpp:10
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
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
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:31
Options.hpp
std::vector
tt::is_a_v
constexpr bool is_a_v
Definition: IsA.hpp:50
StepChoosers::ErrorControl::SafetyFactor
Definition: ErrorControl.hpp:158
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
OptionTags::StepChoosers
Definition: Tags.hpp:151
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
Tags::IsUsingTimeSteppingErrorControlBase
A tag that is true if the ErrorControl step chooser is one of the option-created Events.
Definition: Tags.hpp:256
Registration::Registrar
A template for defining a registrar.
Definition: Registration.hpp:42
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:131
StepChoosers::ErrorControl::MinFactor
Definition: ErrorControl.hpp:150
TimeStepper
Definition: TimeStepper.hpp:47
StepChoosers::UsableFor
UsableFor
Designation for the context in which a step chooser may be used.
Definition: StepChooser.hpp:54
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
StepChoosers::ErrorControl::MaxFactor
Definition: ErrorControl.hpp:143
Tags::StepperError
Tag for the stepper error measure.
Definition: Tags.hpp:110
Time
Definition: Time.hpp:29
Tags::TimeStepper
Tag for a TimeStepper of type StepperType.
Definition: Tags.hpp:206
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:121
StepChoosers::Tags::PreviousStepError
Definition: ErrorControl.hpp:40
StepChoosers::ErrorControl::RelativeTolerance
Definition: ErrorControl.hpp:137
tt::is_iterable_v
constexpr bool is_iterable_v
Definition: IsIterable.hpp:51
Tags::IsUsingTimeSteppingErrorControl
Determines at runtime whether the ErrorControl step chooser has been selected.
Definition: ErrorControl.hpp:298
std::unique_ptr
StepChooser
Definition: StepChooser.hpp:98
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: ReadSpecThirdOrderPiecewisePolynomial.hpp:13