AdamsBashforthN.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 /// \file
5 /// Defines class AdamsBashforthN
6 
7 #pragma once
8 
9 #include <algorithm>
10 #include <boost/iterator/transform_iterator.hpp>
11 #include <cstddef>
12 #include <iosfwd>
13 #include <iterator>
14 #include <limits>
15 #include <map>
16 #include <pup.h>
17 #include <tuple>
18 #include <type_traits>
19 #include <vector>
20 
21 #include "ErrorHandling/Assert.hpp"
22 #include "ErrorHandling/Error.hpp"
24 #include "Options/Options.hpp"
26 #include "Time/EvolutionOrdering.hpp"
27 #include "Time/Time.hpp"
28 #include "Time/TimeId.hpp"
29 #include "Time/TimeSteppers/TimeStepper.hpp" // IWYU pragma: keep
32 #include "Utilities/Gsl.hpp"
34 #include "Utilities/Overloader.hpp"
35 #include "Utilities/TMPL.hpp"
36 
37 /// \cond
38 namespace TimeSteppers {
39 template <typename LocalVars, typename RemoteVars, typename CouplingResult>
40 class BoundaryHistory; // IWYU pragma: keep
41 template <typename Vars, typename DerivVars>
42 class History;
43 } // namespace TimeSteppers
44 /// \endcond
45 
46 namespace TimeSteppers {
47 
48 /// \ingroup TimeSteppersGroup
49 ///
50 /// An Nth Adams-Bashforth time stepper.
51 class AdamsBashforthN : public LtsTimeStepper::Inherit {
52  public:
53  static constexpr const size_t maximum_order = 8;
54 
55  struct Order {
56  using type = size_t;
57  static constexpr OptionString help = {"Convergence order"};
58  static type lower_bound() noexcept { return 1; }
59  static type upper_bound() noexcept { return maximum_order; }
60  };
61  using options = tmpl::list<Order>;
62  static constexpr OptionString help = {
63  "An Adams-Bashforth Nth order time-stepper."};
64 
65  AdamsBashforthN() = default;
66  explicit AdamsBashforthN(size_t order) noexcept;
67  AdamsBashforthN(const AdamsBashforthN&) noexcept = default;
68  AdamsBashforthN& operator=(const AdamsBashforthN&) noexcept = default;
69  AdamsBashforthN(AdamsBashforthN&&) noexcept = default;
70  AdamsBashforthN& operator=(AdamsBashforthN&&) noexcept = default;
71  ~AdamsBashforthN() noexcept override = default;
72 
73  template <typename Vars, typename DerivVars>
74  void update_u(gsl::not_null<Vars*> u,
76  const TimeDelta& time_step) const noexcept;
77 
78  template <typename Vars, typename DerivVars>
79  void dense_update_u(gsl::not_null<Vars*> u,
80  const History<Vars, DerivVars>& history,
81  double time) const noexcept;
82 
83  // This is defined as a separate type alias to keep the doxygen page
84  // width somewhat under control.
85  template <typename LocalVars, typename RemoteVars, typename Coupling>
86  using BoundaryHistoryType =
87  BoundaryHistory<LocalVars, RemoteVars,
89 
90  /*!
91  * An explanation of the computation being performed by this
92  * function:
93  * \f$\newcommand\tL{t^L}\newcommand\tR{t^R}\newcommand\tU{\tilde{t}\!}
94  * \newcommand\mat{\mathbf}\f$
95  *
96  * Suppose the local and remote sides of the interface are evaluated
97  * at times \f$\ldots, \tL_{-1}, \tL_0, \tL_1, \ldots\f$ and
98  * \f$\ldots, \tR_{-1}, \tR_0, \tR_1, \ldots\f$, respectively, with
99  * the starting location of the numbering arbitrary in each case.
100  * Let the step we wish to calculate the effect of be the step from
101  * \f$\tL_{m_S}\f$ to \f$\tL_{m_S+1}\f$. We call the sequence
102  * produced from the union of the local and remote time sequences
103  * \f$\ldots, \tU_{-1}, \tU_0, \tU_1, \ldots\f$. For example, one
104  * possible sequence of times is:
105  * \f{equation}
106  * \begin{aligned}
107  * \text{Local side:} \\ \text{Union times:} \\ \text{Remote side:}
108  * \end{aligned}
109  * \cdots
110  * \begin{gathered}
111  * \, \\ \tU_1 \\ \tR_5
112  * \end{gathered}
113  * \leftarrow \Delta \tU_1 \rightarrow
114  * \begin{gathered}
115  * \tL_4 \\ \tU_2 \\ \,
116  * \end{gathered}
117  * \leftarrow \Delta \tU_2 \rightarrow
118  * \begin{gathered}
119  * \, \\ \tU_3 \\ \tR_6
120  * \end{gathered}
121  * \leftarrow \Delta \tU_3 \rightarrow
122  * \begin{gathered}
123  * \, \\ \tU_4 \\ \tR_7
124  * \end{gathered}
125  * \leftarrow \Delta \tU_4 \rightarrow
126  * \begin{gathered}
127  * \tL_5 \\ \tU_5 \\ \,
128  * \end{gathered}
129  * \cdots
130  * \f}
131  * We call the indices of the step's start and end times in the
132  * union time sequence \f$n_S\f$ and \f$n_E\f$, respectively. We
133  * define \f$n^L_m\f$ to be the union-time index corresponding to
134  * \f$\tL_m\f$ and \f$m^L_n\f$ to be the index of the last local
135  * time not later than \f$\tU_n\f$ and similarly for the remote
136  * side. So for the above example, \f$n^L_4 = 2\f$ and \f$m^R_2 =
137  * 5\f$, and if we wish to compute the step from \f$\tL_4\f$ to
138  * \f$\tL_5\f$ we would have \f$m_S = 4\f$, \f$n_S = 2\f$, and
139  * \f$n_E = 5\f$.
140  *
141  * If we wish to evaluate the change over this step to \f$k\f$th
142  * order, we can write the change in the value as a linear
143  * combination of the values of the coupling between the elements at
144  * unequal times:
145  * \f{equation}
146  * \mat{F}_{m_S} =
147  * \mspace{-10mu}
148  * \sum_{q^L = m_S-(k-1)}^{m_S}
149  * \,
150  * \sum_{q^R = m^R_{n_S}-(k-1)}^{m^R_{n_E-1}}
151  * \mspace{-10mu}
152  * \mat{D}_{q^Lq^R}
153  * I_{q^Lq^R},
154  * \f}
155  * where \f$\mat{D}_{q^Lq^R}\f$ is the coupling function evaluated
156  * between data from \f$\tL_{q^L}\f$ and \f$\tR_{q^R}\f$. The
157  * coefficients can be written as the sum of three terms,
158  * \f{equation}
159  * I_{q^Lq^R} = I^E_{q^Lq^R} + I^R_{q^Lq^R} + I^L_{q^Lq^R},
160  * \f}
161  * which can be interpreted as a contribution from equal-time
162  * evaluations and contributions related to the remote and local
163  * evaluation times. These are given by
164  * \f{align}
165  * I^E_{q^Lq^R} &=
166  * \mspace{-10mu}
167  * \sum_{n=n_S}^{\min\left\{n_E, n^L+k\right\}-1}
168  * \mspace{-10mu}
169  * \tilde{\alpha}_{n,n-n^L} \Delta \tU_n
170  * &&\text{if $\tL_{q^L} = \tR_{q^R}$, otherwise 0}
171  * \\
172  * I^R_{q^Lq^R} &=
173  * \ell_{q^L - m_S + k}\!\left(
174  * \tU_{n^R}; \tL_{m_S - (k-1)}, \ldots, \tL_{m_S}\right)
175  * \mspace{-10mu}
176  * \sum_{n=\max\left\{n_S, n^R\right\}}
177  * ^{\min\left\{n_E, n^R+k\right\}-1}
178  * \mspace{-10mu}
179  * \tilde{\alpha}_{n,n-n^R} \Delta \tU_n
180  * &&\text{if $\tR_{q^R}$ is not in $\{\tL_{\vphantom{|}\cdots}\}$,
181  * otherwise 0}
182  * \\
183  * I^L_{q^Lq^R} &=
184  * \mspace{-10mu}
185  * \sum_{n=\max\left\{n_S, n^R\right\}}
186  * ^{\min\left\{n_E, n^L+k, n^R_{q^R+k}\right\}-1}
187  * \mspace{-10mu}
188  * \ell_{q^R - m^R_n + k}\!\left(\tU_{n^L};
189  * \tR_{m^R_n - (k-1)}, \ldots, \tR_{m^R_n}\right)
190  * \tilde{\alpha}_{n,n-n^L} \Delta \tU_n
191  * &&\text{if $\tL_{q^L}$ is not in $\{\tR_{\vphantom{|}\cdots}\}$,
192  * otherwise 0,}
193  * \f}
194  * where for brevity we write \f$n^L = n^L_{q^L}\f$ and \f$n^R =
195  * n^R_{q^R}\f$, and where \f$\ell_a(t; x_1, \ldots, x_k)\f$ a
196  * Lagrange interpolating polynomial and \f$\tilde{\alpha}_{nj}\f$
197  * is the \f$j\f$th coefficient for an Adams-Bashforth step over the
198  * union times from step \f$n\f$ to step \f$n+1\f$.
199  */
200  template <typename LocalVars, typename RemoteVars, typename Coupling>
201  std::result_of_t<const Coupling&(LocalVars, RemoteVars)>
203  const Coupling& coupling,
205  history,
206  const TimeDelta& time_step) const noexcept;
207 
208  template <typename LocalVars, typename RemoteVars, typename Coupling>
209  std::result_of_t<const Coupling&(LocalVars, RemoteVars)>
210  boundary_dense_output(
211  const Coupling& coupling,
213  double time) const noexcept;
214 
215  size_t number_of_past_steps() const noexcept override;
216 
217  double stable_step() const noexcept override;
218 
219  TimeId next_time_id(const TimeId& current_id,
220  const TimeDelta& time_step) const noexcept override;
221 
222  template <typename Vars, typename DerivVars>
223  bool can_change_step_size(
224  const TimeId& time_id,
225  const TimeSteppers::History<Vars, DerivVars>& history) const noexcept;
226 
228 
229  explicit AdamsBashforthN(CkMigrateMessage* /*unused*/) noexcept {}
230 
231  // clang-tidy: do not pass by non-const reference
232  void pup(PUP::er& p) noexcept override; // NOLINT
233 
234  private:
235  friend bool operator==(const AdamsBashforthN& lhs,
236  const AdamsBashforthN& rhs) noexcept;
237 
238  // Some of the private methods take a parameter of type "Delta" or
239  // "TimeType". Delta is expected to be a TimeDelta or an
240  // ApproximateTimeDelta, and TimeType is expected to be a Time or an
241  // ApproximateTime. The former cases will detect and optimize the
242  // constant-time-step case, while the latter are necessary for dense
243  // output.
244 
245  template <typename Vars, typename DerivVars, typename Delta>
246  void update_u_impl(gsl::not_null<Vars*> u,
247  const History<Vars, DerivVars>& history,
248  const Delta& time_step) const noexcept;
249 
250  template <typename LocalVars, typename RemoteVars, typename Coupling,
251  typename TimeType>
252  std::result_of_t<const Coupling&(LocalVars, RemoteVars)> boundary_impl(
253  const Coupling& coupling,
255  const TimeType& end_time) const noexcept;
256 
257  /// Get coefficients for a time step. Arguments are an iterator
258  /// pair to past times, oldest to newest, and the time step to take.
259  template <typename Iterator, typename Delta>
260  static std::vector<double> get_coefficients(const Iterator& times_begin,
261  const Iterator& times_end,
262  const Delta& step) noexcept;
263 
264  static std::vector<double> get_coefficients_impl(
265  const std::vector<double>& steps) noexcept;
266 
267  static std::vector<double> variable_coefficients(
268  const std::vector<double>& steps) noexcept;
269 
270  static std::vector<double> constant_coefficients(size_t order) noexcept;
271 
272  struct ApproximateTimeDelta;
273 
274  // Time-like interface to a double used for dense output
275  struct ApproximateTime {
277  double value() const noexcept { return time; }
278 
279  // Only the operators that are actually used are defined.
280  friend ApproximateTimeDelta operator-(const ApproximateTime& a,
281  const Time& b) noexcept {
282  return {a.value() - b.value()};
283  }
284 
285  friend bool operator<(const Time& a, const ApproximateTime& b) noexcept {
286  return a.value() < b.value();
287  }
288 
289  friend bool operator<(const ApproximateTime& a, const Time& b) noexcept {
290  return a.value() < b.value();
291  }
292 
293  friend std::ostream& operator<<(std::ostream& s,
294  const ApproximateTime& t) noexcept {
295  return s << t.value();
296  }
297  };
298 
299  // TimeDelta-like interface to a double used for dense output
300  struct ApproximateTimeDelta {
302  double value() const noexcept { return delta; }
303  bool is_positive() const noexcept { return delta > 0.; }
304 
305  // Only the operators that are actually used are defined.
306  friend bool operator<(const ApproximateTimeDelta& a,
307  const ApproximateTimeDelta& b) noexcept {
308  return a.value() < b.value();
309  }
310 
311  friend double operator/(
312  const TimeDelta& a,
313  const AdamsBashforthN::ApproximateTimeDelta& b) noexcept {
314  return a.value() / b.value();
315  }
316  };
317 
318  size_t order_ = 3;
319 };
320 
321 bool operator!=(const AdamsBashforthN& lhs,
322  const AdamsBashforthN& rhs) noexcept;
323 
324 template <typename Vars, typename DerivVars>
325 void AdamsBashforthN::update_u(
326  const gsl::not_null<Vars*> u,
327  const gsl::not_null<History<Vars, DerivVars>*> history,
328  const TimeDelta& time_step) const noexcept {
329  update_u_impl(u, *history, time_step);
330  history->mark_unneeded(history->begin() + 1);
331 }
332 
333 template <typename Vars, typename DerivVars>
334 void AdamsBashforthN::dense_update_u(const gsl::not_null<Vars*> u,
335  const History<Vars, DerivVars>& history,
336  const double time) const noexcept {
337  const ApproximateTimeDelta time_step{
338  time - history[history.size() - 1].value()};
339  update_u_impl(u, history, time_step);
340 }
341 
342 template <typename Vars, typename DerivVars, typename Delta>
343 void AdamsBashforthN::update_u_impl(const gsl::not_null<Vars*> u,
344  const History<Vars, DerivVars>& history,
345  const Delta& time_step) const noexcept {
346  ASSERT(history.size() <= order_,
347  "Length of history (" << history.size() << ") "
348  << "should not exceed target order (" << order_ << ")");
349 
350  const auto& coefficients =
351  get_coefficients(history.begin(), history.end(), time_step);
352 
353  const auto do_update =
354  [u, &time_step, &coefficients, &history](auto order) noexcept {
355  *u += time_step.value() * constexpr_sum<order>(
356  [order, &coefficients, &history](auto i) noexcept {
357  return coefficients[order - 1 - i] *
358  (history.begin() +
359  static_cast<
360  typename History<Vars, DerivVars>::difference_type>(i))
361  .derivative();
362  });
363  };
364 
365  switch (history.size()) {
366  case 1:
368  break;
369  case 2:
371  break;
372  case 3:
374  break;
375  case 4:
377  break;
378  case 5:
380  break;
381  case 6:
383  break;
384  case 7:
386  break;
387  case 8:
389  break;
390  default:
391  ERROR("Bad amount of history data: " << history.size());
392  }
393 }
394 
395 template <typename LocalVars, typename RemoteVars, typename Coupling>
398  const Coupling& coupling,
400  history,
401  const TimeDelta& time_step) const noexcept {
402  auto result = boundary_impl(coupling, *history,
403  *(history->local_end() - 1) + time_step);
404 
405  // We know that the local side will step at end_time, so the step
406  // containing that time will be the next step, which is not
407  // currently in the history. We therefore know we won't need the
408  // oldest value for the next step.
409  history->local_mark_unneeded(history->local_begin() + 1);
410  // We don't know whether the remote side will step at end_time, so
411  // we have to be conservative and assume it will not. If it does we
412  // will ignore the first value in the next call to this function.
413  history->remote_mark_unneeded(
414  history->remote_end() -
415  static_cast<typename decltype(history->remote_begin())::difference_type>(
416  history->local_size() + 1));
417 
418  return result;
419 }
420 
421 template <typename LocalVars, typename RemoteVars, typename Coupling>
423 AdamsBashforthN::boundary_dense_output(
424  const Coupling& coupling,
426  const double time) const noexcept {
427  return boundary_impl(coupling, history, ApproximateTime{time});
428 }
429 
430 template <typename LocalVars, typename RemoteVars, typename Coupling,
431  typename TimeType>
433 AdamsBashforthN::boundary_impl(
434  const Coupling& coupling,
436  const TimeType& end_time) const noexcept {
437  // Might be different from order_ during self-start.
438  const auto current_order = history.local_size();
439 
440  ASSERT(current_order <= order_,
441  "Local history is too long for target order (" << current_order
442  << " should not exceed " << order_ << ")");
443  ASSERT(history.remote_size() >= current_order,
444  "Remote history is too short (" << history.remote_size()
445  << " should be at least " << current_order << ")");
446 
447  // Start and end of the step we are trying to take
448  const Time start_time = *(history.local_end() - 1);
449  const auto time_step = end_time - start_time;
450 
451  // If a remote evaluation is done at the start of the step then that
452  // is part of the history for the first union step. When we did
453  // history cleanup at the end of the previous step we didn't know we
454  // were going to get this point so we kept an extra remote history
455  // value.
456  const bool remote_aligned_at_step_start =
457  history.remote_size() > current_order and
458  *(history.remote_begin() +
459  static_cast<typename decltype(history.remote_begin())::difference_type>(
460  current_order)) == start_time;
461  const auto remote_begin =
462  history.remote_begin() + (remote_aligned_at_step_start ? 1 : 0);
463 
464  // Result variable. We evaluate the coupling only for the
465  // structure. This evaluation may be expensive, but by choosing the
466  // most recent times on both sides we should guarantee that it is a
467  // result we need later, so this will serve to get it into the
468  // coupling cache so we don't have to compute it when we actually use it.
469  auto accumulated_change =
470  make_with_value<std::result_of_t<const Coupling&(LocalVars, RemoteVars)>>(
471  history.coupling(coupling, history.local_end() - 1,
472  history.remote_end() - 1),
473  0.);
474 
475  if (history.local_size() ==
476  static_cast<size_t>(history.remote_end() - remote_begin) and
477  std::equal(history.local_begin(), history.local_end(), remote_begin)) {
478  // No local time-stepping going on.
479  const auto coefficients =
480  get_coefficients(history.local_begin(), history.local_end(), time_step);
481 
482  auto local_it = history.local_begin();
483  auto remote_it = remote_begin;
484  for (auto coefficients_it = coefficients.rbegin();
485  coefficients_it != coefficients.rend();
486  ++coefficients_it, ++local_it, ++remote_it) {
487  accumulated_change +=
488  *coefficients_it * history.coupling(coupling, local_it, remote_it);
489  }
490  accumulated_change *= time_step.value();
491 
492  return accumulated_change;
493  }
494 
495  ASSERT(current_order == order_,
496  "Cannot perform local time-stepping while self-starting.");
497 
498  // Avoid billions of casts
499  const auto order_s = static_cast<typename BoundaryHistoryType<
500  LocalVars, RemoteVars, Coupling>::remote_iterator::difference_type>(
501  order_);
502 
503  const evolution_less<> less{time_step.is_positive()};
504 
505  ASSERT(std::is_sorted(history.local_begin(), history.local_end(), less),
506  "Local history not in order");
507  ASSERT(std::is_sorted(remote_begin, history.remote_end(), less),
508  "Remote history not in order");
509  ASSERT(not less(start_time, *(remote_begin + (order_s - 1))),
510  "Remote history does not extend far enough back");
511  ASSERT(less(*(history.remote_end() - 1), end_time),
512  "Please supply only older data: " << *(history.remote_end() - 1)
513  << " is not before " << end_time);
514 
515  // Union of times of all step boundaries on any side.
516  const auto union_times = [&history, &remote_begin, &less]() noexcept {
517  std::vector<Time> ret;
518  ret.reserve(history.local_size() + history.remote_size());
519  std::set_union(history.local_begin(), history.local_end(), remote_begin,
520  history.remote_end(), std::back_inserter(ret), less);
521  return ret;
522  }();
523 
524  using UnionIter = typename decltype(union_times)::const_iterator;
525 
526  // Find the union times iterator for a given time.
527  const auto union_step = [&union_times, &less](const Time& t) noexcept {
528  return std::lower_bound(union_times.cbegin(), union_times.cend(), t, less);
529  };
530 
531  // The union time index for the step start.
532  const auto union_step_start = union_step(start_time);
533 
534  // min(union_times.end(), it + order_s) except being careful not
535  // to create out-of-range iterators.
536  const auto advance_within_step =
537  [order_s, &union_times](const UnionIter& it) noexcept {
538  return union_times.end() - it >
539  static_cast<typename decltype(union_times)::difference_type>(
540  order_s)
541  ? it + static_cast<typename decltype(
542  union_times)::difference_type>(order_s)
543  : union_times.end();
544  };
545 
546  // Calculating the Adams-Bashforth coefficients is somewhat
547  // expensive, so we cache them. ab_coefs(it, step) returns the
548  // coefficients used to step from *it to *it + step.
549  auto ab_coefs = make_overloader(
551  std::map>([order_s](
552  const std::tuple<UnionIter, TimeDelta>& args) noexcept {
553  return get_coefficients(
554  std::get<0>(args) -
555  static_cast<typename UnionIter::difference_type>(order_s - 1),
556  std::get<0>(args) + 1, std::get<1>(args));
557  }),
559  std::map>([order_s](
560  const std::tuple<UnionIter, ApproximateTimeDelta>& args) noexcept {
561  return get_coefficients(
562  std::get<0>(args) -
563  static_cast<typename UnionIter::difference_type>(order_s - 1),
564  std::get<0>(args) + 1, std::get<1>(args));
565  }));
566 
567  // The value of the coefficient of `evaluation_step` when doing
568  // a standard Adams-Bashforth integration over the union times
569  // from `step` to `step + 1`.
570  const auto base_summand = [&ab_coefs, &end_time, &union_times](
571  const UnionIter& step, const UnionIter& evaluation_step) noexcept {
572  if (step + 1 != union_times.end()) {
573  const TimeDelta step_size = *(step + 1) - *step;
574  return step_size.value() *
575  ab_coefs(std::make_tuple(
576  step, step_size))[static_cast<size_t>(step - evaluation_step)];
577  } else {
578  const auto step_size = end_time - *step;
579  return step_size.value() *
580  ab_coefs(std::make_tuple(
581  step, step_size))[static_cast<size_t>(step - evaluation_step)];
582  }
583  };
584 
585  for (auto local_evaluation_step = history.local_begin();
586  local_evaluation_step != history.local_end();
587  ++local_evaluation_step) {
588  const auto union_local_evaluation_step = union_step(*local_evaluation_step);
589  for (auto remote_evaluation_step = remote_begin;
590  remote_evaluation_step != history.remote_end();
591  ++remote_evaluation_step) {
592  double deriv_coef = 0.;
593 
594  if (*local_evaluation_step == *remote_evaluation_step) {
595  // The two elements stepped at the same time. This gives a
596  // standard Adams-Bashforth contribution to each segment
597  // making up the current step.
598  const auto union_step_upper_bound =
599  advance_within_step(union_local_evaluation_step);
600  for (auto step = union_step_start;
601  step < union_step_upper_bound;
602  ++step) {
603  deriv_coef += base_summand(step, union_local_evaluation_step);
604  }
605  } else {
606  // In this block we consider a coupling evaluation that is not
607  // performed at equal times on the two sides of the mortar.
608 
609  // Makes an iterator with a map to give time as a double.
610  const auto make_lagrange_iterator = [](const auto& it) noexcept {
611  return boost::make_transform_iterator(
612  it, [](const Time& t) noexcept { return t.value(); });
613  };
614 
615  const auto union_remote_evaluation_step =
616  union_step(*remote_evaluation_step);
617  const auto union_step_lower_bound =
618  std::max(union_step_start, union_remote_evaluation_step);
619 
620  // Compute the contribution to an interpolation over the local
621  // times to `remote_evaluation_step->value()`, which we will
622  // use as the coupling value for that time. If there is an
623  // actual evaluation at that time then skip this because the
624  // Lagrange polynomial will be zero.
625  if (not std::binary_search(history.local_begin(), history.local_end(),
626  *remote_evaluation_step, less)) {
627  const auto union_step_upper_bound =
628  advance_within_step(union_remote_evaluation_step);
629  for (auto step = union_step_lower_bound;
630  step < union_step_upper_bound;
631  ++step) {
632  deriv_coef += base_summand(step, union_remote_evaluation_step);
633  }
634  deriv_coef *=
635  lagrange_polynomial(make_lagrange_iterator(local_evaluation_step),
636  remote_evaluation_step->value(),
637  make_lagrange_iterator(history.local_begin()),
638  make_lagrange_iterator(history.local_end()));
639  }
640 
641  // Same qualitative calculation as the previous block, but
642  // interpolating over the remote times. This case is somewhat
643  // more complicated because the latest remote time that can be
644  // used varies for the different segments making up the step.
645  if (not std::binary_search(remote_begin, history.remote_end(),
646  *local_evaluation_step, less)) {
647  auto union_step_upper_bound =
648  advance_within_step(union_local_evaluation_step);
649  if (history.remote_end() - remote_evaluation_step > order_s) {
650  union_step_upper_bound = std::min(
651  union_step_upper_bound,
652  union_step(*(remote_evaluation_step + order_s)));
653  }
654 
655  auto control_points = make_lagrange_iterator(
656  remote_evaluation_step - remote_begin >= order_s
657  ? remote_evaluation_step - (order_s - 1)
658  : remote_begin);
659  for (auto step = union_step_lower_bound;
660  step < union_step_upper_bound;
661  ++step, ++control_points) {
662  deriv_coef +=
663  base_summand(step, union_local_evaluation_step) *
665  make_lagrange_iterator(remote_evaluation_step),
666  local_evaluation_step->value(), control_points,
667  control_points +
668  static_cast<typename decltype(
669  control_points)::difference_type>(order_s));
670  }
671  }
672  }
673 
674  if (deriv_coef != 0.) {
675  // Skip the (potentially expensive) coupling calculation if
676  // the coefficient is zero.
677  accumulated_change +=
678  deriv_coef * history.coupling(coupling, local_evaluation_step,
679  remote_evaluation_step);
680  }
681  } // for remote_evaluation_step
682  } // for local_evaluation_step
683 
684  return accumulated_change;
685 }
686 
687 template <typename Vars, typename DerivVars>
688 bool AdamsBashforthN::can_change_step_size(
689  const TimeId& time_id,
690  const TimeSteppers::History<Vars, DerivVars>& history) const noexcept {
691  // We need to forbid local time-stepping before initialization is
692  // complete. The self-start procedure itself should never consider
693  // changing the step size, but we need to wait during the main
694  // evolution until the self-start history has been replaced with
695  // "real" values.
696  const evolution_less<Time> less{time_id.time_runs_forward()};
697  return history.size() == 0 or
698  (less(history.back(), time_id.time()) and
699  std::is_sorted(history.begin(), history.end(), less));
700 }
701 
702 template <typename Iterator, typename Delta>
703 std::vector<double> AdamsBashforthN::get_coefficients(
704  const Iterator& times_begin, const Iterator& times_end,
705  const Delta& step) noexcept {
706  ASSERT(times_begin != times_end, "No history provided");
707  std::vector<double> steps;
708  // This may be slightly more space than we need, but we can't get
709  // the exact amount without iterating through the iterators, which
710  // is not necessarily cheap depending on the iterator type.
711  steps.reserve(maximum_order);
712  for (auto t = times_begin; std::next(t) != times_end; ++t) {
713  steps.push_back((*std::next(t) - *t) / step);
714  }
715  steps.push_back(1.);
716  return get_coefficients_impl(steps);
717 }
718 } // namespace TimeSteppers
local_iterator local_begin() const noexcept
Access to the sequence of times on the indicated side.
Definition: BoundaryHistory.hpp:86
remote_iterator remote_begin() const noexcept
Access to the sequence of times on the indicated side.
Definition: BoundaryHistory.hpp:93
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:35
size_t remote_size() const noexcept
Access to the sequence of times on the indicated side.
Definition: BoundaryHistory.hpp:101
double lagrange_polynomial(const Iterator &index_point, double x, const Iterator &control_points_begin, const Iterator &control_points_end)
Evaluate the jth Lagrange interpolating polynomial with the given control points where j is the index...
Definition: LagrangePolynomial.hpp:16
Definition: AdamsBashforthN.hpp:55
Overloader< Fs... > make_overloader(Fs... fs)
Create Overloader<Fs...>, see Overloader for details.
Definition: Overloader.hpp:109
T signaling_NaN(T... args)
The time in a simulation. Times can be safely compared for exact equality as long as they do not belo...
Definition: Time.hpp:31
remote_iterator remote_end() const noexcept
Access to the sequence of times on the indicated side.
Definition: BoundaryHistory.hpp:96
A unique identifier for the temporal state of an integrated system.
Definition: TimeId.hpp:25
Defines classes and functions for making classes creatable from input files.
const Time & time() const noexcept
Time of the current substep.
Definition: TimeId.hpp:61
Defines macros to allow serialization of abstract template base classes.
Defines class TimeId.
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
size_t local_size() const noexcept
Access to the sequence of times on the indicated side.
Definition: BoundaryHistory.hpp:100
const_reference back() const noexcept
These return the past times. The other data can be accessed through HistoryIterator methods...
Definition: History.hpp:83
History data used by a TimeStepper.
Definition: History.hpp:25
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
Defines Time and TimeDelta.
double value() const noexcept
Approximate numerical value of the Time.
Definition: Time.hpp:51
bool is_positive() const noexcept
Test if the interval is oriented towards larger time.
Definition: Time.hpp:136
History data used by a TimeStepper for boundary integration.
Definition: BoundaryHistory.hpp:31
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:27
Defines class CachedFunction.
const_iterator begin() const noexcept
These iterators directly return the Time of the past values. The other data can be accessed through t...
Definition: History.hpp:61
local_iterator local_end() const noexcept
Access to the sequence of times on the indicated side.
Definition: BoundaryHistory.hpp:89
Represents an interval of time within a single slab.
Definition: Time.hpp:108
Defines functions lagrange_polynomial.
Define simple functions for constant expressions.
auto make_cached_function(Function function, PassedMapArgs... map_args) noexcept
Construct a CachedFunction wrapping the given function.
Definition: CachedFunction.hpp:53
Defines macro ASSERT.
Wraps the template metaprogramming library used (brigand)
const CouplingResult & coupling(Coupling &&c, const local_iterator &local, const remote_iterator &remote) const noexcept
Evaluate the coupling function at the given local and remote history entries. The coupling function w...
Definition: BoundaryHistory.hpp:151
Ordering functors that reverse their order when time runs backwards.
Definition: EvolutionOrdering.hpp:16
Defines functions and classes from the GSL.
const_iterator end() const noexcept
These iterators directly return the Time of the past values. The other data can be accessed through t...
Definition: History.hpp:66
std::result_of_t< const Coupling &(LocalVars, RemoteVars)> compute_boundary_delta(const Coupling &coupling, gsl::not_null< BoundaryHistoryType< LocalVars, RemoteVars, Coupling > *> history, const TimeDelta &time_step) const noexcept
Holds classes that take time steps.
Definition: BoundaryHistory.hpp:23
Defines macro ERROR.
double value() const noexcept
Approximate numerical length of the interval.
Definition: Time.hpp:131
An Nth Adams-Bashforth time stepper.
Definition: AdamsBashforthN.hpp:51
bool operator<(const Slab &a, const Slab &b) noexcept
Slab comparison operators give the time ordering. Overlapping unequal slabs should not be compared (a...
Definition: Slab.hpp:117
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
Defines make_with_value.