SelfStartActions.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <string>
8 #include <tuple>
9 
13 // IWYU pragma: no_include "DataStructures/Tensor/Tensor.hpp"
14 #include "Parallel/GotoAction.hpp" // IWYU pragma: keep
15 #include "Time/Actions/AdvanceTime.hpp" // IWYU pragma: keep
16 #include "Time/Slab.hpp"
17 #include "Time/Tags.hpp" // IWYU pragma: keep // for item_type<Tags::TimeStep>
18 #include "Time/Time.hpp"
19 #include "Time/TimeId.hpp"
20 #include "Utilities/Gsl.hpp"
21 #include "Utilities/TMPL.hpp"
23 #include "Utilities/TypeTraits.hpp"
24 
25 /// \cond
26 namespace Parallel {
27 template <typename Metavariables>
28 class ConstGlobalCache;
29 } // namespace Parallel
30 // IWYU pragma: no_forward_declare db::DataBox
31 /// \endcond
32 
33 /// \ingroup TimeGroup
34 /// Definition of the integrator self-starting procedure.
35 ///
36 /// The self-start procedure generates \f$N\f$ function values of
37 /// accuracy \f$O\left((\Delta t)^N\right)\f$, where \f$N\f$ is the
38 /// negative of the initial slab number. To generate values, it
39 /// requires the time stepper to be a multistep integrator that will
40 /// produce an order-\f$k\f$-accurate result given \f$k-1\f$ history
41 /// values.
42 ///
43 /// If the integrator is started from analytic history data or
44 /// requires no history (such as for a substep integrator), then the
45 /// initial slab number can be set to zero and no self-start steps
46 /// will be taken.
47 ///
48 /// \details
49 /// To self-start a multistep integrator, the function is integrated
50 /// repeatedly with increasing accuracy. A first order integrator
51 /// (Euler's method) requires no history values, so it can be used,
52 /// starting from the initial point, to generate a
53 /// first-order-accurate value at a later time. We then reset to the
54 /// start conditions and use the new "history" value (at a discarded
55 /// point later in time) to take two steps with a second-order method.
56 /// These values are second-order-accurate despite the history only
57 /// being first-order because the calculation of the change in the
58 /// value multiplies the previous derivatives by a factor of
59 /// \f$\Delta t\f$. The time and value are then again reset to their
60 /// starting values and we start again at third order, and so on.
61 ///
62 /// The choice of performing the low-order integrations in the same
63 /// direction as the main integration makes this a _forward_
64 /// self-start procedure, as opposed to a _backward_ procedure that
65 /// produces values for times before the start time. The primary
66 /// advantage of the forward version is that the solution is
67 /// guaranteed to exist after the start time, but not before. It also
68 /// makes bookkeeping easier, as the reset after each order increase
69 /// is to the initial state, rather than to a time one step further
70 /// back for each order. It does have the disadvantage, however, of
71 /// leaving a non-monotonic history at the end of the procedure, which
72 /// the main evolution loop must be able to handle.
73 ///
74 /// Each time the state is reset the slab number is increased by one.
75 /// This ensures that the evaluations are considered to be ordered in
76 /// their evaluation order, even though they are not monotonic in
77 /// time. When the slab number reaches zero the initialization
78 /// procedure is complete and history appropriate for use for an
79 /// integrator of order \f$N+1\f$ has been generated.
80 ///
81 /// The self-start procedure performs all its evaluations before the
82 /// end of the first time step of the main evolution. It is important
83 /// that none of the early steps fall at the same time as the
84 /// self-start history values, so the main evolution should not
85 /// decrease its step size on the first step after the procedure.
86 /// Additionally, the history times will not be monotonicly increasing
87 /// until \f$N\f$ steps have been taken. The local-time-stepping
88 /// calculations require monotonic time, so local time-stepping should
89 /// not be initiated until the self-start values have expired from the
90 /// history. These restrictions on step-size changing are checked in
91 /// the TimeStepper::can_change_step_size method.
92 namespace SelfStart {
93 /// Self-start tags
94 namespace Tags {
95 /// \ingroup TimeGroup
96 /// The initial value of a quantity. The contents are stored in a
97 /// tuple to avoid putting duplicate tensors into the DataBox.
98 template <typename Tag>
100  static std::string name() noexcept {
101  return "InitialValue(" + Tag::name() + ")";
102  }
103  using tag = Tag;
105 };
106 } // namespace Tags
107 
108 /// Self-start actions
109 namespace Actions {
110 namespace detail {
111 template <typename System, bool HasPrimitives>
112 struct vars_to_save_impl {
113  using type = tmpl::list<typename System::variables_tag>;
114 };
115 
116 template <typename System>
117 struct vars_to_save_impl<System, true> {
118  using type = tmpl::list<typename System::variables_tag,
119  typename System::primitive_variables_tag>;
120 };
121 
122 template <typename System>
123 using vars_to_save = typename vars_to_save_impl<
124  System, System::has_primitive_and_conservative_vars>::type;
125 } // namespace detail
126 
127 /// \ingroup ActionsGroup
128 /// \ingroup TimeGroup
129 /// Prepares the evolution for time-stepper self-starting.
130 ///
131 /// Stores the initial values of the variables and time step and sets
132 /// an appropriate step for self-starting.
133 ///
134 /// \details The self-start procedure must take place within one slab,
135 /// and we want to avoid the end of the slab so that we don't have to
136 /// deal with another action advancing the slab on us. There will be
137 /// problems if the main evolution tries to evaluate at a time that
138 /// the self-start procedure used before the self-start version falls
139 /// out of the history, so we have to make sure that does not happen.
140 /// We can't do that by making sure our steps are large enough to keep
141 /// ahead because of the slab restriction, so instead we have to make
142 /// the self-start step smaller to ensure no collisions. The easiest
143 /// way to do that is to fit the entire procedure before the first
144 /// real step, so we pick an initialization time step of
145 /// \f$\Delta t/(N+1)\f$, for \f$\Delta t\f$ the initial evolution
146 /// time step and \f$N\f$ the number of history points to be
147 /// generated.
148 ///
149 /// Uses:
150 /// - ConstGlobalCache: nothing
151 /// - DataBox:
152 /// - Tags::TimeId
153 /// - Tags::TimeStep
154 /// - variables_tag
155 /// - primitive_variables_tag if the system has primitives
156 ///
157 /// DataBox changes:
158 /// - Adds:
159 /// - SelfStart::Tags::InitialValue<Tags::TimeStep>
160 /// - SelfStart::Tags::InitialValue<variables_tag>
161 /// - SelfStart::Tags::InitialValue<primitive_variables_tag> if the system
162 /// has primitives
163 /// - Removes: nothing
164 /// - Modifies: Tags::TimeStep
165 struct Initialize {
166  private:
167  template <typename TagsToSave>
168  struct StoreInitialValues;
169 
170  template <typename... TagsToSave>
171  struct StoreInitialValues<tmpl::list<TagsToSave...>> {
172  template <typename Box>
173  static auto apply(Box box) noexcept {
174  return db::create_from<
177  std::move(box), std::make_tuple(db::get<TagsToSave>(box))...);
178  }
179  };
180 
181  public:
182  template <typename DbTags, typename... InboxTags, typename Metavariables,
183  typename ArrayIndex, typename ActionList,
184  typename ParallelComponent>
185  static auto apply(db::DataBox<DbTags>& box,
186  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
188  const ArrayIndex& /*array_index*/,
189  const ActionList /*meta*/,
190  const ParallelComponent* const /*meta*/) noexcept {
191  using system = typename Metavariables::system;
192 
193  const TimeDelta initial_step = db::get<::Tags::TimeStep>(box);
194  // The slab number increments each time a new point is generated
195  // until it reaches zero.
196  const auto values_needed =
197  -db::get<::Tags::Next<::Tags::TimeId>>(box).slab_number();
198  const TimeDelta self_start_step = initial_step / (values_needed + 1);
199 
200  auto new_box = StoreInitialValues<tmpl::push_back<
201  detail::vars_to_save<system>, ::Tags::TimeStep>>::apply(std::move(box));
202 
203  db::mutate<::Tags::TimeStep>(
204  make_not_null(&new_box), [&self_start_step](
205  const gsl::not_null<
207  time_step) noexcept {
208  *time_step = self_start_step;
209  });
210 
211  return std::make_tuple(std::move(new_box));
212  }
213 };
214 
215 /// \ingroup ActionsGroup
216 /// \ingroup TimeGroup
217 /// Terminates the self-start phase if the required order has been
218 /// reached.
219 ///
220 /// Uses:
221 /// - ConstGlobalCache: nothing
222 /// - DataBox: Tags::Next<Tags::TimeId>
223 ///
224 /// DataBox changes:
225 /// - Adds: nothing
226 /// - Removes: nothing
227 /// - Modifies: nothing
228 template <typename ExitTag>
230  template <typename DbTags, typename... InboxTags, typename Metavariables,
231  typename ArrayIndex, typename ActionList,
232  typename ParallelComponent>
233  static auto apply(db::DataBox<DbTags>& box,
234  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
236  const ArrayIndex& /*array_index*/,
237  const ActionList /*meta*/,
238  const ParallelComponent* const /*meta*/) noexcept {
239  // The self start procedure begins with slab number
240  // -number_of_past_steps and counts up. When we reach 0 we should
241  // start the evolution proper. The first thing the evolution loop
242  // will do is update the time, so here we need to check if the
243  // next time should be the first real step.
244  return std::tuple<db::DataBox<DbTags>&&, bool, size_t>(
245  std::move(box), false,
246  db::get<::Tags::Next<::Tags::TimeId>>(box).slab_number() == 0
247  ? tmpl::index_of<ActionList, ::Actions::Label<ExitTag>>::value
248  : tmpl::index_of<ActionList, CheckForCompletion>::value + 1);
249  // Once we have full support for phases this action should
250  // terminate the phase:
251  // return std::tuple<db::DataBox<DbTags>&&, bool>(
252  // std::move(box),
253  // db::get<::Tags::Next<::Tags::TimeId>>(box).slab_number() == 0);
254  }
255 };
256 
257 /// \ingroup ActionsGroup
258 /// \ingroup TimeGroup
259 /// If we have taken enough steps for this order, set the next time to
260 /// the start time and increment the slab number
261 ///
262 /// Uses:
263 /// - ConstGlobalCache: nothing
264 /// - DataBox:
265 /// - Tags::HistoryEvolvedVariables<variables_tag, dt_variables_tag>
266 /// - Tags::Time
267 /// - Tags::TimeId
268 /// - Tags::TimeStep
269 ///
270 /// DataBox changes:
271 /// - Adds: nothing
272 /// - Removes: nothing
273 /// - Modifies: Tags::Next<Tags::TimeId> if there is an order increase
275  template <typename DbTags, typename... InboxTags, typename Metavariables,
276  typename ArrayIndex, typename ActionList,
277  typename ParallelComponent>
278  static auto apply(db::DataBox<DbTags>& box,
279  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
281  const ArrayIndex& /*array_index*/,
282  const ActionList /*meta*/,
283  const ParallelComponent* const /*meta*/) noexcept {
284  using variables_tag = typename Metavariables::system::variables_tag;
285 
286  const auto& time = db::get<::Tags::Time>(box);
287  const auto& time_step = db::get<::Tags::TimeStep>(box);
288  const auto& history = db::get<::Tags::HistoryEvolvedVariables<
289  variables_tag, db::add_tag_prefix<::Tags::dt, variables_tag>>>(box);
290 
291  const Time required_time =
292  (time_step.is_positive() ? time.slab().start() : time.slab().end()) +
293  (history.size() + 1) * time_step;
294  const bool done_with_order = time == required_time;
295 
296  if (done_with_order) {
297  db::mutate<::Tags::Next<::Tags::TimeId>>(
298  make_not_null(&box),
300  next_time_id,
301  const db::item_type<::Tags::TimeId>& current_time_id) noexcept {
302  const Slab slab = current_time_id.time().slab();
303  *next_time_id =
304  TimeId(current_time_id.time_runs_forward(),
305  current_time_id.slab_number() + 1,
306  current_time_id.time_runs_forward() ? slab.start()
307  : slab.end());
308  },
309  db::get<::Tags::TimeId>(box));
310  }
311 
312  return std::forward_as_tuple(std::move(box));
313  }
314 };
315 
316 /// \ingroup ActionsGroup
317 /// \ingroup TimeGroup
318 /// Jumps to the start of the self-start algorithm (skipping taking a
319 /// step from the last point) if the generation of the points for the
320 /// current order is complete.
321 ///
322 /// Uses:
323 /// - ConstGlobalCache: nothing
324 /// - DataBox:
325 /// - Tags::Next<Tags::TimeId>
326 /// - SelfStart::Tags::InitialValue<variables_tag>
327 /// - SelfStart::Tags::InitialValue<primitive_variables_tag> if the system
328 /// has primitives
329 ///
330 /// DataBox changes:
331 /// - Adds: nothing
332 /// - Removes: nothing
333 /// - Modifies:
334 /// - variables_tag if there is an order increase
335 /// - primitive_variables_tag if there is an order increase and the system
336 /// has primitives
337 template <typename RestartTag>
339  template <typename DbTags, typename... InboxTags, typename Metavariables,
340  typename ArrayIndex, typename ActionList,
341  typename ParallelComponent>
342  static auto apply(db::DataBox<DbTags>& box,
343  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
345  const ArrayIndex& /*array_index*/,
346  const ActionList /*meta*/,
347  const ParallelComponent* const /*meta*/) noexcept {
348  using system = typename Metavariables::system;
349 
350  constexpr size_t restart_index =
351  tmpl::index_of<ActionList, ::Actions::Label<RestartTag>>::value + 1;
352  constexpr size_t continue_index =
353  tmpl::index_of<ActionList, StartNextOrderIfReady>::value + 1;
354 
355  const bool done_with_order =
356  db::get<::Tags::Next<::Tags::TimeId>>(box).is_at_slab_boundary();
357 
358  if (done_with_order) {
359  tmpl::for_each<detail::vars_to_save<system>>([&box](auto tag) noexcept {
360  using Tag = tmpl::type_from<decltype(tag)>;
361  db::mutate<Tag>(
362  make_not_null(&box),
363  [](const gsl::not_null<db::item_type<Tag>*> value,
365  initial_value) noexcept {
366  *value = get<0>(initial_value);
367  },
368  db::get<Tags::InitialValue<Tag>>(box));
369  });
370  }
371 
372  return std::tuple<db::DataBox<DbTags>&&, bool, size_t>(
373  std::move(box), false,
374  done_with_order ? restart_index : continue_index);
375  }
376 };
377 
378 /// \ingroup ActionsGroup
379 /// \ingroup TimeGroup
380 /// Cleans up after the self-start procedure
381 ///
382 /// Resets the time step to that requested for the evolution and
383 /// removes temporary self-start data.
384 ///
385 /// Uses:
386 /// - ConstGlobalCache: nothing
387 /// - DataBox: SelfStart::Tags::InitialValue<Tags::TimeStep>
388 ///
389 /// DataBox changes:
390 /// - Adds: nothing
391 /// - Removes:
392 /// - All SelfStart::Tags::InitialValue tags
393 /// - Modifies: Tags::TimeStep
394 struct Cleanup {
395  private:
396  template <typename T>
397  struct is_a_initial_value : tt::is_a<Tags::InitialValue, T> {};
398 
399  public:
400  template <typename DbTags, typename... InboxTags, typename Metavariables,
401  typename ArrayIndex, typename ActionList,
402  typename ParallelComponent>
403  static auto apply(db::DataBox<DbTags>& box,
404  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
406  const ArrayIndex& /*array_index*/,
407  const ActionList /*meta*/,
408  const ParallelComponent* const /*meta*/) noexcept {
409  using initial_step_tag = Tags::InitialValue<::Tags::TimeStep>;
410 
411  // Reset the time step to the value requested by the user. The
412  // variables were reset in StartNextOrderIfReady.
413  db::mutate<::Tags::TimeStep>(
414  make_not_null(&box),
416  const db::item_type<initial_step_tag>& initial_step) noexcept {
417  *time_step = get<0>(initial_step);
418  },
419  db::get<initial_step_tag>(box));
420 
421  using remove_tags = tmpl::filter<DbTags, is_a_initial_value<tmpl::_1>>;
422  return std::make_tuple(db::create_from<remove_tags>(std::move(box)));
423  }
424 };
425 } // namespace Actions
426 
427 namespace detail {
428 struct PhaseStart;
429 struct PhaseEnd;
430 } // namespace detail
431 
432 /// \ingroup TimeGroup
433 /// The list of actions required to self-start an integrator.
434 ///
435 /// \tparam ComputeRhs Action or list of actions computing and
436 /// recording the system derivative.
437 /// \tparam UpdateVariables Action or list of actions updating the
438 /// evolved variables (but not the time).
439 ///
440 /// \see SelfStart
441 // clang-format off
442 template <typename ComputeRhs, typename UpdateVariables>
443 using self_start_procedure = tmpl::flatten<tmpl::list<
449  ComputeRhs,
451  UpdateVariables,
455 // clang-format on
456 } // namespace SelfStart
Prefix< DataBox_detail::dispatch_add_tag_prefix_impl< Prefix, Tag, Args... >, Args... > add_tag_prefix
Wrap Tag in Prefix<_, Args...>, also wrapping variables tags if Tag is a Tags::Variables.
Definition: DataBoxTag.hpp:533
Tag for step size.
Definition: Tags.hpp:45
Defines class tuples::TaggedTuple.
Terminates the self-start phase if the required order has been reached.
Definition: SelfStartActions.hpp:229
Cleans up after the self-start procedure.
Definition: SelfStartActions.hpp:394
Defines action AdvanceTime.
The time in a simulation. Times can be safely compared for exact equality as long as they do not belo...
Definition: Time.hpp:31
Labels a location in the action list that can be jumped to using Goto.
Definition: GotoAction.hpp:37
A unique identifier for the temporal state of an integrated system.
Definition: TimeId.hpp:25
Define prefixes for DataBox tags.
Tags for the DataBox inherit from this type.
Definition: DataBoxTag.hpp:65
A chunk of time. Every element must reach slab boundaries exactly, no matter how it actually takes ti...
Definition: Slab.hpp:29
Defines class TimeId.
Contains functions that forward to Charm++ parallel functions.
Definition: Abort.hpp:13
tmpl::flatten< tmpl::list< Tags... > > RemoveTags
List of Tags to remove from the DataBox.
Definition: DataBox.hpp:1220
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args)
Apply the function f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1595
Defines Time and TimeDelta.
Definition: Determinant.hpp:11
Prepares the evolution for time-stepper self-starting.
Definition: SelfStartActions.hpp:165
Prefix for TimeStepper history.
Definition: Tags.hpp:75
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Defines classes and functions used for manipulating DataBox&#39;s.
Jumps to the start of the self-start algorithm (skipping taking a step from the last point) if the ge...
Definition: SelfStartActions.hpp:338
constexpr auto create_from(db::DataBox< TagsList > &&box, Args &&... args) noexcept
Create a new DataBox from an existing one adding or removing items and compute items.
Definition: DataBox.hpp:1414
Jumps to a Label.
Definition: GotoAction.hpp:68
tmpl::flatten< tmpl::list< Tags... > > AddSimpleTags
List of Tags to add to the DataBox.
Definition: DataBox.hpp:1227
Represents an interval of time within a single slab.
Definition: Time.hpp:108
Defines class Slab.
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
Definition: DataBoxTag.hpp:29
A Charm++ chare that caches constant data once per Charm++ node.
Definition: ConstGlobalCache.hpp:76
Check if type T is a template specialization of U
Definition: TypeTraits.hpp:536
If we have taken enough steps for this order, set the next time to the start time and increment the s...
Definition: SelfStartActions.hpp:274
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:1211
Definition of the integrator self-starting procedure.
Definition: SelfStartActions.hpp:92
tmpl::flatten< tmpl::list< SelfStart::Actions::Initialize, ::Actions::Label< detail::PhaseStart >, SelfStart::Actions::CheckForCompletion< detail::PhaseEnd >, ::Actions::AdvanceTime, SelfStart::Actions::CheckForOrderIncrease, ComputeRhs, SelfStart::Actions::StartNextOrderIfReady< detail::PhaseStart >, UpdateVariables, ::Actions::Goto< detail::PhaseStart >, ::Actions::Label< detail::PhaseEnd >, SelfStart::Actions::Cleanup > > self_start_procedure
The list of actions required to self-start an integrator.
Definition: SelfStartActions.hpp:454
The Poisson equation formulated as a set of coupled first-order PDEs.
Definition: FirstOrderSystem.hpp:55
Wraps the template metaprogramming library used (brigand)
Advance time one substep.
Definition: AdvanceTime.hpp:45
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that is returned by the Tag. If it is a base tag then a TagList must be passed as a seco...
Definition: DataBoxTag.hpp:410
Defines functions and classes from the GSL.
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:863
Marks an item as being a prefix to another tag.
Definition: DataBoxTag.hpp:112
Definition: SolvePoissonProblem.hpp:38
The initial value of a quantity. The contents are stored in a tuple to avoid putting duplicate tensor...
Definition: SelfStartActions.hpp:99
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
Defines type traits, some of which are future STL type_traits header.
Defines tags related to Time quantities.
Definition: ComputeTimeDerivative.hpp:28
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12