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