Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <optional>
8 : #include <string>
9 : #include <tuple>
10 :
11 : #include "DataStructures/DataBox/DataBox.hpp"
12 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
13 : #include "DataStructures/DataBox/Prefixes.hpp"
14 : #include "DataStructures/DataBox/Tag.hpp"
15 : #include "DataStructures/DataBox/TagName.hpp"
16 : #include "Parallel/AlgorithmExecution.hpp"
17 : #include "ParallelAlgorithms/Actions/Goto.hpp"
18 : #include "ParallelAlgorithms/Actions/MutateApply.hpp"
19 : #include "ParallelAlgorithms/Actions/TerminatePhase.hpp"
20 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
21 : #include "Time/AdvanceTime.hpp"
22 : #include "Time/Slab.hpp"
23 : #include "Time/Tags/HistoryEvolvedVariables.hpp"
24 : #include "Time/Tags/TimeStep.hpp"
25 : #include "Time/Time.hpp"
26 : #include "Time/TimeStepId.hpp"
27 : #include "Utilities/Gsl.hpp"
28 : #include "Utilities/TMPL.hpp"
29 : #include "Utilities/TaggedTuple.hpp"
30 : #include "Utilities/TypeTraits/IsA.hpp"
31 :
32 : /// \cond
33 : namespace Parallel {
34 : template <typename Metavariables>
35 : class GlobalCache;
36 : } // namespace Parallel
37 : namespace Tags {
38 : struct TimeStepId;
39 : template <typename StepperInterface>
40 : struct TimeStepper;
41 : } // namespace Tags
42 : class TimeStepper;
43 : namespace TimeSteppers::Tags {
44 : struct FixedOrder;
45 : } // namespace TimeSteppers::Tags
46 : /// \endcond
47 :
48 : /// \ingroup TimeGroup
49 : /// Definition of the integrator self-starting procedure.
50 : ///
51 : /// The self-start procedure generates \f$N\f$ function values of
52 : /// accuracy \f$O\left((\Delta t)^N\right)\f$, where \f$N\f$ is the
53 : /// negative of the initial slab number. To generate values, it
54 : /// requires the time stepper to be a multistep integrator that will
55 : /// produce an order-\f$k\f$-accurate result given \f$k-1\f$ history
56 : /// values.
57 : ///
58 : /// If the integrator is started from analytic history data or
59 : /// requires no history (such as for a substep integrator), then the
60 : /// initial slab number can be set to zero and no self-start steps
61 : /// will be taken.
62 : ///
63 : /// \details
64 : /// To self-start a multistep integrator, the function is integrated
65 : /// repeatedly with increasing accuracy. A first order integrator
66 : /// (Euler's method) requires no history values, so it can be used,
67 : /// starting from the initial point, to generate a
68 : /// first-order-accurate value at a later time. We then reset to the
69 : /// start conditions and use the new "history" value (at a discarded
70 : /// point later in time) to take two steps with a second-order method.
71 : /// These values are second-order-accurate despite the history only
72 : /// being first-order because the calculation of the change in the
73 : /// value multiplies the previous derivatives by a factor of
74 : /// \f$\Delta t\f$. The time and value are then again reset to their
75 : /// starting values and we start again at third order, and so on.
76 : ///
77 : /// The choice of performing the low-order integrations in the same
78 : /// direction as the main integration makes this a _forward_
79 : /// self-start procedure, as opposed to a _backward_ procedure that
80 : /// produces values for times before the start time. The primary
81 : /// advantage of the forward version is that the solution is
82 : /// guaranteed to exist after the start time, but not before. It also
83 : /// makes bookkeeping easier, as the reset after each order increase
84 : /// is to the initial state, rather than to a time one step further
85 : /// back for each order. It does have the disadvantage, however, of
86 : /// leaving a non-monotonic history at the end of the procedure, which
87 : /// the main evolution loop must be able to handle.
88 : ///
89 : /// Each time the state is reset the slab number is increased by one.
90 : /// This ensures that the evaluations are considered to be ordered in
91 : /// their evaluation order, even though they are not monotonic in
92 : /// time. When the slab number reaches zero the initialization
93 : /// procedure is complete and history appropriate for use for an
94 : /// integrator of order \f$N+1\f$ has been generated.
95 : ///
96 : /// The self-start procedure performs all its evaluations before the
97 : /// end of the first time step of the main evolution. It is important
98 : /// that none of the early steps fall at the same time as the
99 : /// self-start history values, so the main evolution should not
100 : /// decrease its step size on the first step after the procedure.
101 : /// Additionally, the history times will not be monotonicly increasing
102 : /// until \f$N\f$ steps have been taken. The local-time-stepping
103 : /// calculations require monotonic time, so local time-stepping should
104 : /// not be initiated until the self-start values have expired from the
105 : /// history. These restrictions on step-size changing are checked in
106 : /// the TimeStepper::can_change_step_size method.
107 1 : namespace SelfStart {
108 : /// Self-start tags
109 1 : namespace Tags {
110 : /// \ingroup TimeGroup
111 : /// The initial value of a quantity. The contents are stored in a
112 : /// tuple to avoid putting duplicate tensors into the DataBox.
113 : template <typename Tag>
114 1 : struct InitialValue : db::PrefixTag, db::SimpleTag {
115 0 : using tag = Tag;
116 0 : using type = std::tuple<typename Tag::type>;
117 : };
118 : } // namespace Tags
119 :
120 : /// Self-start actions
121 1 : namespace Actions {
122 : namespace detail {
123 : template <typename System, bool HasPrimitives>
124 : struct vars_to_save_impl {
125 : using type = tmpl::flatten<tmpl::list<typename System::variables_tag>>;
126 : };
127 :
128 : template <typename System>
129 : struct vars_to_save_impl<System, true> {
130 : using type =
131 : tmpl::flatten<tmpl::list<typename System::variables_tag,
132 : typename System::primitive_variables_tag>>;
133 : };
134 :
135 : template <typename System>
136 : using vars_to_save = typename vars_to_save_impl<
137 : System, System::has_primitive_and_conservative_vars>::type;
138 : } // namespace detail
139 :
140 : /// \ingroup ActionsGroup
141 : /// \ingroup TimeGroup
142 : /// Prepares the evolution for time-stepper self-starting.
143 : ///
144 : /// Stores the initial values of the variables and time step and sets
145 : /// an appropriate step for self-starting.
146 : ///
147 : /// \details The self-start procedure must take place within one slab,
148 : /// and we want to avoid the end of the slab so that we don't have to
149 : /// deal with another action advancing the slab on us. There will be
150 : /// problems if the main evolution tries to evaluate at a time that
151 : /// the self-start procedure used before the self-start version falls
152 : /// out of the history, so we have to make sure that does not happen.
153 : /// We can't do that by making sure our steps are large enough to keep
154 : /// ahead because of the slab restriction, so instead we have to make
155 : /// the self-start step smaller to ensure no collisions. The easiest
156 : /// way to do that is to fit the entire procedure before the first
157 : /// real step, so we pick an initialization time step of
158 : /// \f$\Delta t/(N+1)\f$, for \f$\Delta t\f$ the initial evolution
159 : /// time step and \f$N\f$ the number of history points to be
160 : /// generated.
161 : ///
162 : /// The original `Tags::TimeStep` is temporarily stored separately
163 : /// during the self-start procedure, and restored to its original
164 : /// value at the conclusion of the self-start procedure in preparation
165 : /// for the main evolution, so that the initial time steps during the
166 : /// evolution are appropriately set according to `Initialization`
167 : /// phase values.
168 : ///
169 : /// Uses:
170 : /// - GlobalCache: nothing
171 : /// - DataBox:
172 : /// - Tags::TimeStepId
173 : /// - Tags::TimeStep
174 : /// - variables_tag
175 : /// - primitive_variables_tag if the system has primitives
176 : ///
177 : /// DataBox changes:
178 : /// - Adds:
179 : /// - SelfStart::Tags::InitialValue<Tags::TimeStep>
180 : /// - SelfStart::Tags::InitialValue<variables_tag>
181 : /// - SelfStart::Tags::InitialValue<primitive_variables_tag> if the system
182 : /// has primitives
183 : /// - Removes: nothing
184 : /// - Modifies: Tags::TimeStep
185 : template <typename System>
186 1 : struct Initialize {
187 : private:
188 : template <typename TagsToSave>
189 0 : struct StoreInitialValues;
190 :
191 : template <typename... TagsToSave>
192 0 : struct StoreInitialValues<tmpl::list<TagsToSave...>> {
193 0 : using simple_tags = tmpl::list<Tags::InitialValue<TagsToSave>...>;
194 :
195 : template <typename Box>
196 0 : static void apply(Box& box) {
197 : ::Initialization::mutate_assign<simple_tags>(
198 : make_not_null(&box), std::make_tuple(db::get<TagsToSave>(box))...);
199 : }
200 : };
201 :
202 : public:
203 0 : using simple_tags = typename StoreInitialValues<tmpl::push_back<
204 : detail::vars_to_save<System>, ::Tags::TimeStep>>::simple_tags;
205 :
206 : template <typename DbTags, typename... InboxTags, typename Metavariables,
207 : typename ArrayIndex, typename ActionList,
208 : typename ParallelComponent>
209 0 : static Parallel::iterable_action_return_t apply(
210 : db::DataBox<DbTags>& box,
211 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
212 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
213 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
214 : const ParallelComponent* const /*meta*/) {
215 : const TimeDelta initial_step = db::get<::Tags::TimeStep>(box);
216 : // The slab number increments each time a new point is generated
217 : // until it reaches zero. The multistep integrators all need
218 : // `order` control points at distinct times, even if some of those
219 : // are not counted in the reported required past steps.
220 : // (Predictor stages can't line up with history points.) This
221 : // doesn't count the value as the initial time, hence the "- 1".
222 : const auto values_needed =
223 : db::get<::Tags::Next<::Tags::TimeStepId>>(box).slab_number() == 0
224 : ? 0
225 : : get<TimeSteppers::Tags::FixedOrder>(
226 : db::get<::Tags::TimeStepper<TimeStepper>>(box).order()) -
227 : 1;
228 :
229 : // Decrease the step so that the generated history will be
230 : // entirely before the first step. This ensures we will not
231 : // generate any duplicate times in the history as we start the
232 : // real evolution and that the starting procedure does not require
233 : // any more information (such as function-of-time values) than the
234 : // first real step.
235 : const TimeDelta self_start_step = initial_step / (values_needed + 1);
236 :
237 : StoreInitialValues<tmpl::push_back<detail::vars_to_save<System>,
238 : ::Tags::TimeStep>>::apply(box);
239 : db::mutate<::Tags::TimeStep>(
240 : [&self_start_step](const gsl::not_null<::TimeDelta*> time_step) {
241 : *time_step = self_start_step;
242 : },
243 : make_not_null(&box));
244 :
245 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
246 : }
247 : };
248 :
249 : /// \ingroup ActionsGroup
250 : /// \ingroup TimeGroup
251 : /// Resets the state for the next iteration if the current order is
252 : /// complete, and exits the self-start loop if the required order has
253 : /// been reached.
254 : ///
255 : /// Uses:
256 : /// - GlobalCache: nothing
257 : /// - DataBox: Tags::Next<Tags::TimeStepId>
258 : ///
259 : /// DataBox changes:
260 : /// - Adds: nothing
261 : /// - Removes: nothing
262 : /// - Modifies: nothing
263 : template <typename ExitTag, typename System>
264 1 : struct CheckForCompletion {
265 : template <typename DbTags, typename... InboxTags, typename Metavariables,
266 : typename ArrayIndex, typename ActionList,
267 : typename ParallelComponent>
268 0 : static Parallel::iterable_action_return_t apply(
269 : db::DataBox<DbTags>& box,
270 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
271 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
272 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
273 : const ParallelComponent* const /*meta*/) {
274 : const bool done_with_order =
275 : db::get<::Tags::Next<::Tags::TimeStepId>>(box).is_at_slab_boundary();
276 :
277 : if (done_with_order) {
278 : tmpl::for_each<detail::vars_to_save<System>>([&box](auto tag) {
279 : using Tag = tmpl::type_from<decltype(tag)>;
280 : db::mutate<Tag>(
281 : [](const gsl::not_null<typename Tag::type*> value,
282 : const std::tuple<typename Tag::type>& initial_value) {
283 : *value = get<0>(initial_value);
284 : },
285 : make_not_null(&box), db::get<Tags::InitialValue<Tag>>(box));
286 : });
287 : }
288 :
289 : // The self start procedure begins with slab number
290 : // -number_of_past_steps and counts up. When we reach 0 we should
291 : // start the evolution proper. The first thing the evolution loop
292 : // will do is update the time, so here we need to check if the
293 : // next time should be the first real step.
294 : if (db::get<::Tags::Next<::Tags::TimeStepId>>(box).slab_number() == 0) {
295 : return {Parallel::AlgorithmExecution::Continue,
296 : tmpl::index_of<ActionList, ::Actions::Label<ExitTag>>::value};
297 : }
298 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
299 : }
300 : };
301 :
302 : /// \ingroup ActionsGroup
303 : /// \ingroup TimeGroup
304 : /// If we have taken enough steps for this order, set the next time to
305 : /// the start time and increment the slab number
306 : ///
307 : /// Uses:
308 : /// - GlobalCache: nothing
309 : /// - DataBox:
310 : /// - Tags::HistoryEvolvedVariables
311 : /// - Tags::TimeStepId
312 : /// - Tags::TimeStep
313 : ///
314 : /// DataBox changes:
315 : /// - Adds: nothing
316 : /// - Removes: nothing
317 : /// - Modifies: Tags::Next<Tags::TimeStepId> if there is an order increase
318 1 : struct CheckForOrderIncrease {
319 : template <typename DbTags, typename... InboxTags, typename Metavariables,
320 : typename ArrayIndex, typename ActionList,
321 : typename ParallelComponent>
322 0 : static Parallel::iterable_action_return_t apply(
323 : db::DataBox<DbTags>& box,
324 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
325 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
326 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
327 : const ParallelComponent* const /*meta*/) { // NOLINT const
328 : const auto& time_step_id = db::get<::Tags::TimeStepId>(box);
329 : const auto& time = time_step_id.step_time();
330 : const auto& time_step = db::get<::Tags::TimeStep>(box);
331 : using history_tags = ::Tags::get_all_history_tags<DbTags>;
332 : const size_t history_integration_order =
333 : db::get<tmpl::front<history_tags>>(box).integration_order();
334 :
335 : // This is correct for both AdamsBashforth and AdamsMoultonPc.
336 : // The latter doesn't need as many history values, but it can't
337 : // use the first step of the previous order because it would line
338 : // up with the predictor stage of the first step for the current
339 : // order.
340 : const Time required_time =
341 : (time_step.is_positive() ? time.slab().start() : time.slab().end()) +
342 : history_integration_order * time_step;
343 : const bool done_with_order =
344 : time_step_id.substep() == 0 and time == required_time;
345 :
346 : if (done_with_order) {
347 : db::mutate<::Tags::Next<::Tags::TimeStepId>>(
348 : [](const gsl::not_null<::TimeStepId*> next_time_id,
349 : const ::TimeStepId& current_time_id) {
350 : const Slab slab = current_time_id.step_time().slab();
351 : *next_time_id =
352 : TimeStepId(current_time_id.time_runs_forward(),
353 : current_time_id.slab_number() + 1,
354 : current_time_id.time_runs_forward() ? slab.start()
355 : : slab.end());
356 : },
357 : make_not_null(&box), db::get<::Tags::TimeStepId>(box));
358 : tmpl::for_each<history_tags>([&box,
359 : &history_integration_order](auto tag_v) {
360 : using history_tag = typename decltype(tag_v)::type;
361 : db::mutate<history_tag>(
362 : [&history_integration_order](
363 : const gsl::not_null<typename history_tag::type*>
364 : mutable_history) {
365 : ASSERT(mutable_history->integration_order() ==
366 : history_integration_order,
367 : "Using multiple histories of different integration "
368 : "orders. When using multiple histories sharing the same "
369 : "time stepper, all histories must maintain identical "
370 : "integration order.");
371 : mutable_history->integration_order(
372 : mutable_history->integration_order() + 1);
373 : },
374 : make_not_null(&box));
375 : });
376 : }
377 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
378 : }
379 : };
380 :
381 : /// \ingroup ActionsGroup
382 : /// \ingroup TimeGroup
383 : /// Cleans up after the self-start procedure
384 : ///
385 : /// Resets the time step to that requested for the evolution and
386 : /// removes temporary self-start data.
387 : ///
388 : /// Uses:
389 : /// - GlobalCache: nothing
390 : /// - DataBox: SelfStart::Tags::InitialValue<Tags::TimeStep>
391 : ///
392 : /// DataBox changes:
393 : /// - Adds: nothing
394 : /// - Removes:
395 : /// - All SelfStart::Tags::InitialValue tags
396 : /// - Modifies: Tags::TimeStep
397 1 : struct Cleanup {
398 : private:
399 : template <typename T>
400 0 : struct is_a_initial_value : tt::is_a<Tags::InitialValue, T> {};
401 :
402 0 : using initial_step_tag = Tags::InitialValue<::Tags::TimeStep>;
403 :
404 : public:
405 : template <typename DbTags, typename... InboxTags, typename Metavariables,
406 : typename ArrayIndex, typename ActionList,
407 : typename ParallelComponent>
408 0 : static Parallel::iterable_action_return_t apply(
409 : db::DataBox<DbTags>& box,
410 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
411 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
412 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
413 : const ParallelComponent* const /*meta*/) {
414 : // Reset the time step to the value requested by the user. The
415 : // variables were reset in CheckForCompletion.
416 : db::mutate<::Tags::TimeStep>(
417 : [](const gsl::not_null<::TimeDelta*> time_step,
418 : const std::tuple<::TimeDelta>& initial_step) {
419 : *time_step = get<0>(initial_step);
420 : },
421 : make_not_null(&box), db::get<initial_step_tag>(box));
422 : using remove_tags = tmpl::filter<DbTags, is_a_initial_value<tmpl::_1>>;
423 : // reset each tag to default constructed values to reduce memory usage (Data
424 : // structures like `DataVector`s and `Tensor`s have negligible memory usage
425 : // when default-constructed). This tends to be better for build time than
426 : // constructing more DataBox types with and without this handful of tags.
427 : tmpl::for_each<remove_tags>([&box](auto tag_v) {
428 : using tag = typename decltype(tag_v)::type;
429 : db::mutate<tag>(
430 : [](auto tag_value) {
431 : *tag_value = std::decay_t<decltype(*tag_value)>{};
432 : },
433 : make_not_null(&box));
434 : });
435 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
436 : }
437 : };
438 : } // namespace Actions
439 :
440 : namespace detail {
441 : struct PhaseStart;
442 : struct PhaseEnd;
443 : } // namespace detail
444 :
445 : /// \ingroup TimeGroup
446 : /// The list of actions required to self-start an integrator.
447 : ///
448 : /// \tparam StepActions List of actions computing and recording the
449 : /// system derivative and updating the evolved variables (but not the
450 : /// time).
451 : ///
452 : /// \see SelfStart
453 : template <typename StepActions, typename System>
454 1 : using self_start_procedure = tmpl::flatten<tmpl::list<
455 : // clang-format off
456 : SelfStart::Actions::Initialize<System>,
457 : ::Actions::Label<detail::PhaseStart>,
458 : SelfStart::Actions::CheckForCompletion<detail::PhaseEnd, System>,
459 : ::Actions::MutateApply<AdvanceTime>,
460 : SelfStart::Actions::CheckForOrderIncrease,
461 : StepActions,
462 : ::Actions::Goto<detail::PhaseStart>,
463 : ::Actions::Label<detail::PhaseEnd>,
464 : SelfStart::Actions::Cleanup,
465 : ::Actions::MutateApply<AdvanceTime>,
466 : Parallel::Actions::TerminatePhase>>;
467 : // clang-format on
468 : } // namespace SelfStart
|
