Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <boost/container/static_vector.hpp>
7 : #include <cstddef>
8 : #include <optional>
9 : #include <pup.h>
10 : #include <type_traits>
11 : #include <utility>
12 :
13 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
14 : #include "DataStructures/DataBox/Prefixes.hpp"
15 : #include "DataStructures/MathWrapper.hpp"
16 : #include "DataStructures/StaticDeque.hpp"
17 : #include "Time/EvolutionOrdering.hpp"
18 : #include "Time/TimeStepId.hpp"
19 : #include "Utilities/ContainsAllocations.hpp"
20 : #include "Utilities/ErrorHandling/Assert.hpp"
21 : #include "Utilities/ErrorHandling/Error.hpp"
22 : #include "Utilities/Gsl.hpp"
23 : #include "Utilities/Serialization/PupBoost.hpp"
24 : #include "Utilities/Serialization/PupStlCpp17.hpp"
25 : #include "Utilities/StdHelpers.hpp"
26 : #include "Utilities/StlBoilerplate.hpp"
27 : #include "Utilities/TMPL.hpp"
28 :
29 : namespace TimeSteppers {
30 :
31 : /// Largest number of past steps supported by the time stepper
32 : /// `History`. Corresponds to the `number_of_past_steps()` method of
33 : /// `TimeStepper`. AdamsBashforth with order 8 has the largest
34 : /// requirement.
35 1 : constexpr size_t history_max_past_steps = 7;
36 : /// Largest number of substeps supported by the time stepper
37 : /// `History`. Corresponds to the `number_of_substeps()` and
38 : /// `number_of_substeps_for_error()` methods of `TimeStepper`.
39 : /// DormandPrince5 with an error estimate and Rk5Owren have the
40 : /// largest requirement.
41 1 : constexpr size_t history_max_substeps = 6;
42 :
43 : /// \ingroup TimeSteppersGroup
44 : /// Entry in the time-stepper history, in type-erased form.
45 : ///
46 : /// The history access classes do not provide mutable references to
47 : /// these structs, so they cannot be used to modify history data.
48 : ///
49 : /// This struct mirrors the typed `StepRecord` struct. See that for
50 : /// details.
51 : ///
52 : /// \tparam T One of the types in \ref MATH_WRAPPER_TYPES
53 : template <typename T>
54 1 : struct UntypedStepRecord {
55 0 : TimeStepId time_step_id;
56 0 : std::optional<T> value;
57 0 : T derivative;
58 :
59 : static_assert(tmpl::list_contains_v<tmpl::list<MATH_WRAPPER_TYPES>, T>);
60 : };
61 :
62 : template <typename T>
63 0 : bool operator==(const UntypedStepRecord<T>& a, const UntypedStepRecord<T>& b);
64 :
65 : template <typename T>
66 0 : bool operator!=(const UntypedStepRecord<T>& a, const UntypedStepRecord<T>& b);
67 :
68 : #if defined(__GNUC__) and not defined(__clang__)
69 : #pragma GCC diagnostic push
70 : #pragma GCC diagnostic ignored "-Wnon-virtual-dtor"
71 : #endif // defined(__GNUC__) and not defined(__clang__)
72 : /// \ingroup TimeSteppersGroup
73 : /// Access to the history data used by a TimeStepper in type-erased
74 : /// form. Obtain an instance with `History::untyped()`.
75 : ///
76 : /// The methods mirror similar ones in `History`. See that class for
77 : /// details.
78 : ///
79 : /// \tparam T One of the types in \ref MATH_WRAPPER_TYPES
80 : template <typename T>
81 1 : class ConstUntypedHistory
82 : : public stl_boilerplate::RandomAccessSequence<
83 : ConstUntypedHistory<T>, const UntypedStepRecord<T>, false> {
84 : static_assert(tmpl::list_contains_v<tmpl::list<MATH_WRAPPER_TYPES>, T>);
85 :
86 : protected:
87 0 : ConstUntypedHistory() = default;
88 0 : ConstUntypedHistory(const ConstUntypedHistory&) = default;
89 0 : ConstUntypedHistory(ConstUntypedHistory&&) = default;
90 0 : ConstUntypedHistory& operator=(const ConstUntypedHistory&) = default;
91 0 : ConstUntypedHistory& operator=(ConstUntypedHistory&&) = default;
92 0 : ~ConstUntypedHistory() = default;
93 :
94 : public:
95 0 : using WrapperType = T;
96 :
97 : /// \cond
98 : class UntypedSubsteps
99 : : public stl_boilerplate::RandomAccessSequence<
100 : UntypedSubsteps, const UntypedStepRecord<T>, true> {
101 : public:
102 : UntypedSubsteps() = delete;
103 : UntypedSubsteps(const UntypedSubsteps&) = delete;
104 : UntypedSubsteps(UntypedSubsteps&&) = default;
105 : UntypedSubsteps& operator=(const UntypedSubsteps&) = delete;
106 : UntypedSubsteps& operator=(UntypedSubsteps&&) = default;
107 : ~UntypedSubsteps() = default;
108 :
109 : using WrapperType = T;
110 :
111 : size_t size() const;
112 : size_t max_size() const;
113 :
114 : const UntypedStepRecord<T>& operator[](const size_t index) const;
115 :
116 : private:
117 : friend ConstUntypedHistory;
118 : explicit UntypedSubsteps(const ConstUntypedHistory& history);
119 :
120 : gsl::not_null<const ConstUntypedHistory*> history_;
121 : };
122 : /// \endcond
123 :
124 0 : virtual size_t integration_order() const = 0;
125 :
126 0 : virtual size_t size() const = 0;
127 0 : virtual size_t max_size() const = 0;
128 :
129 0 : virtual const UntypedStepRecord<T>& operator[](size_t index) const = 0;
130 :
131 0 : virtual const UntypedStepRecord<T>& operator[](
132 : const TimeStepId& id) const = 0;
133 :
134 0 : virtual bool at_step_start() const = 0;
135 :
136 0 : UntypedSubsteps substeps() const;
137 :
138 : private:
139 0 : friend UntypedSubsteps;
140 : virtual const boost::container::static_vector<UntypedStepRecord<WrapperType>,
141 : history_max_substeps>&
142 0 : substep_values() const = 0;
143 : };
144 :
145 : /// \ingroup TimeSteppersGroup
146 : /// Mutable access to the history data used by a TimeStepper in
147 : /// type-erased form. Obtain an instance with `History::untyped()`.
148 : ///
149 : /// Data cannot be inserted or modified through the type-erased
150 : /// interface. The only mutability exposed is the ability to delete
151 : /// data.
152 : ///
153 : /// The methods mirror similar ones in `History`. See that class for
154 : /// details.
155 : ///
156 : /// \tparam T One of the types in \ref MATH_WRAPPER_TYPES
157 : template <typename T>
158 1 : class MutableUntypedHistory : public ConstUntypedHistory<T> {
159 : static_assert(tmpl::list_contains_v<tmpl::list<MATH_WRAPPER_TYPES>, T>);
160 :
161 : protected:
162 0 : MutableUntypedHistory() = default;
163 0 : MutableUntypedHistory(const MutableUntypedHistory&) = default;
164 0 : MutableUntypedHistory(MutableUntypedHistory&&) = default;
165 0 : MutableUntypedHistory& operator=(const MutableUntypedHistory&) = default;
166 0 : MutableUntypedHistory& operator=(MutableUntypedHistory&&) = default;
167 0 : ~MutableUntypedHistory() = default;
168 :
169 : public:
170 0 : virtual void discard_value(const TimeStepId& id_to_discard) const = 0;
171 :
172 0 : virtual void pop_front() const = 0;
173 :
174 0 : virtual void clear_substeps() const = 0;
175 : };
176 : #if defined(__GNUC__) and not defined(__clang__)
177 : #pragma GCC diagnostic pop
178 : #endif // defined(__GNUC__) and not defined(__clang__)
179 :
180 : /// \ingroup TimeSteppersGroup
181 : /// Data in an entry of the time-stepper history.
182 : ///
183 : /// The `value` field may be empty if the value has been discarded to
184 : /// save memory. See `History::discard_value` and
185 : /// `ConstUntypedHistory::discard_value`.
186 : template <typename Vars>
187 1 : struct StepRecord {
188 0 : using DerivVars = db::prefix_variables<::Tags::dt, Vars>;
189 :
190 0 : TimeStepId time_step_id;
191 0 : std::optional<Vars> value;
192 0 : DerivVars derivative;
193 :
194 : // NOLINTNEXTLINE(google-runtime-references)
195 0 : void pup(PUP::er& p);
196 :
197 0 : std::ostream& print(std::ostream& os) const;
198 : };
199 :
200 : template <typename Vars>
201 : void StepRecord<Vars>::pup(PUP::er& p) {
202 : p | time_step_id;
203 : p | value;
204 : p | derivative;
205 : }
206 :
207 : template <typename Vars>
208 : std::ostream& StepRecord<Vars>::print(std::ostream& os) const {
209 : using ::operator<<;
210 : os << "TimeStepId: " << time_step_id << "\n";
211 : os << "Value: " << value << "\n";
212 : os << "Derivative: " << derivative << "\n";
213 : return os;
214 : }
215 :
216 : template <typename Vars>
217 0 : bool operator==(const StepRecord<Vars>& a, const StepRecord<Vars>& b) {
218 : return a.time_step_id == b.time_step_id and a.value == b.value and
219 : a.derivative == b.derivative;
220 : }
221 :
222 : template <typename Vars>
223 0 : bool operator!=(const StepRecord<Vars>& a, const StepRecord<Vars>& b) {
224 : return not(a == b);
225 : }
226 :
227 : /// \cond
228 : template <typename Vars>
229 : class History;
230 : /// \endcond
231 :
232 : namespace History_detail {
233 : // Find a record from either the steps or substeps. This takes the
234 : // arguments it does so it can be used to find both const and
235 : // non-const values from both typed and untyped histories.
236 : template <typename History>
237 : decltype(auto) find_record(History&& history, const TimeStepId& id) {
238 : const size_t substep = id.substep();
239 : if (substep == 0) {
240 : for (size_t i = 0; i < history.size(); ++i) {
241 : auto& record = history[i];
242 : if (record.time_step_id == id) {
243 : return record;
244 : }
245 : }
246 : ERROR(id << " not present");
247 : } else {
248 : ASSERT(substep - 1 < history.substeps().size(), id << " not present");
249 : auto& record = history.substeps()[substep - 1];
250 : ASSERT(record.time_step_id == id, id << " not present");
251 : return record;
252 : }
253 : }
254 :
255 : #if defined(__GNUC__) and not defined(__clang__)
256 : #pragma GCC diagnostic push
257 : #pragma GCC diagnostic ignored "-Wnon-virtual-dtor"
258 : #endif // defined(__GNUC__) and not defined(__clang__)
259 : template <typename UntypedBase>
260 : class UntypedAccessCommon : public UntypedBase {
261 : protected:
262 : UntypedAccessCommon() = delete;
263 : UntypedAccessCommon(const UntypedAccessCommon&) = delete;
264 : UntypedAccessCommon(UntypedAccessCommon&&) = default;
265 : UntypedAccessCommon& operator=(const UntypedAccessCommon&) = delete;
266 : // Can't move-assign non-owning DataVectors.
267 : UntypedAccessCommon& operator=(UntypedAccessCommon&&) = delete;
268 : ~UntypedAccessCommon() = default;
269 :
270 : public:
271 : using WrapperType = typename UntypedBase::WrapperType;
272 :
273 : size_t integration_order() const override;
274 :
275 : size_t size() const override;
276 : size_t max_size() const override;
277 :
278 : const UntypedStepRecord<WrapperType>& operator[](
279 : const size_t index) const override;
280 : const UntypedStepRecord<WrapperType>& operator[](
281 : const TimeStepId& id) const override;
282 :
283 : bool at_step_start() const override;
284 :
285 : protected:
286 : template <typename Vars>
287 : UntypedAccessCommon(const History<Vars>& history) {
288 : reinitialize(history);
289 : }
290 :
291 : template <typename Vars>
292 : void reinitialize(const History<Vars>& history) const {
293 : // This class is presenting a view. Regenerating its internal
294 : // representation after a change to the parent structure is an
295 : // implementation detail.
296 : auto* const mutable_this = const_cast<UntypedAccessCommon*>(this);
297 : mutable_this->integration_order_ = history.integration_order();
298 : mutable_this->step_values_.clear();
299 : for (const StepRecord<Vars>& record : history) {
300 : mutable_this->step_values_.push_back(make_untyped(record));
301 : }
302 : mutable_this->substep_values_.clear();
303 : for (const StepRecord<Vars>& record : history.substeps()) {
304 : mutable_this->substep_values_.push_back(make_untyped(record));
305 : }
306 : }
307 :
308 : private:
309 : const boost::container::static_vector<UntypedStepRecord<WrapperType>,
310 : history_max_substeps>&
311 : substep_values() const override;
312 :
313 : template <typename Vars>
314 : static UntypedStepRecord<WrapperType> make_untyped(
315 : const StepRecord<Vars>& record) {
316 : // This class only exposes these records as const references, so
317 : // it is OK if we break non-allocating references to the original
318 : // data by moving out of the MathWrappers since you can't modify
319 : // through them.
320 : return {record.time_step_id,
321 : record.value.has_value() ? std::optional{const_cast<WrapperType&&>(
322 : *make_math_wrapper(*record.value))}
323 : : std::nullopt,
324 : const_cast<WrapperType&&>(*make_math_wrapper(record.derivative))};
325 : }
326 :
327 : size_t integration_order_{};
328 : // static_vector never reallocates, so storing non-owning
329 : // DataVectors is safe.
330 : boost::container::static_vector<UntypedStepRecord<WrapperType>,
331 : history_max_past_steps + 2>
332 : step_values_{};
333 : boost::container::static_vector<UntypedStepRecord<WrapperType>,
334 : history_max_substeps>
335 : substep_values_{};
336 : };
337 : #if defined(__GNUC__) and not defined(__clang__)
338 : #pragma GCC diagnostic pop
339 : #endif // defined(__GNUC__) and not defined(__clang__)
340 : } // namespace History_detail
341 :
342 : /// \ingroup TimeSteppersGroup
343 : /// The past-time data used by TimeStepper classes to update the
344 : /// evolved variables.
345 : ///
346 : /// This class exposes an STL-like container interface for accessing
347 : /// the step (not substep) `StepRecord` data. The records can be
348 : /// freely modified through this interface, although modifying the
349 : /// `time_step_id` field is generally inadvisable.
350 : ///
351 : /// This class is designed to minimize the number of memory
352 : /// allocations performed during a step, and so caches discarded
353 : /// entries for reuse if they contain dynamic allocations. During
354 : /// steady-state operation, this class will perform no heap
355 : /// allocations. If `Vars` and the associated `DerivVars` do not
356 : /// allocate internally, then this class will perform no heap
357 : /// allocations under any circumstances.
358 : template <typename Vars>
359 1 : class History
360 : : public stl_boilerplate::RandomAccessSequence<History<Vars>,
361 : StepRecord<Vars>, false> {
362 : public:
363 0 : using DerivVars = db::prefix_variables<::Tags::dt, Vars>;
364 :
365 0 : History() = default;
366 0 : explicit History(const size_t integration_order)
367 : : integration_order_(integration_order) {}
368 0 : History(const History& other);
369 0 : History(History&&) = default;
370 0 : History& operator=(const History& other);
371 0 : History& operator=(History&&) = default;
372 :
373 : /// The wrapped type presented by the type-erased history. One of
374 : /// the types in \ref MATH_WRAPPER_TYPES.
375 1 : using UntypedVars = math_wrapper_type<Vars>;
376 :
377 : /// \cond
378 : // This wrapper around UntypedAccessCommon exists because we want
379 : // the special members of that to be protected, and we want this to
380 : // be final.
381 : class ConstUntypedAccess final : public History_detail::UntypedAccessCommon<
382 : ConstUntypedHistory<UntypedVars>> {
383 : public:
384 : ConstUntypedAccess() = delete;
385 : ConstUntypedAccess(const ConstUntypedAccess&) = delete;
386 : ConstUntypedAccess(ConstUntypedAccess&&) = default;
387 : ConstUntypedAccess& operator=(const ConstUntypedAccess&) = delete;
388 : ConstUntypedAccess& operator=(ConstUntypedAccess&&) = default;
389 : ~ConstUntypedAccess() = default;
390 :
391 : private:
392 : friend History;
393 :
394 : explicit ConstUntypedAccess(const History& history)
395 : : History_detail::UntypedAccessCommon<ConstUntypedHistory<UntypedVars>>(
396 : history) {}
397 : };
398 :
399 : class MutableUntypedAccess final : public History_detail::UntypedAccessCommon<
400 : MutableUntypedHistory<UntypedVars>> {
401 : public:
402 : MutableUntypedAccess() = delete;
403 : MutableUntypedAccess(const MutableUntypedAccess&) = delete;
404 : MutableUntypedAccess(MutableUntypedAccess&&) = default;
405 : MutableUntypedAccess& operator=(const MutableUntypedAccess&) = delete;
406 : MutableUntypedAccess& operator=(MutableUntypedAccess&&) = default;
407 : ~MutableUntypedAccess() = default;
408 :
409 : void discard_value(const TimeStepId& id_to_discard) const override {
410 : history_->discard_value(id_to_discard);
411 : this->reinitialize(*history_);
412 : }
413 :
414 : void pop_front() const override {
415 : history_->pop_front();
416 : this->reinitialize(*history_);
417 : }
418 :
419 : void clear_substeps() const override {
420 : history_->clear_substeps();
421 : this->reinitialize(*history_);
422 : }
423 :
424 : private:
425 : friend History;
426 :
427 : explicit MutableUntypedAccess(const gsl::not_null<History*> history)
428 : : History_detail::UntypedAccessCommon<
429 : MutableUntypedHistory<UntypedVars>>(*history),
430 : history_(history) {}
431 :
432 : private:
433 : gsl::not_null<History*> history_;
434 : };
435 :
436 : class ConstSubsteps : public stl_boilerplate::RandomAccessSequence<
437 : ConstSubsteps, const StepRecord<Vars>, true> {
438 : public:
439 : ConstSubsteps() = delete;
440 : ConstSubsteps(const ConstSubsteps&) = delete;
441 : ConstSubsteps(ConstSubsteps&&) = default;
442 : ConstSubsteps& operator=(const ConstSubsteps&) = delete;
443 : ConstSubsteps& operator=(ConstSubsteps&&) = default;
444 : ~ConstSubsteps() = default;
445 :
446 : size_t size() const { return history_->substep_values_.size(); }
447 : static constexpr size_t max_size() { return history_max_substeps; }
448 :
449 : const StepRecord<Vars>& operator[](const size_t index) const {
450 : ASSERT(index < size(),
451 : "Requested substep " << index << " but only have " << size());
452 : return history_->substep_values_[index];
453 : }
454 :
455 : private:
456 : friend History;
457 : explicit ConstSubsteps(const History& history) : history_(&history) {}
458 :
459 : gsl::not_null<const History*> history_;
460 : };
461 :
462 : class MutableSubsteps
463 : : public stl_boilerplate::RandomAccessSequence<MutableSubsteps,
464 : StepRecord<Vars>, true> {
465 : public:
466 : MutableSubsteps() = delete;
467 : MutableSubsteps(const MutableSubsteps&) = delete;
468 : MutableSubsteps(MutableSubsteps&&) = default;
469 : MutableSubsteps& operator=(const MutableSubsteps&) = delete;
470 : MutableSubsteps& operator=(MutableSubsteps&&) = default;
471 : ~MutableSubsteps() = default;
472 :
473 : size_t size() const { return history_->substep_values_.size(); }
474 : static constexpr size_t max_size() { return history_max_substeps; }
475 :
476 : StepRecord<Vars>& operator[](const size_t index) const {
477 : ASSERT(index < size(),
478 : "Requested substep " << index << " but only have " << size());
479 : return history_->substep_values_[index];
480 : }
481 :
482 : private:
483 : friend History;
484 : explicit MutableSubsteps(const gsl::not_null<History*> history)
485 : : history_(history) {}
486 :
487 : gsl::not_null<History*> history_;
488 : };
489 : /// \endcond
490 :
491 : /// Immutable, type-erased access to the history. This method
492 : /// returns a class derived from `ConstUntypedHistory<UntypedVars>`.
493 : /// Any modifications to the History class invalidate the object
494 : /// returned by this function.
495 1 : ConstUntypedAccess untyped() const { return ConstUntypedAccess(*this); }
496 :
497 : /// Mutable, type-erased access to the history. This method returns
498 : /// a class derived from `MutableUntypedHistory<UntypedVars>`. Any
499 : /// modifications to the History class invalidate the object
500 : /// returned by this function, except for modifications performed
501 : /// through the object itself.
502 1 : MutableUntypedAccess untyped() { return MutableUntypedAccess(this); }
503 :
504 : /// Get or set the order the time stepper is running at. Many time
505 : /// steppers expect this to have a particular value. This has no
506 : /// effect on the storage of past data.
507 : /// @{
508 1 : size_t integration_order() const { return integration_order_; }
509 1 : void integration_order(const size_t new_integration_order) {
510 : integration_order_ = new_integration_order;
511 : }
512 : /// @}
513 :
514 : /// Type and value used to indicate that a record is to be created
515 : /// without the `value` field set.
516 : /// @{
517 1 : struct NoValue {};
518 0 : static constexpr NoValue no_value{};
519 : /// @}
520 :
521 : /// Insert a new entry into the history. It will be inserted as a
522 : /// step or substep as appropriate.
523 : ///
524 : /// The supplied `time_step_id` must be later than the current
525 : /// latest entry, and if the substep of `time_step_id` is nonzero,
526 : /// the id must be consistent with the existing data for the current
527 : /// step.
528 : ///
529 : /// If the constant `History::no_value` is passed, the created
530 : /// record will not have its `value` field set. This is useful for
531 : /// histories other than the history of the primary evolution
532 : /// integration, such as the implicit portion of an IMEX evolution.
533 : /// These other uses adjust the result of the main history, and so
534 : /// do not require values (which are only used for the zeroth-order
535 : /// terms).
536 : ///
537 : /// The `value` and `derivative` data will be copied into cached
538 : /// allocations, if any are available. That is, only a copy, not a
539 : /// memory allocation, will be performed when possible.
540 : /// @{
541 1 : void insert(const TimeStepId& time_step_id, const Vars& value,
542 : const DerivVars& derivative);
543 1 : void insert(const TimeStepId& time_step_id, NoValue /*unused*/,
544 : const DerivVars& derivative);
545 : /// @}
546 :
547 : /// Insert a new entry in the history by modifying the fields of the
548 : /// record directly, instead of passing a value to copy.
549 : ///
550 : /// The (optional) \p value_inserter must be a functor callable with
551 : /// single argument of type `gsl::not_null<Vars*>` and the \p
552 : /// derivative_inserter must be callable with a single argument of
553 : /// type `gsl::not_null<DerivVars*>`. The passed objects will be
554 : /// created from a cached memory allocation if one is available, and
555 : /// will otherwise be default-constructed.
556 : ///
557 : /// All the restrictions on valid values \p time_step_id for the
558 : /// `insert` method apply here as well.
559 : /// @{
560 : template <typename ValueFunc, typename DerivativeFunc>
561 1 : void insert_in_place(const TimeStepId& time_step_id,
562 : ValueFunc&& value_inserter,
563 : DerivativeFunc&& derivative_inserter);
564 : template <typename DerivativeFunc>
565 1 : void insert_in_place(const TimeStepId& time_step_id, NoValue /*unused*/,
566 : DerivativeFunc&& derivative_inserter);
567 : /// @}
568 :
569 : /// Insert data at the start of the history, similar to `push_front`
570 : /// on some STL containers. This can be useful when initializing a
571 : /// multistep integrator.
572 : ///
573 : /// This should only be used for initialization, and cannot be used
574 : /// for substep data. The supplied `time_step_id` must be earlier
575 : /// than the current first entry.
576 : /// @{
577 1 : void insert_initial(TimeStepId time_step_id, Vars value,
578 : DerivVars derivative);
579 1 : void insert_initial(TimeStepId time_step_id, NoValue /*unused*/,
580 : DerivVars derivative);
581 : /// @}
582 :
583 : /// The number of stored step (not substep) entries.
584 1 : size_t size() const { return step_values_.size(); }
585 : /// The maximum number of step (not substep) entries that can be
586 : /// stored in the history. This number is not very large, but will
587 : /// be sufficient for running a time stepper requiring
588 : /// `history_max_past_steps` past step values.
589 1 : static constexpr size_t max_size() {
590 : // In addition to the past steps, we must store the start and end
591 : // of the current step.
592 : return history_max_past_steps + 2;
593 : }
594 :
595 : /// Access the `StepRecord` for a step (not substep).
596 : /// @{
597 1 : const StepRecord<Vars>& operator[](const size_t index) const {
598 : return step_values_[index];
599 : }
600 1 : StepRecord<Vars>& operator[](const size_t index) {
601 : return step_values_[index];
602 : }
603 : /// @}
604 :
605 : /// Access the `StepRecord` for a step or substep with a given
606 : /// TimeStepId. It is an error if there is no entry with that id.
607 : /// @{
608 1 : const StepRecord<Vars>& operator[](const TimeStepId& id) const {
609 : return History_detail::find_record(*this, id);
610 : }
611 1 : StepRecord<Vars>& operator[](const TimeStepId& id) {
612 : return History_detail::find_record(*this, id);
613 : }
614 : /// @}
615 :
616 : /// Get the value at the latest step or substep in the history, even
617 : /// if the value in that record has been discarded. This is not
618 : /// available if `undo_latest` has been called since the last
619 : /// insertion.
620 : ///
621 : /// This function exists to allow access to the previous value after
622 : /// a potentially bad step has been taken without restricting how
623 : /// the time steppers can manage their history.
624 1 : const Vars& latest_value() const;
625 :
626 : /// Get the record for the start of the step containing the passed
627 : /// time. If the time matches a step time exactly, the record for
628 : /// that step is returned.
629 1 : const StepRecord<Vars>& step_start(double time) const;
630 :
631 : /// Check whether we are at the start of a step, i.e, the most
632 : /// recent entry in the history is not a substep.
633 1 : bool at_step_start() const;
634 :
635 : /// Container view of the `StepRecord`s for the history substeps.
636 : /// These methods return classes providing an STL-like container
637 : /// interface to the substep data. These containers are
638 : /// zero-indexed, so `substeps()[0]` will have a substep value of 1.
639 : ///
640 : /// These containers do not have methods to modify their sizes. Use
641 : /// the methods on the `History` class for those operations.
642 : /// @{
643 1 : ConstSubsteps substeps() const { return ConstSubsteps(*this); }
644 1 : MutableSubsteps substeps() { return MutableSubsteps(this); }
645 : /// @}
646 :
647 : /// Clear the `value` in the indicated record. It is an error if
648 : /// there is no record with the passed TimeStepId. Any memory
649 : /// allocations will be cached for future reuse.
650 1 : void discard_value(const TimeStepId& id_to_discard);
651 :
652 : /// Drop the oldest step (not substep) entry in the history. Any
653 : /// memory allocations will be cached for future reuse.
654 1 : void pop_front();
655 :
656 : /// Drop the newest step or substep entry in the history. Any
657 : /// memory allocations will be cached for future reuse.
658 1 : void undo_latest();
659 :
660 : /// Remove all substep entries from the history. Any memory
661 : /// allocations will be cached for future reuse.
662 1 : void clear_substeps();
663 :
664 : /// Remove all (step and substep) entries from the history. Any
665 : /// memory allocations will be cached for future reuse.
666 1 : void clear();
667 :
668 : /// Release any cached memory allocations.
669 1 : void shrink_to_fit();
670 :
671 : /// Apply \p func to `make_not_null(&e)` for `e` every `derivative`
672 : /// and valid `*value` in records held by the history.
673 : template <typename F>
674 1 : void map_entries(F&& func);
675 :
676 : // NOLINTNEXTLINE(google-runtime-references)
677 0 : void pup(PUP::er& p);
678 :
679 0 : std::ostream& print(std::ostream& os) const;
680 :
681 : private:
682 0 : void discard_value(gsl::not_null<std::optional<Vars>*> value);
683 0 : void cache_allocations(gsl::not_null<StepRecord<Vars>*> record);
684 :
685 : template <typename ValueFunc, typename DerivativeFunc>
686 0 : void insert_impl(const TimeStepId& time_step_id, ValueFunc&& value_inserter,
687 : DerivativeFunc&& derivative_inserter);
688 :
689 : template <typename InsertedVars>
690 0 : void insert_initial_impl(TimeStepId time_step_id, InsertedVars value,
691 : DerivVars derivative);
692 :
693 0 : size_t integration_order_{0};
694 :
695 0 : StaticDeque<StepRecord<Vars>, max_size()> step_values_{};
696 : boost::container::static_vector<StepRecord<Vars>, history_max_substeps>
697 0 : substep_values_{};
698 :
699 : // If the algorithm undoes a substep, it may want to reset the
700 : // variables to the previous value. We don't want the time stepper
701 : // to have to keep track of this when discarding values, so we hang
702 : // onto the last value if it gets discarded.
703 0 : std::optional<Vars> latest_value_if_discarded_{};
704 :
705 : // Memory allocations available for reuse.
706 : boost::container::static_vector<Vars, max_size() + history_max_substeps>
707 0 : vars_allocation_cache_{};
708 : boost::container::static_vector<DerivVars, max_size() + history_max_substeps>
709 0 : deriv_vars_allocation_cache_{};
710 : };
711 :
712 : // Don't copy the allocation caches.
713 : template <typename Vars>
714 : History<Vars>::History(const History& other)
715 : : integration_order_(other.integration_order_),
716 : step_values_(other.step_values_),
717 : substep_values_(other.substep_values_),
718 : latest_value_if_discarded_(other.latest_value_if_discarded_) {}
719 :
720 : // Don't copy the allocation caches.
721 : template <typename Vars>
722 : History<Vars>& History<Vars>::operator=(const History& other) {
723 : integration_order_ = other.integration_order_;
724 : step_values_ = other.step_values_;
725 : substep_values_ = other.substep_values_;
726 : latest_value_if_discarded_ = other.latest_value_if_discarded_;
727 : return *this;
728 : }
729 :
730 : template <typename Vars>
731 : void History<Vars>::insert(const TimeStepId& time_step_id, const Vars& value,
732 : const DerivVars& derivative) {
733 : insert_impl(
734 : time_step_id,
735 : [&](const gsl::not_null<Vars*> record_value) { *record_value = value; },
736 : [&](const gsl::not_null<DerivVars*> record_derivative) {
737 : *record_derivative = derivative;
738 : });
739 : }
740 :
741 : template <typename Vars>
742 : void History<Vars>::insert(const TimeStepId& time_step_id, NoValue /*unused*/,
743 : const DerivVars& derivative) {
744 : insert_impl(time_step_id, NoValue{},
745 : [&](const gsl::not_null<DerivVars*> record_derivative) {
746 : *record_derivative = derivative;
747 : });
748 : }
749 :
750 : template <typename Vars>
751 : template <typename ValueFunc, typename DerivativeFunc>
752 : void History<Vars>::insert_in_place(const TimeStepId& time_step_id,
753 : ValueFunc&& value_inserter,
754 : DerivativeFunc&& derivative_inserter) {
755 : insert_impl(time_step_id, std::forward<ValueFunc>(value_inserter),
756 : std::forward<DerivativeFunc>(derivative_inserter));
757 : }
758 :
759 : template <typename Vars>
760 : template <typename DerivativeFunc>
761 : void History<Vars>::insert_in_place(const TimeStepId& time_step_id,
762 : NoValue /*unused*/,
763 : DerivativeFunc&& derivative_inserter) {
764 : insert_impl(time_step_id, NoValue{},
765 : std::forward<DerivativeFunc>(derivative_inserter));
766 : }
767 :
768 : template <typename Vars>
769 : void History<Vars>::insert_initial(TimeStepId time_step_id, Vars value,
770 : DerivVars derivative) {
771 : insert_initial_impl(std::move(time_step_id), std::move(value),
772 : std::move(derivative));
773 : }
774 :
775 : template <typename Vars>
776 : void History<Vars>::insert_initial(TimeStepId time_step_id, NoValue /*unused*/,
777 : DerivVars derivative) {
778 : insert_initial_impl(std::move(time_step_id), NoValue{},
779 : std::move(derivative));
780 : }
781 :
782 : template <typename Vars>
783 : const Vars& History<Vars>::latest_value() const {
784 : const auto& latest_record =
785 : at_step_start() ? this->back() : substeps().back();
786 : if (latest_record.value.has_value()) {
787 : return *latest_record.value;
788 : } else {
789 : ASSERT(latest_value_if_discarded_.has_value(),
790 : "Latest value unavailable. The latest insertion was undone.");
791 : return *latest_value_if_discarded_;
792 : }
793 : }
794 :
795 : template <typename Vars>
796 : const StepRecord<Vars>& History<Vars>::step_start(const double time) const {
797 : // Search starting at the end to handle self-start correctly (and
798 : // because the result under usual use is one of the last two
799 : // entries).
800 : const auto first_step = this->begin();
801 : auto step = this->end();
802 : ASSERT(step != first_step, "History is empty");
803 : --step;
804 : const evolution_less<double> before{step->time_step_id.time_runs_forward()};
805 : while (before(time, step->time_step_id.step_time().value())) {
806 : ASSERT(step != first_step,
807 : "Start of step at time " << time << " is before start of history "
808 : << first_step->time_step_id);
809 : --step;
810 : }
811 : return *step;
812 : }
813 :
814 : template <typename Vars>
815 : bool History<Vars>::at_step_start() const {
816 : ASSERT(not this->empty(), "History is empty");
817 : return substep_values_.empty() or
818 : substep_values_.back().time_step_id < this->back().time_step_id;
819 : }
820 :
821 : template <typename Vars>
822 : void History<Vars>::discard_value(const TimeStepId& id_to_discard) {
823 : auto& latest_record = at_step_start() ? this->back() : substeps().back();
824 : if (latest_record.time_step_id == id_to_discard) {
825 : discard_value(&latest_value_if_discarded_);
826 : latest_value_if_discarded_ = std::move(latest_record.value);
827 : latest_record.value.reset();
828 : } else {
829 : discard_value(&History_detail::find_record(*this, id_to_discard).value);
830 : }
831 : }
832 :
833 : template <typename Vars>
834 : void History<Vars>::pop_front() {
835 : ASSERT(not this->empty(), "History is empty");
836 : ASSERT(substeps().empty() or substeps().front().time_step_id.step_time() !=
837 : this->front().time_step_id.step_time(),
838 : "Cannot remove a step with substeps. Call clear_substeps() first.");
839 : cache_allocations(&step_values_.front());
840 : step_values_.pop_front();
841 : }
842 :
843 : template <typename Vars>
844 : void History<Vars>::undo_latest() {
845 : if (at_step_start()) {
846 : cache_allocations(&step_values_.back());
847 : step_values_.pop_back();
848 : } else {
849 : cache_allocations(&substep_values_.back());
850 : substep_values_.pop_back();
851 : }
852 : discard_value(&latest_value_if_discarded_);
853 : }
854 :
855 : template <typename Vars>
856 : void History<Vars>::clear_substeps() {
857 : for (auto& record : substep_values_) {
858 : cache_allocations(&record);
859 : }
860 : substep_values_.clear();
861 : }
862 :
863 : template <typename Vars>
864 : void History<Vars>::clear() {
865 : clear_substeps();
866 : while (not this->empty()) {
867 : pop_front();
868 : }
869 : discard_value(&latest_value_if_discarded_);
870 : }
871 :
872 : template <typename Vars>
873 : void History<Vars>::shrink_to_fit() {
874 : vars_allocation_cache_.clear();
875 : vars_allocation_cache_.shrink_to_fit();
876 : deriv_vars_allocation_cache_.clear();
877 : deriv_vars_allocation_cache_.shrink_to_fit();
878 : }
879 :
880 : template <typename Vars>
881 : template <typename F>
882 : void History<Vars>::map_entries(F&& func) {
883 : for (auto& record : *this) {
884 : func(make_not_null(&record.derivative));
885 : if (record.value.has_value()) {
886 : func(make_not_null(&*record.value));
887 : }
888 : }
889 : for (auto& record : this->substeps()) {
890 : func(make_not_null(&record.derivative));
891 : if (record.value.has_value()) {
892 : func(make_not_null(&*record.value));
893 : }
894 : }
895 : }
896 :
897 : template <typename Vars>
898 : void History<Vars>::pup(PUP::er& p) {
899 : p | integration_order_;
900 : p | step_values_;
901 : p | substep_values_;
902 : p | latest_value_if_discarded_;
903 :
904 : // Don't serialize the allocation cache.
905 : }
906 :
907 : template <typename Vars>
908 : std::ostream& History<Vars>::print(std::ostream& os) const {
909 : using ::operator<<;
910 : os << "Integration order: " << integration_order_ << "\n";
911 : os << "Step values:\n";
912 : for (auto& record : *this) {
913 : record.print(os);
914 : }
915 : os << "Substep values:\n";
916 : for (auto& record : substep_values_) {
917 : record.print(os);
918 : }
919 : os << "Latest value if discarded: " << latest_value_if_discarded_ << "\n";
920 : return os;
921 : }
922 :
923 : // Doxygen is confused by this function for some reason.
924 : /// \cond
925 : template <typename Vars>
926 : void History<Vars>::discard_value(
927 : const gsl::not_null<std::optional<Vars>*> value) {
928 : if (not value->has_value()) {
929 : return;
930 : }
931 : // If caching doesn't save anything, don't allocate memory for the cache.
932 : if (contains_allocations(**value)) {
933 : vars_allocation_cache_.emplace_back(std::move(**value));
934 : }
935 : value->reset();
936 : }
937 : /// \endcond
938 :
939 : template <typename Vars>
940 : void History<Vars>::cache_allocations(
941 : const gsl::not_null<StepRecord<Vars>*> record) {
942 : discard_value(&record->value);
943 : // If caching doesn't save anything, don't allocate memory for the cache.
944 : if (contains_allocations(record->derivative)) {
945 : deriv_vars_allocation_cache_.emplace_back(std::move(record->derivative));
946 : }
947 : }
948 :
949 : template <typename Vars>
950 : template <typename ValueFunc, typename DerivativeFunc>
951 : void History<Vars>::insert_impl(const TimeStepId& time_step_id,
952 : ValueFunc&& value_inserter,
953 : DerivativeFunc&& derivative_inserter) {
954 : ASSERT(this->empty() or time_step_id > this->back().time_step_id,
955 : "New entry at " << time_step_id
956 : << " must be later than previous entry at "
957 : << this->back().time_step_id);
958 : discard_value(&latest_value_if_discarded_);
959 : StepRecord<Vars> record{};
960 : record.time_step_id = time_step_id;
961 : if constexpr (not std::is_same_v<ValueFunc, NoValue>) {
962 : if (vars_allocation_cache_.empty()) {
963 : record.value.emplace();
964 : } else {
965 : record.value.emplace(std::move(vars_allocation_cache_.back()));
966 : vars_allocation_cache_.pop_back();
967 : }
968 : std::forward<ValueFunc>(value_inserter)(make_not_null(&*record.value));
969 : }
970 : if (not deriv_vars_allocation_cache_.empty()) {
971 : record.derivative = std::move(deriv_vars_allocation_cache_.back());
972 : deriv_vars_allocation_cache_.pop_back();
973 : }
974 : std::forward<DerivativeFunc>(derivative_inserter)(
975 : make_not_null(&record.derivative));
976 :
977 : const size_t substep = time_step_id.substep();
978 : if (substep == 0) {
979 : step_values_.push_back(std::move(record));
980 : } else {
981 : ASSERT(not this->empty(), "Cannot insert substep into empty history.");
982 : ASSERT(time_step_id.step_time() == this->back().time_step_id.step_time(),
983 : "Cannot insert substep " << time_step_id << " of different step "
984 : << this->back().time_step_id);
985 : ASSERT(substep == substeps().size() + 1,
986 : "Cannot insert substep " << substep << " following "
987 : << substeps().size());
988 : ASSERT(substep_values_.size() < substep_values_.max_size(),
989 : "Cannot insert new substep because the History is full.");
990 : substep_values_.push_back(std::move(record));
991 : }
992 : }
993 :
994 : template <typename Vars>
995 : template <typename InsertedVars>
996 : void History<Vars>::insert_initial_impl(TimeStepId time_step_id,
997 : InsertedVars value,
998 : DerivVars derivative) {
999 : ASSERT(vars_allocation_cache_.empty(),
1000 : "insert_initial should only be used for initialization");
1001 : ASSERT(deriv_vars_allocation_cache_.empty(),
1002 : "insert_initial should only be used for initialization");
1003 : ASSERT(time_step_id.substep() == 0, "Cannot use insert_initial for substeps");
1004 : ASSERT(this->empty() or time_step_id < this->front().time_step_id,
1005 : "New initial entry at " << time_step_id
1006 : << " must be earlier than previous entry at "
1007 : << this->front().time_step_id);
1008 :
1009 : if constexpr (std::is_same_v<InsertedVars, Vars>) {
1010 : step_values_.push_front(StepRecord<Vars>{
1011 : std::move(time_step_id), std::move(value), std::move(derivative)});
1012 : } else {
1013 : static_assert(std::is_same_v<InsertedVars, NoValue>);
1014 : (void)value;
1015 : step_values_.push_front(StepRecord<Vars>{
1016 : std::move(time_step_id), std::nullopt, std::move(derivative)});
1017 : }
1018 : }
1019 :
1020 : template <typename Vars>
1021 0 : bool operator==(const History<Vars>& a, const History<Vars>& b) {
1022 : return a.integration_order() == b.integration_order() and
1023 : a.size() == b.size() and std::equal(a.begin(), a.end(), b.begin()) and
1024 : a.substeps() == b.substeps();
1025 : }
1026 :
1027 : template <typename Vars>
1028 0 : bool operator!=(const History<Vars>& a, const History<Vars>& b) {
1029 : return not(a == b);
1030 : }
1031 :
1032 : template <typename Vars>
1033 0 : std::ostream& operator<<(std::ostream& os, const History<Vars>& history) {
1034 : return history.print(os);
1035 : }
1036 :
1037 : /// \ingroup TimeSteppersGroup
1038 : /// Initialize a History object based on the contents of another,
1039 : /// applying a transformation to each value and derivative.
1040 : ///
1041 : /// The transformation functions can either take a value from the
1042 : /// source history and return a value for the destination history or
1043 : /// take a `gsl::not_null` value from the destination history and a
1044 : /// value from the source history to initialize it with. For the sake
1045 : /// of implementation simplicity, either both transformers must mutate
1046 : /// or both must produce values.
1047 : ///
1048 : /// An overload applying the same transformation to the values and
1049 : /// derivatives is provided for convenience.
1050 : ///
1051 : /// \see transform_mutate
1052 : /// @{
1053 : template <typename DestVars, typename SourceVars, typename ValueTransformer,
1054 : typename DerivativeTransformer>
1055 1 : void transform(const gsl::not_null<History<DestVars>*> dest,
1056 : const History<SourceVars>& source,
1057 : ValueTransformer&& value_transformer,
1058 : DerivativeTransformer&& derivative_transformer) {
1059 : dest->clear_substeps();
1060 : dest->clear();
1061 : dest->integration_order(source.integration_order());
1062 : if (source.empty()) {
1063 : return;
1064 : }
1065 : auto pre_substep_end = source.end();
1066 : if (not source.substeps().empty() and
1067 : source.back().time_step_id > source.substeps().back().time_step_id) {
1068 : --pre_substep_end;
1069 : }
1070 :
1071 : const auto transform_record =
1072 : [&derivative_transformer, &dest, &value_transformer](
1073 : const typename History<SourceVars>::value_type& record) {
1074 : if constexpr (std::is_invocable_v<ValueTransformer,
1075 : gsl::not_null<DestVars*>,
1076 : const SourceVars&>) {
1077 : if (record.value.has_value()) {
1078 : dest->insert_in_place(
1079 : record.time_step_id,
1080 : [&](const auto result) {
1081 : value_transformer(result, *record.value);
1082 : },
1083 : [&](const auto result) {
1084 : derivative_transformer(result, record.derivative);
1085 : });
1086 : } else {
1087 : dest->insert_in_place(
1088 : record.time_step_id, History<DestVars>::no_value,
1089 : [&](const auto result) {
1090 : derivative_transformer(result, record.derivative);
1091 : });
1092 : }
1093 : } else {
1094 : static_assert(
1095 : std::is_invocable_v<ValueTransformer, const SourceVars&>,
1096 : "Transform function must either be callable to mutate entries "
1097 : "or return the transformed state by value.");
1098 : if (record.value.has_value()) {
1099 : dest->insert(record.time_step_id, value_transformer(*record.value),
1100 : derivative_transformer(record.derivative));
1101 : } else {
1102 : dest->insert(record.time_step_id, History<DestVars>::no_value,
1103 : derivative_transformer(record.derivative));
1104 : }
1105 : }
1106 : };
1107 :
1108 : auto copying_step = source.begin();
1109 : for (; copying_step != pre_substep_end; ++copying_step) {
1110 : transform_record(*copying_step);
1111 : }
1112 : for (const auto& record : source.substeps()) {
1113 : transform_record(record);
1114 : }
1115 : if (pre_substep_end != source.end()) {
1116 : transform_record(*pre_substep_end);
1117 : }
1118 : }
1119 :
1120 : template <typename DestVars, typename SourceVars, typename Transformer>
1121 1 : void transform(const gsl::not_null<History<DestVars>*> dest,
1122 : const History<SourceVars>& source, Transformer&& transformer) {
1123 : transform(dest, source, transformer, transformer);
1124 : }
1125 : /// @}
1126 :
1127 : /// \ingroup TimeSteppersGroup
1128 : /// Combine two History objects by applying a transformation to each
1129 : /// value and derivative.
1130 : ///
1131 : /// The transformers must accept a `gsl::not_null` value from the \p
1132 : /// dest history and a const reference to value from the \p source
1133 : /// history. They will be applied to each corresponding (non-nullopt)
1134 : /// pair of values from the two histories. It is an error if the
1135 : /// entries in the two histories do not match.
1136 : ///
1137 : /// An overload applying the same transformation to the values and
1138 : /// derivatives is provided for convenience.
1139 : ///
1140 : /// \see transform
1141 : /// @{
1142 : template <typename DestVars, typename SourceVars, typename ValueTransformer,
1143 : typename DerivativeTransformer>
1144 1 : void transform_mutate(const gsl::not_null<History<DestVars>*> dest,
1145 : const History<SourceVars>& source,
1146 : ValueTransformer&& value_transformer,
1147 : DerivativeTransformer&& derivative_transformer) {
1148 : ASSERT(dest->integration_order() == source.integration_order(),
1149 : "Attempting to combine histories with integration orders "
1150 : << dest->integration_order() << " and "
1151 : << source.integration_order());
1152 :
1153 : const auto transform_record =
1154 : [&derivative_transformer, &value_transformer](
1155 : typename History<DestVars>::value_type& dest_record,
1156 : const typename History<SourceVars>::value_type& source_record) {
1157 : ASSERT(dest_record.time_step_id == source_record.time_step_id,
1158 : "Entries to combine do not match: "
1159 : << dest_record.time_step_id << " "
1160 : << source_record.time_step_id);
1161 : ASSERT(dest_record.value.has_value() == source_record.value.has_value(),
1162 : "Only one of the entries to combine has a value at "
1163 : << dest_record.time_step_id);
1164 : if (dest_record.value.has_value()) {
1165 : value_transformer(make_not_null(&*dest_record.value),
1166 : *source_record.value);
1167 : }
1168 : derivative_transformer(make_not_null(&dest_record.derivative),
1169 : source_record.derivative);
1170 : };
1171 :
1172 : ASSERT(dest->size() == source.size(),
1173 : "Attempting to combine histories with sizes "
1174 : << dest->size() << " and " << source.size());
1175 : {
1176 : auto dest_entry = dest->begin();
1177 : auto source_entry = source.begin();
1178 : while (dest_entry != dest->end()) {
1179 : transform_record(*dest_entry, *source_entry);
1180 : ++dest_entry;
1181 : ++source_entry;
1182 : }
1183 : }
1184 :
1185 : {
1186 : auto dest_substeps = dest->substeps();
1187 : const auto source_substeps = source.substeps();
1188 : ASSERT(dest_substeps.size() == source_substeps.size(),
1189 : "Attempting to combine histories with substep sizes "
1190 : << dest_substeps.size() << " and " << source_substeps.size());
1191 : auto dest_entry = dest_substeps.begin();
1192 : auto source_entry = source_substeps.begin();
1193 : while (dest_entry != dest_substeps.end()) {
1194 : transform_record(*dest_entry, *source_entry);
1195 : ++dest_entry;
1196 : ++source_entry;
1197 : }
1198 : }
1199 : }
1200 :
1201 : template <typename DestVars, typename SourceVars, typename Transformer>
1202 1 : void transform_mutate(const gsl::not_null<History<DestVars>*> dest,
1203 : const History<SourceVars>& source,
1204 : Transformer&& transformer) {
1205 : transform_mutate(dest, source, transformer, transformer);
1206 : }
1207 : /// @}
1208 : } // namespace TimeSteppers
|
