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