SpECTRE Documentation Coverage Report
Current view: top level - Time/TimeSteppers - AdamsBashforthN.hpp Hit Total Coverage
Commit: b1342d46f40e2d46bbd11d0cef68fd973031a24b Lines: 4 54 7.4 %
Date: 2020-09-24 20:24:42
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14