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
|
