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