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" // IWYU pragma: keep
18 : #include "ParallelAlgorithms/Actions/TerminatePhase.hpp" // IWYU pragma: keep
19 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
20 : #include "Time/Actions/AdvanceTime.hpp" // IWYU pragma: keep
21 : #include "Time/Slab.hpp"
22 : #include "Time/Tags/HistoryEvolvedVariables.hpp"
23 : #include "Time/Tags/TimeStep.hpp"
24 : #include "Time/Time.hpp"
25 : #include "Time/TimeStepId.hpp"
26 : #include "Utilities/Gsl.hpp"
27 : #include "Utilities/TMPL.hpp"
28 : #include "Utilities/TaggedTuple.hpp"
29 : #include "Utilities/TypeTraits/IsA.hpp"
30 :
31 : // IWYU pragma: no_include "DataStructures/Tensor/Tensor.hpp"
32 :
33 : /// \cond
34 : namespace Parallel {
35 : template <typename Metavariables>
36 : class GlobalCache;
37 : } // namespace Parallel
38 : namespace Tags {
39 : struct TimeStepId;
40 : template <typename StepperInterface>
41 : struct TimeStepper;
42 : } // namespace Tags
43 : class TimeStepper;
44 : /// \endcond
45 :
46 : /// \ingroup TimeGroup
47 : /// Definition of the integrator self-starting procedure.
48 : ///
49 : /// The self-start procedure generates \f$N\f$ function values of
50 : /// accuracy \f$O\left((\Delta t)^N\right)\f$, where \f$N\f$ is the
51 : /// negative of the initial slab number. To generate values, it
52 : /// requires the time stepper to be a multistep integrator that will
53 : /// produce an order-\f$k\f$-accurate result given \f$k-1\f$ history
54 : /// values.
55 : ///
56 : /// If the integrator is started from analytic history data or
57 : /// requires no history (such as for a substep integrator), then the
58 : /// initial slab number can be set to zero and no self-start steps
59 : /// will be taken.
60 : ///
61 : /// \details
62 : /// To self-start a multistep integrator, the function is integrated
63 : /// repeatedly with increasing accuracy. A first order integrator
64 : /// (Euler's method) requires no history values, so it can be used,
65 : /// starting from the initial point, to generate a
66 : /// first-order-accurate value at a later time. We then reset to the
67 : /// start conditions and use the new "history" value (at a discarded
68 : /// point later in time) to take two steps with a second-order method.
69 : /// These values are second-order-accurate despite the history only
70 : /// being first-order because the calculation of the change in the
71 : /// value multiplies the previous derivatives by a factor of
72 : /// \f$\Delta t\f$. The time and value are then again reset to their
73 : /// starting values and we start again at third order, and so on.
74 : ///
75 : /// The choice of performing the low-order integrations in the same
76 : /// direction as the main integration makes this a _forward_
77 : /// self-start procedure, as opposed to a _backward_ procedure that
78 : /// produces values for times before the start time. The primary
79 : /// advantage of the forward version is that the solution is
80 : /// guaranteed to exist after the start time, but not before. It also
81 : /// makes bookkeeping easier, as the reset after each order increase
82 : /// is to the initial state, rather than to a time one step further
83 : /// back for each order. It does have the disadvantage, however, of
84 : /// leaving a non-monotonic history at the end of the procedure, which
85 : /// the main evolution loop must be able to handle.
86 : ///
87 : /// Each time the state is reset the slab number is increased by one.
88 : /// This ensures that the evaluations are considered to be ordered in
89 : /// their evaluation order, even though they are not monotonic in
90 : /// time. When the slab number reaches zero the initialization
91 : /// procedure is complete and history appropriate for use for an
92 : /// integrator of order \f$N+1\f$ has been generated.
93 : ///
94 : /// The self-start procedure performs all its evaluations before the
95 : /// end of the first time step of the main evolution. It is important
96 : /// that none of the early steps fall at the same time as the
97 : /// self-start history values, so the main evolution should not
98 : /// decrease its step size on the first step after the procedure.
99 : /// Additionally, the history times will not be monotonicly increasing
100 : /// until \f$N\f$ steps have been taken. The local-time-stepping
101 : /// calculations require monotonic time, so local time-stepping should
102 : /// not be initiated until the self-start values have expired from the
103 : /// history. These restrictions on step-size changing are checked in
104 : /// the TimeStepper::can_change_step_size method.
105 1 : namespace SelfStart {
106 : /// Self-start tags
107 1 : namespace Tags {
108 : /// \ingroup TimeGroup
109 : /// The initial value of a quantity. The contents are stored in a
110 : /// tuple to avoid putting duplicate tensors into the DataBox.
111 : template <typename Tag>
112 1 : struct InitialValue : db::PrefixTag, db::SimpleTag {
113 0 : using tag = Tag;
114 0 : using type = std::tuple<typename Tag::type>;
115 : };
116 : } // namespace Tags
117 :
118 : /// Self-start actions
119 1 : namespace Actions {
120 : namespace detail {
121 : template <typename System, bool HasPrimitives>
122 : struct vars_to_save_impl {
123 : using type = tmpl::flatten<tmpl::list<typename System::variables_tag>>;
124 : };
125 :
126 : template <typename System>
127 : struct vars_to_save_impl<System, true> {
128 : using type =
129 : tmpl::flatten<tmpl::list<typename System::variables_tag,
130 : typename System::primitive_variables_tag>>;
131 : };
132 :
133 : template <typename System>
134 : using vars_to_save = typename vars_to_save_impl<
135 : System, System::has_primitive_and_conservative_vars>::type;
136 : } // namespace detail
137 :
138 : /// \ingroup ActionsGroup
139 : /// \ingroup TimeGroup
140 : /// Prepares the evolution for time-stepper self-starting.
141 : ///
142 : /// Stores the initial values of the variables and time step and sets
143 : /// an appropriate step for self-starting.
144 : ///
145 : /// \details The self-start procedure must take place within one slab,
146 : /// and we want to avoid the end of the slab so that we don't have to
147 : /// deal with another action advancing the slab on us. There will be
148 : /// problems if the main evolution tries to evaluate at a time that
149 : /// the self-start procedure used before the self-start version falls
150 : /// out of the history, so we have to make sure that does not happen.
151 : /// We can't do that by making sure our steps are large enough to keep
152 : /// ahead because of the slab restriction, so instead we have to make
153 : /// the self-start step smaller to ensure no collisions. The easiest
154 : /// way to do that is to fit the entire procedure before the first
155 : /// real step, so we pick an initialization time step of
156 : /// \f$\Delta t/(N+1)\f$, for \f$\Delta t\f$ the initial evolution
157 : /// time step and \f$N\f$ the number of history points to be
158 : /// generated.
159 : ///
160 : /// The original `Tags::TimeStep` and `Tags::Next<Tags::TimeStep>` are
161 : /// temporarily stored separately during the self-start procedure, and restored
162 : /// to their original values at the conclusion of the self-start procedure in
163 : /// preparation for the main evolution, so that the initial time steps during
164 : /// the evolution are appropriately set according to `Initialization` phase
165 : /// values.
166 : ///
167 : /// Uses:
168 : /// - GlobalCache: nothing
169 : /// - DataBox:
170 : /// - Tags::TimeStepId
171 : /// - Tags::TimeStep
172 : /// - variables_tag
173 : /// - primitive_variables_tag if the system has primitives
174 : ///
175 : /// DataBox changes:
176 : /// - Adds:
177 : /// - SelfStart::Tags::InitialValue<Tags::TimeStep>
178 : /// - SelfStart::Tags::InitialValue<variables_tag>
179 : /// - SelfStart::Tags::InitialValue<primitive_variables_tag> if the system
180 : /// has primitives
181 : /// - Removes: nothing
182 : /// - Modifies: Tags::TimeStep
183 : template <typename System>
184 1 : struct Initialize {
185 : private:
186 : template <typename TagsToSave>
187 0 : struct StoreInitialValues;
188 :
189 : template <typename... TagsToSave>
190 0 : struct StoreInitialValues<tmpl::list<TagsToSave...>> {
191 0 : using simple_tags = tmpl::list<Tags::InitialValue<TagsToSave>...>;
192 :
193 : template <typename Box>
194 0 : static void apply(Box& box) {
195 : ::Initialization::mutate_assign<simple_tags>(
196 : make_not_null(&box), std::make_tuple(db::get<TagsToSave>(box))...);
197 : }
198 : };
199 :
200 : public:
201 0 : using simple_tags = typename StoreInitialValues<
202 : tmpl::push_back<detail::vars_to_save<System>, ::Tags::TimeStep,
203 : ::Tags::Next<::Tags::TimeStep>>>::simple_tags;
204 :
205 : template <typename DbTags, typename... InboxTags, typename Metavariables,
206 : typename ArrayIndex, typename ActionList,
207 : typename ParallelComponent>
208 0 : static Parallel::iterable_action_return_t apply(
209 : db::DataBox<DbTags>& box,
210 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
211 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
212 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
213 : const ParallelComponent* const /*meta*/) {
214 : const TimeDelta initial_step = db::get<::Tags::TimeStep>(box);
215 : // The slab number increments each time a new point is generated
216 : // until it reaches zero. The multistep integrators all need
217 : // `order` control points at distinct times, even if some of those
218 : // are not counted in the reported required past steps.
219 : // (Predictor stages can't line up with history points.) This
220 : // doesn't count the value as the initial time, hence the "- 1".
221 : const auto values_needed =
222 : db::get<::Tags::Next<::Tags::TimeStepId>>(box).slab_number() == 0
223 : ? 0
224 : : db::get<::Tags::TimeStepper<TimeStepper>>(box).order() - 1;
225 :
226 : // Decrease the step so that the generated history will be
227 : // entirely before the first step. This ensures we will not
228 : // generate any duplicate times in the history as we start the
229 : // real evolution and that the starting procedure does not require
230 : // any more information (such as function-of-time values) than the
231 : // first real step.
232 : const TimeDelta self_start_step = initial_step / (values_needed + 1);
233 :
234 : StoreInitialValues<
235 : tmpl::push_back<detail::vars_to_save<System>, ::Tags::TimeStep,
236 : ::Tags::Next<::Tags::TimeStep>>>::apply(box);
237 : db::mutate<::Tags::TimeStep, ::Tags::Next<::Tags::TimeStep>>(
238 : [&self_start_step](const gsl::not_null<::TimeDelta*> time_step,
239 : const gsl::not_null<::TimeDelta*> next_time_step) {
240 : *time_step = self_start_step;
241 : *next_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 : // The time steppers do cleanup internally when taking a
374 : // step, but since we are resetting to the start instead
375 : // of taking a step we have to do it here. If we
376 : // skipped this the extra value would harmlessly fall
377 : // off the end after a few steps, but that would
378 : // unnecessarily increase the size of the data held by
379 : // the History object. A different way to handle this
380 : // would be to call shrink_to_fit() on the history a few
381 : // steps after the start of the evolution.
382 : for (const auto& record : *mutable_history) {
383 : mutable_history->discard_value(record.time_step_id);
384 : }
385 : },
386 : make_not_null(&box));
387 : });
388 : }
389 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
390 : }
391 : };
392 :
393 : /// \ingroup ActionsGroup
394 : /// \ingroup TimeGroup
395 : /// Cleans up after the self-start procedure
396 : ///
397 : /// Resets the time step to that requested for the evolution and
398 : /// removes temporary self-start data.
399 : ///
400 : /// Uses:
401 : /// - GlobalCache: nothing
402 : /// - DataBox: SelfStart::Tags::InitialValue<Tags::TimeStep>
403 : ///
404 : /// DataBox changes:
405 : /// - Adds: nothing
406 : /// - Removes:
407 : /// - All SelfStart::Tags::InitialValue tags
408 : /// - Modifies: Tags::TimeStep
409 1 : struct Cleanup {
410 : private:
411 : template <typename T>
412 0 : struct is_a_initial_value : tt::is_a<Tags::InitialValue, T> {};
413 :
414 0 : using initial_step_tag = Tags::InitialValue<::Tags::TimeStep>;
415 :
416 : public:
417 : template <typename DbTags, typename... InboxTags, typename Metavariables,
418 : typename ArrayIndex, typename ActionList,
419 : typename ParallelComponent>
420 0 : static Parallel::iterable_action_return_t apply(
421 : db::DataBox<DbTags>& box,
422 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
423 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
424 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
425 : const ParallelComponent* const /*meta*/) {
426 : // Reset the time step to the value requested by the user. The
427 : // variables were reset in CheckForCompletion.
428 : db::mutate<::Tags::TimeStep, ::Tags::Next<::Tags::TimeStep>>(
429 : [](const gsl::not_null<::TimeDelta*> time_step,
430 : const gsl::not_null<::TimeDelta*> next_time_step,
431 : const std::tuple<::TimeDelta>& initial_step,
432 : const std::tuple<::TimeDelta>& initial_next_step) {
433 : *time_step = get<0>(initial_step);
434 : *next_time_step = get<0>(initial_next_step);
435 : },
436 : make_not_null(&box), db::get<initial_step_tag>(box),
437 : db::get<Tags::InitialValue<::Tags::Next<::Tags::TimeStep>>>(box));
438 : using remove_tags = tmpl::filter<DbTags, is_a_initial_value<tmpl::_1>>;
439 : // reset each tag to default constructed values to reduce memory usage (Data
440 : // structures like `DataVector`s and `Tensor`s have negligible memory usage
441 : // when default-constructed). This tends to be better for build time than
442 : // constructing more DataBox types with and without this handful of tags.
443 : tmpl::for_each<remove_tags>([&box](auto tag_v) {
444 : using tag = typename decltype(tag_v)::type;
445 : db::mutate<tag>(
446 : [](auto tag_value) {
447 : *tag_value = std::decay_t<decltype(*tag_value)>{};
448 : },
449 : make_not_null(&box));
450 : });
451 : using history_tags = ::Tags::get_all_history_tags<DbTags>;
452 : tmpl::for_each<history_tags>([&box](auto tag_v) {
453 : using tag = typename decltype(tag_v)::type;
454 : ASSERT(db::get<tag>(box).integration_order() ==
455 : db::get<::Tags::TimeStepper<TimeStepper>>(box).order(),
456 : "Volume history order is: "
457 : << db::get<tag>(box).integration_order()
458 : << " but time stepper requires order: "
459 : << db::get<::Tags::TimeStepper<TimeStepper>>(box).order()
460 : << ". This may indicate that the step size has varied during "
461 : "self-start, which should not be permitted.");
462 : });
463 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
464 : }
465 : };
466 : } // namespace Actions
467 :
468 : namespace detail {
469 : struct PhaseStart;
470 : struct PhaseEnd;
471 : } // namespace detail
472 :
473 : /// \ingroup TimeGroup
474 : /// The list of actions required to self-start an integrator.
475 : ///
476 : /// \tparam StepActions List of actions computing and recording the
477 : /// system derivative and updating the evolved variables (but not the
478 : /// time).
479 : ///
480 : /// \see SelfStart
481 : template <typename StepActions, typename System>
482 1 : using self_start_procedure = tmpl::flatten<tmpl::list<
483 : // clang-format off
484 : SelfStart::Actions::Initialize<System>,
485 : ::Actions::Label<detail::PhaseStart>,
486 : SelfStart::Actions::CheckForCompletion<detail::PhaseEnd, System>,
487 : ::Actions::AdvanceTime,
488 : SelfStart::Actions::CheckForOrderIncrease,
489 : StepActions,
490 : ::Actions::Goto<detail::PhaseStart>,
491 : ::Actions::Label<detail::PhaseEnd>,
492 : SelfStart::Actions::Cleanup,
493 : ::Actions::AdvanceTime,
494 : Parallel::Actions::TerminatePhase>>;
495 : // clang-format on
496 : } // namespace SelfStart
|