Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <cstdint>
8 :
9 : #include "Options/String.hpp"
10 : #include "Time/TimeStepId.hpp"
11 : #include "Time/TimeSteppers/LtsTimeStepper.hpp"
12 : #include "Utilities/Serialization/CharmPupable.hpp"
13 : #include "Utilities/TMPL.hpp"
14 :
15 : /// \cond
16 : class TimeDelta;
17 : namespace PUP {
18 : class er;
19 : } // namespace PUP
20 : namespace TimeSteppers {
21 : template <typename T>
22 : class BoundaryHistoryEvaluator;
23 : class ConstBoundaryHistoryTimes;
24 : template <typename T>
25 : class ConstUntypedHistory;
26 : class MutableBoundaryHistoryTimes;
27 : template <typename T>
28 : class MutableUntypedHistory;
29 : } // namespace TimeSteppers
30 : namespace gsl {
31 : template <class T>
32 : class not_null;
33 : } // namespace gsl
34 : /// \endcond
35 :
36 : namespace TimeSteppers {
37 :
38 : /*!
39 : * \ingroup TimeSteppersGroup
40 : *
41 : * An Nth order Adams-Bashforth time stepper.
42 : *
43 : * The stable step size factors for different orders are given by:
44 : *
45 : * <table class="doxtable">
46 : * <tr>
47 : * <th> %Order </th>
48 : * <th> CFL Factor </th>
49 : * </tr>
50 : * <tr>
51 : * <td> 1 </td>
52 : * <td> 1 </td>
53 : * </tr>
54 : * <tr>
55 : * <td> 2 </td>
56 : * <td> 1 / 2 </td>
57 : * </tr>
58 : * <tr>
59 : * <td> 3 </td>
60 : * <td> 3 / 11 </td>
61 : * </tr>
62 : * <tr>
63 : * <td> 4 </td>
64 : * <td> 3 / 20 </td>
65 : * </tr>
66 : * <tr>
67 : * <td> 5 </td>
68 : * <td> 45 / 551 </td>
69 : * </tr>
70 : * <tr>
71 : * <td> 6 </td>
72 : * <td> 5 / 114 </td>
73 : * </tr>
74 : * <tr>
75 : * <td> 7 </td>
76 : * <td> 945 / 40663 </td>
77 : * </tr>
78 : * <tr>
79 : * <td> 8 </td>
80 : * <td> 945 / 77432 </td>
81 : * </tr>
82 : * </table>
83 : *
84 : * \section lts Local time stepping calculation
85 : *
86 : * \f$\newcommand\tL{t^L}\newcommand\tR{t^R}\newcommand\tU{\tilde{t}\!}
87 : * \newcommand\mat{\mathbf}\f$
88 : *
89 : * Suppose the local and remote sides of the interface are evaluated
90 : * at times \f$\ldots, \tL_{-1}, \tL_0, \tL_1, \ldots\f$ and
91 : * \f$\ldots, \tR_{-1}, \tR_0, \tR_1, \ldots\f$, respectively, with
92 : * the starting location of the numbering arbitrary in each case.
93 : * Let the step we wish to calculate the effect of be the step from
94 : * \f$\tL_{m_S}\f$ to \f$\tL_{m_S+1}\f$. We call the sequence
95 : * produced from the union of the local and remote time sequences
96 : * \f$\ldots, \tU_{-1}, \tU_0, \tU_1, \ldots\f$. For example, one
97 : * possible sequence of times is:
98 : * \f{equation}
99 : * \begin{aligned}
100 : * \text{Local side:} \\ \text{Union times:} \\ \text{Remote side:}
101 : * \end{aligned}
102 : * \cdots
103 : * \begin{gathered}
104 : * \, \\ \tU_1 \\ \tR_5
105 : * \end{gathered}
106 : * \leftarrow \Delta \tU_1 \rightarrow
107 : * \begin{gathered}
108 : * \tL_4 \\ \tU_2 \\ \,
109 : * \end{gathered}
110 : * \leftarrow \Delta \tU_2 \rightarrow
111 : * \begin{gathered}
112 : * \, \\ \tU_3 \\ \tR_6
113 : * \end{gathered}
114 : * \leftarrow \Delta \tU_3 \rightarrow
115 : * \begin{gathered}
116 : * \, \\ \tU_4 \\ \tR_7
117 : * \end{gathered}
118 : * \leftarrow \Delta \tU_4 \rightarrow
119 : * \begin{gathered}
120 : * \tL_5 \\ \tU_5 \\ \,
121 : * \end{gathered}
122 : * \cdots
123 : * \f}
124 : * We call the indices of the step's start and end times in the
125 : * union time sequence \f$n_S\f$ and \f$n_E\f$, respectively. We
126 : * define \f$n^L_m\f$ to be the union-time index corresponding to
127 : * \f$\tL_m\f$ and \f$m^L_n\f$ to be the index of the last local
128 : * time not later than \f$\tU_n\f$ and similarly for the remote
129 : * side. So for the above example, \f$n^L_4 = 2\f$ and \f$m^R_2 =
130 : * 5\f$, and if we wish to compute the step from \f$\tL_4\f$ to
131 : * \f$\tL_5\f$ we would have \f$m_S = 4\f$, \f$n_S = 2\f$, and
132 : * \f$n_E = 5\f$.
133 : *
134 : * If we wish to evaluate the change over this step to \f$k\f$th
135 : * order, we can write the change in the value as a linear
136 : * combination of the values of the coupling between the elements at
137 : * unequal times:
138 : * \f{equation}
139 : * \mat{F}_{m_S} =
140 : * \mspace{-10mu}
141 : * \sum_{q^L = m_S-(k-1)}^{m_S}
142 : * \,
143 : * \sum_{q^R = m^R_{n_S}-(k-1)}^{m^R_{n_E-1}}
144 : * \mspace{-10mu}
145 : * \mat{D}_{q^Lq^R}
146 : * I_{q^Lq^R},
147 : * \f}
148 : * where \f$\mat{D}_{q^Lq^R}\f$ is the coupling function evaluated
149 : * between data from \f$\tL_{q^L}\f$ and \f$\tR_{q^R}\f$. The
150 : * coefficients can be written as the sum of three terms,
151 : * \f{equation}
152 : * I_{q^Lq^R} = I^E_{q^Lq^R} + I^R_{q^Lq^R} + I^L_{q^Lq^R},
153 : * \f}
154 : * which can be interpreted as a contribution from equal-time
155 : * evaluations and contributions related to the remote and local
156 : * evaluation times. These are given by
157 : * \f{align}
158 : * I^E_{q^Lq^R} &=
159 : * \mspace{-10mu}
160 : * \sum_{n=n_S}^{\min\left\{n_E, n^L+k\right\}-1}
161 : * \mspace{-10mu}
162 : * \tilde{\alpha}_{n,n-n^L} \Delta \tU_n
163 : * &&\text{if $\tL_{q^L} = \tR_{q^R}$, otherwise 0}
164 : * \\
165 : * I^R_{q^Lq^R} &=
166 : * \ell_{q^L - m_S + k}\!\left(
167 : * \tU_{n^R}; \tL_{m_S - (k-1)}, \ldots, \tL_{m_S}\right)
168 : * \mspace{-10mu}
169 : * \sum_{n=\max\left\{n_S, n^R\right\}}
170 : * ^{\min\left\{n_E, n^R+k\right\}-1}
171 : * \mspace{-10mu}
172 : * \tilde{\alpha}_{n,n-n^R} \Delta \tU_n
173 : * &&\text{if $\tR_{q^R}$ is not in $\{\tL_{\vphantom{|}\cdots}\}$,
174 : * otherwise 0}
175 : * \\
176 : * I^L_{q^Lq^R} &=
177 : * \mspace{-10mu}
178 : * \sum_{n=\max\left\{n_S, n^R\right\}}
179 : * ^{\min\left\{n_E, n^L+k, n^R_{q^R+k}\right\}-1}
180 : * \mspace{-10mu}
181 : * \ell_{q^R - m^R_n + k}\!\left(\tU_{n^L};
182 : * \tR_{m^R_n - (k-1)}, \ldots, \tR_{m^R_n}\right)
183 : * \tilde{\alpha}_{n,n-n^L} \Delta \tU_n
184 : * &&\text{if $\tL_{q^L}$ is not in $\{\tR_{\vphantom{|}\cdots}\}$,
185 : * otherwise 0,}
186 : * \f}
187 : * where for brevity we write \f$n^L = n^L_{q^L}\f$ and \f$n^R =
188 : * n^R_{q^R}\f$, and where \f$\ell_a(t; x_1, \ldots, x_k)\f$ a
189 : * Lagrange interpolating polynomial and \f$\tilde{\alpha}_{nj}\f$
190 : * is the \f$j\f$th coefficient for an Adams-Bashforth step over the
191 : * union times from step \f$n\f$ to step \f$n+1\f$.
192 : */
193 1 : class AdamsBashforth : public LtsTimeStepper {
194 : public:
195 0 : static constexpr const size_t minimum_order = 1;
196 0 : static constexpr const size_t maximum_order = 8;
197 :
198 0 : struct Order {
199 0 : using type = size_t;
200 0 : static constexpr Options::String help = {"Convergence order"};
201 0 : static type lower_bound() { return 1; }
202 0 : static type upper_bound() { return maximum_order; }
203 : };
204 0 : using options = tmpl::list<Order>;
205 0 : static constexpr Options::String help = {
206 : "An Adams-Bashforth Nth order time-stepper."};
207 :
208 0 : AdamsBashforth() = default;
209 0 : explicit AdamsBashforth(size_t order);
210 0 : AdamsBashforth(const AdamsBashforth&) = default;
211 0 : AdamsBashforth& operator=(const AdamsBashforth&) = default;
212 0 : AdamsBashforth(AdamsBashforth&&) = default;
213 0 : AdamsBashforth& operator=(AdamsBashforth&&) = default;
214 0 : ~AdamsBashforth() override = default;
215 :
216 1 : size_t order() const override;
217 :
218 1 : size_t error_estimate_order() const override;
219 :
220 1 : uint64_t number_of_substeps() const override;
221 :
222 1 : uint64_t number_of_substeps_for_error() const override;
223 :
224 1 : size_t number_of_past_steps() const override;
225 :
226 1 : double stable_step() const override;
227 :
228 1 : bool monotonic() const override;
229 :
230 1 : TimeStepId next_time_id(const TimeStepId& current_id,
231 : const TimeDelta& time_step) const override;
232 :
233 1 : TimeStepId next_time_id_for_error(const TimeStepId& current_id,
234 : const TimeDelta& time_step) const override;
235 :
236 1 : bool neighbor_data_required(
237 : const TimeStepId& next_substep_id,
238 : const TimeStepId& neighbor_data_id) const override;
239 :
240 1 : bool neighbor_data_required(
241 : double dense_output_time,
242 : const TimeStepId& neighbor_data_id) const override;
243 :
244 0 : WRAPPED_PUPable_decl_template(AdamsBashforth); // NOLINT
245 :
246 0 : explicit AdamsBashforth(CkMigrateMessage* /*unused*/) {}
247 :
248 : // clang-tidy: do not pass by non-const reference
249 0 : void pup(PUP::er& p) override; // NOLINT
250 :
251 : private:
252 0 : friend bool operator==(const AdamsBashforth& lhs, const AdamsBashforth& rhs);
253 :
254 : // Some of the private methods take a parameter of type "Delta" or
255 : // "TimeType". Delta is expected to be a TimeDelta or an
256 : // ApproximateTimeDelta, and TimeType is expected to be a Time or an
257 : // ApproximateTime. The former cases will detect and optimize the
258 : // constant-time-step case, while the latter are necessary for dense
259 : // output.
260 : template <typename T>
261 0 : void update_u_impl(gsl::not_null<T*> u,
262 : const MutableUntypedHistory<T>& history,
263 : const TimeDelta& time_step) const;
264 :
265 : template <typename T>
266 0 : bool update_u_impl(gsl::not_null<T*> u, gsl::not_null<T*> u_error,
267 : const MutableUntypedHistory<T>& history,
268 : const TimeDelta& time_step) const;
269 :
270 : template <typename T>
271 0 : bool dense_update_u_impl(gsl::not_null<T*> u,
272 : const ConstUntypedHistory<T>& history,
273 : double time) const;
274 :
275 : template <typename T, typename Delta>
276 0 : void update_u_common(gsl::not_null<T*> u,
277 : const ConstUntypedHistory<T>& history,
278 : const Delta& time_step, size_t order) const;
279 :
280 : template <typename T>
281 0 : bool can_change_step_size_impl(const TimeStepId& time_id,
282 : const ConstUntypedHistory<T>& history) const;
283 :
284 : template <typename T>
285 0 : void add_boundary_delta_impl(
286 : gsl::not_null<T*> result,
287 : const TimeSteppers::MutableBoundaryHistoryTimes& local_times,
288 : const TimeSteppers::MutableBoundaryHistoryTimes& remote_times,
289 : const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
290 : const TimeDelta& time_step) const;
291 :
292 : template <typename T>
293 0 : void boundary_dense_output_impl(
294 : gsl::not_null<T*> result,
295 : const TimeSteppers::ConstBoundaryHistoryTimes& local_times,
296 : const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
297 : const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
298 : const double time) const;
299 :
300 : template <typename T, typename TimeType>
301 0 : void boundary_impl(gsl::not_null<T*> result,
302 : const ConstBoundaryHistoryTimes& local_times,
303 : const ConstBoundaryHistoryTimes& remote_times,
304 : const BoundaryHistoryEvaluator<T>& coupling,
305 : const TimeType& end_time) const;
306 :
307 : TIME_STEPPER_DECLARE_OVERLOADS
308 : LTS_TIME_STEPPER_DECLARE_OVERLOADS
309 :
310 0 : size_t order_ = 3;
311 : };
312 :
313 0 : bool operator!=(const AdamsBashforth& lhs, const AdamsBashforth& rhs);
314 : } // namespace TimeSteppers
|
