ryujin 2.1.1 revision 350e54cc11f3d21282bddcf3e6517944dbc508bf
time_integrator.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2022 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include "time_integrator.h"
9
10namespace ryujin
11{
12 using namespace dealii;
13
14 template <typename StateVector, typename Number>
15 void
16 sadd(StateVector &dst, const Number s, const Number b, const StateVector &src)
17 {
18 auto &dst_U = std::get<0>(dst);
19 auto &src_U = std::get<0>(src);
20 dst_U.sadd(s, b, src_U);
21
22 auto &dst_V = std::get<2>(dst);
23 auto &src_V = std::get<2>(src);
24 dst_V.sadd(s, b, src_V);
25 }
26
27
28 template <typename Description, int dim, typename Number>
29 TimeIntegrator<Description, dim, Number>::TimeIntegrator(
30 const MPIEnsemble &mpi_ensemble,
31 const OfflineData<dim, Number> &offline_data,
32 const HyperbolicModule<Description, dim, Number> &hyperbolic_module,
33 const ParabolicModule<Description, dim, Number> &parabolic_module,
34 const std::string &subsection /*= "TimeIntegrator"*/)
35 : ParameterAcceptor(subsection)
36 , mpi_ensemble_(mpi_ensemble)
37 , offline_data_(&offline_data)
38 , hyperbolic_module_(&hyperbolic_module)
39 , parabolic_module_(&parabolic_module)
40 {
41 cfl_min_ = Number(0.45);
42 add_parameter(
43 "cfl min",
44 cfl_min_,
45 "Minimal admissible relative CFL constant. How this parameter is used "
46 "depends on the chosen CFL recovery strategy");
47
48 cfl_max_ = Number(0.90);
49 add_parameter(
50 "cfl max",
51 cfl_max_,
52 "Maximal admissible relative CFL constant. How this parameter is used "
53 "depends on the chosen CFL recovery strategy");
54
55 cfl_recovery_strategy_ = CFLRecoveryStrategy::cruise_control;
56 add_parameter("cfl recovery strategy",
57 cfl_recovery_strategy_,
58 "CFL/invariant domain violation recovery strategy: none, "
59 "bang bang control, cruise control");
60
61 acceptable_tau_max_ratio_ = Number(2.0);
62 add_parameter("acceptable tau_max ratio",
63 acceptable_tau_max_ratio_,
64 "Maximal acceptable discrepancy between computed tau_max of "
65 "a (sub)step and enforced time-step size tau. If the ratio "
66 "is violated then a restart will be singnalled.");
67
68 if (ParabolicSystem::is_identity)
69 time_stepping_scheme_ = TimeSteppingScheme::erk_33;
70 else
71 time_stepping_scheme_ = TimeSteppingScheme::strang_erk_33_cn;
72 add_parameter("time stepping scheme",
73 time_stepping_scheme_,
74 "Time stepping scheme: ssprk 22, ssprk 33, erk 11, erk 22, "
75 "erk 33, erk 43, erk "
76 "54, strang ssprk 33 cn, strang erk 33 cn, strang erk 43 cn, "
77 "imex 11, imex 22, imex 33");
78 }
79
80
81 template <typename Description, int dim, typename Number>
82 void TimeIntegrator<Description, dim, Number>::prepare()
83 {
84#ifdef DEBUG_OUTPUT
85 std::cout << "TimeIntegrator<dim, Number>::prepare()" << std::endl;
86#endif
87
88 /* Resize temporary storage to appropriate sizes: */
89
90 switch (time_stepping_scheme_) {
91 case TimeSteppingScheme::ssprk_22:
92 temp_.resize(2);
93 efficiency_ = 1.;
94 break;
95 case TimeSteppingScheme::ssprk_33:
96 temp_.resize(2);
97 efficiency_ = 1.;
98 break;
99 case TimeSteppingScheme::erk_11:
100 temp_.resize(1);
101 efficiency_ = 1.;
102 break;
103 case TimeSteppingScheme::erk_22:
104 temp_.resize(2);
105 efficiency_ = 2.;
106 break;
107 case TimeSteppingScheme::erk_33:
108 temp_.resize(3);
109 efficiency_ = 3.;
110 break;
111 case TimeSteppingScheme::erk_43:
112 temp_.resize(4);
113 efficiency_ = 4.;
114 break;
115 case TimeSteppingScheme::erk_54:
116 temp_.resize(5);
117 efficiency_ = 5.;
118 break;
119 case TimeSteppingScheme::strang_ssprk_33_cn:
120 temp_.resize(3);
121 efficiency_ = 2.;
122 break;
123 case TimeSteppingScheme::strang_erk_33_cn:
124 temp_.resize(4);
125 efficiency_ = 6.;
126 break;
127 case TimeSteppingScheme::strang_erk_43_cn:
128 temp_.resize(4);
129 efficiency_ = 8.;
130 break;
131 case TimeSteppingScheme::imex_11:
132 temp_.resize(2);
133 efficiency_ = 1.;
134 break;
135 case TimeSteppingScheme::imex_22:
136 temp_.resize(4);
137 efficiency_ = 2.;
138 break;
139 case TimeSteppingScheme::imex_33:
140 temp_.resize(6);
141 efficiency_ = 3.;
142 break;
143 }
144
145 /* Initialize temporary vectors: */
146
147 for (auto &it : temp_) {
148 Vectors::reinit_state_vector<Description>(it, *offline_data_);
149 }
150
151 /* Reset CFL to starting value, set maximal acceptable tau_max ratio: */
152
153 AssertThrow(cfl_min_ > 0., ExcMessage("cfl min must be a positive value"));
154 AssertThrow(cfl_max_ >= cfl_min_,
155 ExcMessage("cfl max must be greater than or equal to cfl min"));
156
157 AssertThrow(
158 acceptable_tau_max_ratio_ >= 1.0,
159 ExcMessage(
160 "acceptable tau_max ratio must be greater than or equal to 1."));
161
162 hyperbolic_module_->cfl(cfl_max_);
163 hyperbolic_module_->acceptable_tau_max_ratio(acceptable_tau_max_ratio_);
164
165 const auto check_whether_timestepping_makes_sense = [&]() {
166 /*
167 * Make sure the user selects an appropriate time-stepping scheme.
168 */
169
170 switch (time_stepping_scheme_) {
171 case TimeSteppingScheme::ssprk_22:
172 [[fallthrough]];
173 case TimeSteppingScheme::ssprk_33:
174 [[fallthrough]];
175 case TimeSteppingScheme::erk_11:
176 [[fallthrough]];
177 case TimeSteppingScheme::erk_22:
178 [[fallthrough]];
179 case TimeSteppingScheme::erk_33:
180 [[fallthrough]];
181 case TimeSteppingScheme::erk_43:
182 [[fallthrough]];
183 case TimeSteppingScheme::erk_54: {
184 AssertThrow(
185 ParabolicSystem::is_identity,
186 dealii::ExcMessage(
187 "The selected equation consists of a hyperbolic and nontrivial "
188 "parabolic subsystem and requires an IMEX timestepping "
189 "scheme such as »strang erk 33 cn«."));
190 break;
191 }
192 case TimeSteppingScheme::imex_11:
193 [[fallthrough]];
194 case TimeSteppingScheme::imex_22:
195 [[fallthrough]];
196 case TimeSteppingScheme::imex_33:
197 [[fallthrough]];
198 case TimeSteppingScheme::strang_ssprk_33_cn:
199 [[fallthrough]];
200 case TimeSteppingScheme::strang_erk_33_cn:
201 [[fallthrough]];
202 case TimeSteppingScheme::strang_erk_43_cn: {
203 AssertThrow(
204 !ParabolicSystem::is_identity,
205 dealii::ExcMessage(
206 "The selected equation has a trivial parabolic subsystem and "
207 "should not be run with an IMEX timestepping scheme."));
208 break;
209 }
210 }
211 };
212
213 check_whether_timestepping_makes_sense();
214 this->parse_parameters_call_back.connect(
215 check_whether_timestepping_makes_sense);
216 }
217
218
219 /*
220 * -------------------------------------------------------------------------
221 * Prepare state vector:
222 * -------------------------------------------------------------------------
223 */
224
225
226 template <typename Description, int dim, typename Number>
227 void TimeIntegrator<Description, dim, Number>::prepare_state_vector(
228 StateVector &state_vector, Number t) const
229 {
230 if (!ParabolicSystem::is_identity)
231 parabolic_module_->prepare_state_vector(state_vector, t);
232 hyperbolic_module_->prepare_state_vector(state_vector, t);
233 }
234
235
236 /*
237 * -------------------------------------------------------------------------
238 * High level step function implementing various CFLRecoveryStrategy
239 * -------------------------------------------------------------------------
240 */
241
242
243 template <typename Description, int dim, typename Number>
244 Number TimeIntegrator<Description, dim, Number>::step(
245 StateVector &state_vector,
246 Number t,
247 Number t_final /*=std::numeric_limits<Number>::max()*/)
248 {
249 Number tau_max = t_final - t; /* enforces t <= t_final */
250
251#ifdef DEBUG_OUTPUT
252 std::cout << "TimeIntegrator<dim, Number>::step()" << std::endl;
253 std::cout << " enforcing tau_max <= " << tau_max << std::endl;
254#endif
255
256 const auto single_step = [&]() {
257 switch (time_stepping_scheme_) {
258 case TimeSteppingScheme::ssprk_22:
259 return step_ssprk_22(state_vector, t, tau_max);
260 case TimeSteppingScheme::ssprk_33:
261 return step_ssprk_33(state_vector, t, tau_max);
262 case TimeSteppingScheme::erk_11:
263 return step_erk_11(state_vector, t, tau_max);
264 case TimeSteppingScheme::erk_22:
265 return step_erk_22(state_vector, t, tau_max);
266 case TimeSteppingScheme::erk_33:
267 return step_erk_33(state_vector, t, tau_max);
268 case TimeSteppingScheme::erk_43:
269 return step_erk_43(state_vector, t, tau_max);
270 case TimeSteppingScheme::erk_54:
271 return step_erk_54(state_vector, t, tau_max);
272 case TimeSteppingScheme::strang_ssprk_33_cn:
273 return step_strang_ssprk_33_cn(state_vector, t, tau_max);
274 case TimeSteppingScheme::strang_erk_33_cn:
275 return step_strang_erk_33_cn(state_vector, t, tau_max);
276 case TimeSteppingScheme::strang_erk_43_cn:
277 return step_strang_erk_43_cn(state_vector, t, tau_max);
278 case TimeSteppingScheme::imex_11:
279 return step_imex_11(state_vector, t, tau_max);
280 case TimeSteppingScheme::imex_22:
281 return step_imex_22(state_vector, t, tau_max);
282 case TimeSteppingScheme::imex_33:
283 return step_imex_33(state_vector, t, tau_max);
284 default:
285 __builtin_unreachable();
286 }
287 };
288
289 if (cfl_recovery_strategy_ != CFLRecoveryStrategy::none) {
290 hyperbolic_module_->id_violation_strategy_ =
291 IDViolationStrategy::raise_exception;
292 parabolic_module_->id_violation_strategy_ =
293 IDViolationStrategy::raise_exception;
294 hyperbolic_module_->cfl(cfl_max_);
295 }
296
297 try {
298 return single_step();
299
300 } catch (const Restart &restart) {
301
302 AssertThrow(cfl_recovery_strategy_ != CFLRecoveryStrategy::none,
303 dealii::ExcInternalError());
304
305 hyperbolic_module_->id_violation_strategy_ = IDViolationStrategy::warn;
306 parabolic_module_->id_violation_strategy_ = IDViolationStrategy::warn;
307
308 if (cfl_recovery_strategy_ == CFLRecoveryStrategy::bang_bang_control) {
309 /* Retry with cfl_min instead of cfl_max: */
310#ifdef DEBUG_OUTPUT
311 std::cout
312 << " restart with bang bang control: setting cfl to cfl_min"
313 << std::endl;
314#endif
315 hyperbolic_module_->cfl(cfl_min_);
316 }
317
318 if (cfl_recovery_strategy_ == CFLRecoveryStrategy::cruise_control) {
319 /* Retry with the suggested tau_max: */
320#ifdef DEBUG_OUTPUT
321 std::cout
322 << " restart with cruise control: using suggested_tau_max"
323 << std::endl;
324#endif
325 //
326 // Multiply the suggested tau_max value with the efficiency.
327 //
328 // We have to account for the fact that the e Restart exception is
329 // thrown within a substep of the hyperbolic or parabolic module.
330 // This implies that the suggested_tau_max is computed for that
331 // particular substep and not for the full combined method (where
332 // tau_max can be larger). We thus multiply tau_max with the
333 // efficiency factor.
334 //
335 tau_max =
336 std::min(tau_max, efficiency_ * Number(restart.suggested_tau_max));
337 }
338
339 return single_step();
340 }
341 }
342
343
344 /*
345 * -------------------------------------------------------------------------
346 * Concrete implementation of ERK / IMEX time stepping strategies.
347 * -------------------------------------------------------------------------
348 */
349
350
351 template <typename Description, int dim, typename Number>
352 Number TimeIntegrator<Description, dim, Number>::step_ssprk_22(
353 StateVector &state_vector, Number t, Number tau_max)
354 {
355 /* SSP-RK2, see @cite Shu1988, Eq. 2.15. */
356
357#ifdef DEBUG_OUTPUT
358 std::cout << "TimeIntegrator<dim, Number>::step_ssprk_22()" << std::endl;
359#endif
360
361 Assert(efficiency_ == 1., dealii::ExcInternalError());
362
363 /* Step 1: T0 = U_old + tau * L(U_old) at t -> t + tau */
364 Number tau = hyperbolic_module_->template step<0>(
365 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
366
367 /* Step 2: T1 = T0 + tau L(T0) at time t + tau -> t + 2*tau */
368 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
369 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
370
371 /* Step 2: convex combination: T1 = 1/2 U_old + 1/2 T1 at time t + tau */
372 sadd(temp_[1], Number(1.0 / 2.0), Number(1.0 / 2.0), state_vector);
373
374 state_vector.swap(temp_[1]);
375 return tau;
376 }
377
378
379 template <typename Description, int dim, typename Number>
380 Number TimeIntegrator<Description, dim, Number>::step_ssprk_33(
381 StateVector &state_vector, Number t, Number tau_max)
382 {
383 /* SSP-RK3, see @cite Shu1988, Eq. 2.18. */
384
385#ifdef DEBUG_OUTPUT
386 std::cout << "TimeIntegrator<dim, Number>::step_ssprk_33()" << std::endl;
387#endif
388
389 Assert(efficiency_ == 1., dealii::ExcInternalError());
390
391 /* Step 1: T0 = U_old + tau * L(U_old) at time t -> t + tau */
392 Number tau = hyperbolic_module_->template step<0>(
393 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
394
395 /* Step 2: T1 = T0 + tau L(T0) at time t + tau -> t + 2*tau */
396 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
397 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
398
399 /* Step 2: convex combination T1 = 3/4 U_old + 1/4 T1 at time t + 0.5*tau */
400 sadd(temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), state_vector);
401
402 /* Step 3: T0 = T1 + tau L(T1) at time t + 0.5*tau -> t + 1.5*tau */
403 hyperbolic_module_->prepare_state_vector(temp_[1], t + 0.5 * tau);
404 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
405
406 /* Step 3: convex combination: T0 = 1/3 U_old + 2/3 T0 at time t + tau */
407 sadd(temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), state_vector);
408
409 state_vector.swap(temp_[0]);
410 return tau;
411 }
412
413
414 template <typename Description, int dim, typename Number>
415 Number TimeIntegrator<Description, dim, Number>::step_erk_11(
416 StateVector &state_vector, Number /*t*/, Number tau_max)
417 {
418#ifdef DEBUG_OUTPUT
419 std::cout << "TimeIntegrator<dim, Number>::step_erk_11()" << std::endl;
420#endif
421
422 Assert(efficiency_ == 1., dealii::ExcInternalError());
423
424 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
425 Number tau = hyperbolic_module_->template step<0>(
426 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
427
428 state_vector.swap(temp_[0]);
429 return tau;
430 }
431
432
433 template <typename Description, int dim, typename Number>
434 Number TimeIntegrator<Description, dim, Number>::step_erk_22(
435 StateVector &state_vector, Number t, Number tau_max)
436 {
437#ifdef DEBUG_OUTPUT
438 std::cout << "TimeIntegrator<dim, Number>::step_erk_22()" << std::endl;
439#endif
440
441 Assert(efficiency_ == 2., dealii::ExcInternalError());
442
443 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
444 Number tau = hyperbolic_module_->template step<0>(
445 state_vector, {}, {}, temp_[0], Number(.0), tau_max / efficiency_);
446
447 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + tau -> t + 2*tau */
448 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
449 hyperbolic_module_->template step<1>(
450 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
451
452 state_vector.swap(temp_[1]);
453 return efficiency_ * tau;
454 }
455
456
457 template <typename Description, int dim, typename Number>
458 Number TimeIntegrator<Description, dim, Number>::step_erk_33(
459 StateVector &state_vector, Number t, Number tau_max)
460 {
461#ifdef DEBUG_OUTPUT
462 std::cout << "TimeIntegrator<dim, Number>::step_erk_33()" << std::endl;
463#endif
464
465 Assert(efficiency_ == 3., dealii::ExcInternalError());
466
467 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
468 Number tau = hyperbolic_module_->template step<0>(
469 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
470
471 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + 1*tau -> t + 2*tau */
472 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
473 hyperbolic_module_->template step<1>(
474 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
475
476 /*
477 * Step 3: T2 <- {T1, 9/4} and {T0, -2} and {U_old, 3/4}
478 * at time t + 2*tau -> t + 3*tau
479 */
480 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
481 hyperbolic_module_->template step<2>(temp_[1],
482 {{state_vector, temp_[0]}},
483 {{Number(0.75), Number(-2.)}},
484 temp_[2],
485 tau);
486
487 state_vector.swap(temp_[2]);
488 return efficiency_ * tau;
489 }
490
491
492 template <typename Description, int dim, typename Number>
493 Number TimeIntegrator<Description, dim, Number>::step_erk_43(
494 StateVector &state_vector, Number t, Number tau_max)
495 {
496#ifdef DEBUG_OUTPUT
497 std::cout << "TimeIntegrator<dim, Number>::step_erk_43()" << std::endl;
498#endif
499
500 Assert(efficiency_ == 4., dealii::ExcInternalError());
501
502 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
503 Number tau = hyperbolic_module_->template step<0>(
504 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
505
506 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + 1*tau -> t + 2*tau */
507 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
508 hyperbolic_module_->template step<1>(
509 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
510
511 /* Step 3: T2 <- {T1, 2} and {T0, -1} at time t + 2*tau -> t + 3*tau */
512 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
513 hyperbolic_module_->template step<1>(
514 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
515
516 /*
517 * Step 4: T3 <- {T2, 8/3} and {T1,-10/3} and {T0, 5/3}
518 * at time t + 3*tau -> t + 4*tau
519 */
520 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
521 hyperbolic_module_->template step<2>(temp_[2],
522 {{temp_[0], temp_[1]}},
523 {{Number(5. / 3.), Number(-10. / 3.)}},
524 temp_[3],
525 tau);
526
527 state_vector.swap(temp_[3]);
528 return efficiency_ * tau;
529 }
530
531
532 template <typename Description, int dim, typename Number>
533 Number TimeIntegrator<Description, dim, Number>::step_erk_54(
534 StateVector &state_vector, Number t, Number tau_max)
535 {
536#ifdef DEBUG_OUTPUT
537 std::cout << "TimeIntegrator<dim, Number>::step_erk_54()" << std::endl;
538#endif
539
540 Assert(efficiency_ == 5., dealii::ExcInternalError());
541
542 constexpr Number c = 0.2; /* equidistant c_i */
543 constexpr Number a_21 = +0.2;
544 constexpr Number a_31 = +0.26075582269554909;
545 constexpr Number a_32 = +0.13924417730445096;
546 constexpr Number a_41 = -0.25856517872570289;
547 constexpr Number a_42 = +0.91136274166280729;
548 constexpr Number a_43 = -0.05279756293710430;
549 constexpr Number a_51 = +0.21623276431503774;
550 constexpr Number a_52 = +0.51534223099602405;
551 constexpr Number a_53 = -0.81662794199265554;
552 constexpr Number a_54 = +0.88505294668159373;
553 constexpr Number a_61 = -0.10511678454691901; /* aka b_1 */
554 constexpr Number a_62 = +0.87880047152100838; /* aka b_2 */
555 constexpr Number a_63 = -0.58903404061484477; /* aka b_3 */
556 constexpr Number a_64 = +0.46213380485434047; /* aka b_4 */
557 constexpr Number a_65 [[maybe_unused]] = +0.35321654878641495; /* aka b_5 */
558
559 /* Step 1: at time t -> t + 1*tau */
560 Number tau = hyperbolic_module_->template step<0>(
561 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
562
563 /* Step 2: at time t + 1*tau -> t + 2*tau */
564 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
565 hyperbolic_module_->template step<1>(
566 temp_[0], {{state_vector}}, {{(a_31 - a_21) / c}}, temp_[1], tau);
567
568 /* Step 3: at time t + 2*tau -> t + 3*tau */
569 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
570 hyperbolic_module_->template step<2>(
571 temp_[1],
572 {{state_vector, temp_[0]}},
573 {{(a_41 - a_31) / c, (a_42 - a_32) / c}},
574 temp_[2],
575 tau);
576
577 /* Step 4: at time t + 3*tau -> t + 4*tau */
578 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
579 hyperbolic_module_->template step<3>(
580 temp_[2],
581 {{state_vector, temp_[0], temp_[1]}},
582 {{(a_51 - a_41) / c, (a_52 - a_42) / c, (a_53 - a_43) / c}},
583 temp_[3],
584 tau);
585
586 /* Step 5: at time t + 4*tau -> t + 5*tau */
587 hyperbolic_module_->prepare_state_vector(temp_[3], t + 4.0 * tau);
588 hyperbolic_module_->template step<4>(
589 temp_[3],
590 {{state_vector, temp_[0], temp_[1], temp_[2]}},
591 {{(a_61 - a_51) / c,
592 (a_62 - a_52) / c,
593 (a_63 - a_53) / c,
594 (a_64 - a_54) / c}},
595 temp_[4],
596 tau);
597
598 state_vector.swap(temp_[4]);
599 return efficiency_ * tau;
600 }
601
602
603 template <typename Description, int dim, typename Number>
604 Number TimeIntegrator<Description, dim, Number>::step_strang_ssprk_33_cn(
605 StateVector &state_vector, Number t, Number tau_max)
606 {
607 // FIXME: avoid code duplication with step_ssprk_33
608
609#ifdef DEBUG_OUTPUT
610 std::cout << "TimeIntegrator<dim, Number>::step_strang_ssprk_33_cn()"
611 << std::endl;
612#endif
613
614 Assert(efficiency_ == 2., dealii::ExcInternalError());
615
616 /* First explicit SSPRK 3 step with final result in temp_[0]: */
617
618 Number tau = hyperbolic_module_->template step<0>(
619 state_vector, {}, {}, temp_[0], Number(0.0), tau_max / efficiency_);
620
621 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
622 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
623 sadd(temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), state_vector);
624
625 hyperbolic_module_->prepare_state_vector(temp_[1], t + 0.5 * tau);
626 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
627 sadd(temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), state_vector);
628
629 /* Implicit Crank-Nicolson step with final result in temp_[2]: */
630
631 try {
632 parabolic_module_->crank_nicolson_step(temp_[0], t, temp_[2], 2.0 * tau);
633 } catch (Restart &restart) {
634 /* Adjust suggested_tau_max. We multiply with efficiency_ again later */
635 restart.suggested_tau_max /= efficiency_;
636 throw;
637 }
638
639 /* Second SSPRK 3 step with final result in temp_[0]: */
640
641 hyperbolic_module_->prepare_state_vector( temp_[2], t + 1.0 * tau);
642 hyperbolic_module_->template step<0>( temp_[2], {}, {}, temp_[0], tau);
643
644 hyperbolic_module_->prepare_state_vector(temp_[0], t + 2.0 * tau);
645 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
646 sadd(temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), temp_[2]);
647
648 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.5 * tau);
649 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
650 sadd(temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), temp_[2]);
651
652 state_vector.swap(temp_[0]);
653 return efficiency_ * tau;
654 }
655
656
657 template <typename Description, int dim, typename Number>
658 Number TimeIntegrator<Description, dim, Number>::step_strang_erk_33_cn(
659 StateVector &state_vector, Number t, Number tau_max)
660 {
661 // FIXME: refactor to avoid code duplication with step_erk_33
662
663#ifdef DEBUG_OUTPUT
664 std::cout << "TimeIntegrator<dim, Number>::step_strang_erk_33_cn()"
665 << std::endl;
666#endif
667
668 Assert(efficiency_ == 6., dealii::ExcInternalError());
669
670 /* First explicit ERK(3,3,1) step with final result in temp_[2]: */
671
672 Number tau = hyperbolic_module_->template step<0>(
673 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
674
675 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
676 hyperbolic_module_->template step<1>(
677 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
678
679 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
680 hyperbolic_module_->template step<2>(temp_[1],
681 {{state_vector, temp_[0]}},
682 {{Number(0.75), Number(-2.)}},
683 temp_[2],
684 tau);
685
686 /* Implicit Crank-Nicolson step with final result in temp_[3]: */
687
688 try {
689 parabolic_module_->crank_nicolson_step(temp_[2], t, temp_[3], 6.0 * tau);
690 } catch (Restart &restart) {
691 /* Adjust suggested_tau_max. We multiply with efficiency_ again later */
692 restart.suggested_tau_max /= efficiency_;
693 throw;
694 }
695
696 /* Second explicit ERK(3,3,1) 3 step with final result in temp_[2]: */
697
698 hyperbolic_module_->prepare_state_vector(temp_[3], t + 3.0 * tau);
699 hyperbolic_module_->template step<0>(
700 temp_[3], {}, {}, temp_[0], tau);
701
702 hyperbolic_module_->prepare_state_vector(temp_[0], t + 4.0 * tau);
703 hyperbolic_module_->template step<1>(
704 temp_[0], {{ temp_[3]}}, {{Number(-1.)}}, temp_[1], tau);
705
706 hyperbolic_module_->prepare_state_vector(temp_[1], t + 5.0 * tau);
707 hyperbolic_module_->template step<2>(temp_[1],
708 {{ temp_[3], temp_[0]}},
709 {{Number(0.75), Number(-2.)}},
710 temp_[2],
711 tau);
712
713 state_vector.swap(temp_[2]);
714 return efficiency_ * tau;
715 }
716
717
718 template <typename Description, int dim, typename Number>
719 Number TimeIntegrator<Description, dim, Number>::step_strang_erk_43_cn(
720 StateVector &state_vector, Number t, Number tau_max)
721 {
722 // FIXME: refactor to avoid code duplication with step_erk_43
723
724#ifdef DEBUG_OUTPUT
725 std::cout << "TimeIntegrator<dim, Number>::step_strang_erk_43_cn()"
726 << std::endl;
727#endif
728
729 Assert(efficiency_ == 8., dealii::ExcInternalError());
730
731 /* First explicit ERK(4,3,1) step with final result in temp_[3]: */
732
733 Number tau = hyperbolic_module_->template step<0>(
734 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
735
736 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
737 hyperbolic_module_->template step<1>(
738 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
739
740 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
741 hyperbolic_module_->template step<1>(
742 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
743
744 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
745 hyperbolic_module_->template step<2>(temp_[2],
746 {{temp_[0], temp_[1]}},
747 {{Number(5. / 3.), Number(-10. / 3.)}},
748 temp_[3],
749 tau);
750
751 /* Implicit Crank-Nicolson step with final result in temp_[2]: */
752
753 try {
754 parabolic_module_->crank_nicolson_step(temp_[3], t, temp_[2], 8.0 * tau);
755 } catch (Restart &restart) {
756 /* Adjust suggested_tau_max. We multiply with efficiency_ again later */
757 restart.suggested_tau_max /= efficiency_;
758 throw;
759 }
760
761 /* Second explicit ERK(4,3,1) step with final result in temp_[3]: */
762
763 hyperbolic_module_->prepare_state_vector(temp_[2], t + 4.0 * tau);
764 hyperbolic_module_->template step<0>(
765 temp_[2], {}, {}, temp_[0], tau);
766
767 hyperbolic_module_->prepare_state_vector(temp_[0], t + 5.0 * tau);
768 hyperbolic_module_->template step<1>(
769 temp_[0], {{ temp_[2]}}, {{Number(-1.)}}, temp_[1], tau);
770
771 hyperbolic_module_->prepare_state_vector(temp_[1], t + 6.0 * tau);
772 hyperbolic_module_->template step<1>(
773 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
774
775 hyperbolic_module_->prepare_state_vector(temp_[2], t + 7.0 * tau);
776 hyperbolic_module_->template step<2>(temp_[2],
777 {{temp_[0], temp_[1]}},
778 {{Number(5. / 3.), Number(-10. / 3.)}},
779 temp_[3],
780 tau);
781
782 state_vector.swap(temp_[3]);
783 return efficiency_ * tau;
784 }
785
786
787 template <typename Description, int dim, typename Number>
788 Number TimeIntegrator<Description, dim, Number>::step_imex_11(
789 StateVector &state_vector, Number t, Number tau_max)
790 {
791#ifdef DEBUG_OUTPUT
792 std::cout << "TimeIntegrator<dim, Number>::step_imex_11()" << std::endl;
793#endif
794
795 Assert(efficiency_ == 1., dealii::ExcInternalError());
796
797 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
798 Number tau = hyperbolic_module_->template step<0>(
799 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
800
801 /* Implicit step 1: T1 <- {T0, 1} at time t -> t + tau */
802 parabolic_module_->template backward_euler_step<0>(
803 temp_[0], t, {}, {}, temp_[1], 1.0 * tau);
804
805 state_vector.swap(temp_[1]);
806 return tau;
807 }
808
809
810 template <typename Description, int dim, typename Number>
811 Number TimeIntegrator<Description, dim, Number>::step_imex_22(
812 StateVector &state_vector, Number t, Number tau_max)
813 {
814#ifdef DEBUG_OUTPUT
815 std::cout << "TimeIntegrator<dim, Number>::step_imex_22()" << std::endl;
816#endif
817
818 Assert(efficiency_ == 2., dealii::ExcInternalError());
819
820 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
821 Number tau = hyperbolic_module_->template step<0>(
822 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
823
824 /* Implicit step 1: T1 <- {T0, 1} at time t -> t + tau */
825 parabolic_module_->template backward_euler_step<0>(
826 temp_[0], t, {}, {}, temp_[1], tau);
827
828 /* Explicit step 2: T2 <- {T1, 2} and {U_old, -1} at t + tau -> t + 2 tau */
829 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.0 * tau);
830 hyperbolic_module_->template step<1>(
831 temp_[1], {{state_vector}}, {{Number(-1.)}}, temp_[2], tau);
832
833 /* Implicit step 2: T3 <- {T2, 0} and {U_old, 1} at t + tau -> t + 2 tau */
834 parabolic_module_->template backward_euler_step<1>(temp_[2],
835 t + 1.0 * tau,
836 {{state_vector}},
837 {{Number(1.)}},
838 temp_[3],
839 tau);
840
841 state_vector.swap(temp_[3]);
842 return efficiency_ * tau;
843 }
844
845
846 template <typename Description, int dim, typename Number>
847 Number TimeIntegrator<Description, dim, Number>::step_imex_33(
848 StateVector &state_vector, Number t, Number tau_max)
849 {
850#ifdef DEBUG_OUTPUT
851 std::cout << "TimeIntegrator<dim, Number>::step_imex_33()" << std::endl;
852#endif
853
854 Assert(efficiency_ == 3., dealii::ExcInternalError());
855
856 /* IMEX(3, 3; 1), see @cite ErnGuermond2023, Sec. 4.3. */
857
858 const Number gamma = Number(0.5) + std::sqrt(Number(3.0)) / Number(6.0);
859
860 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
861 Number tau = hyperbolic_module_->template step<0>(
862 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
863
864 /* Implicit step 1: T1 <- {U_old, 1 - 3*gamma} at time t -> t + tau */
865 parabolic_module_->template backward_euler_step<1>(
866 temp_[0],
867 t,
868 {{state_vector}},
869 {{Number(1. - 3. * gamma)}},
870 temp_[1],
871 tau);
872
873 /* Explicit step 2: T2 <- {U_old, -1} and {T1, 2} at time t -> t + 2 tau */
874 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.0 * tau);
875 hyperbolic_module_->template step<1>(
876 temp_[1], {{state_vector}}, {{Number(-1.)}}, temp_[2], tau);
877
878 /*
879 * Implicit step 2:
880 * T3 <- {U_old, 6*gamma-1} and {T1, 2-9*gamma} at t -> t + * 2 tau
881 */
882 parabolic_module_->template backward_euler_step<2>(
883 temp_[2],
884 t + tau,
885 {{state_vector, temp_[1]}},
886 {{Number(6. * gamma - 1.), Number(2. - 9 * gamma)}},
887 temp_[3],
888 tau);
889
890 /* Explicit step 3: T4 <- {U_old, 3 / 4} and {T1, -2} at t -> t + 3 tau */
891 hyperbolic_module_->prepare_state_vector(temp_[3], t + 2. * tau);
892 hyperbolic_module_->template step<2>(temp_[3],
893 {{state_vector, temp_[1]}},
894 {{Number(0.75), Number(-2.)}},
895 temp_[4],
896 tau);
897
898 /* Implicit step 3: */
899 parabolic_module_->template backward_euler_step<3>(
900 temp_[4],
901 t + 2. * tau,
902 {{state_vector, temp_[1], temp_[3]}},
903 {{Number(0.75 - 3. * gamma),
904 Number(6. * gamma - 2.),
905 Number(9. / 4. - 3. * gamma)}},
906 temp_[5],
907 tau);
908
909 state_vector.swap(temp_[5]);
910 return efficiency_ * tau;
911 }
912
913} /* namespace ryujin */
ryujin::TimeIntegrator::step_imex_33
Number step_imex_33(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:847
ryujin::TimeIntegrator::TimeIntegrator
TimeIntegrator(const MPIEnsemble &mpi_ensemble, const OfflineData< dim, Number > &offline_data, const HyperbolicModule< Description, dim, Number > &hyperbolic_module, const ParabolicModule< Description, dim, Number > &parabolic_module, const std::string &subsection="/TimeIntegrator")
Definition: time_integrator.template.h:29
ryujin::TimeIntegrator::step_strang_ssprk_33_cn
Number step_strang_ssprk_33_cn(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:604
ryujin::TimeIntegrator::step_erk_22
Number step_erk_22(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:434
ryujin::TimeIntegrator::step_imex_11
Number step_imex_11(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:788
ryujin::TimeIntegrator::step_strang_erk_43_cn
Number step_strang_erk_43_cn(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:719
ryujin::TimeIntegrator::step_erk_54
Number step_erk_54(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:533
ryujin::TimeIntegrator::step_erk_33
Number step_erk_33(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:458
ryujin::TimeIntegrator::step_erk_43
Number step_erk_43(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:493
ryujin::TimeIntegrator::step_imex_22
Number step_imex_22(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:811
ryujin::TimeIntegrator::StateVector
typename View::StateVector StateVector
Definition: time_integrator.h:267
ryujin::TimeIntegrator::step_erk_11
Number step_erk_11(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:415
ryujin::TimeIntegrator::step_ssprk_22
Number step_ssprk_22(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:352
ryujin::TimeIntegrator::step_ssprk_33
Number step_ssprk_33(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:380
ryujin::TimeIntegrator::prepare_state_vector
void prepare_state_vector(StateVector &state_vector, Number t) const
Definition: time_integrator.template.h:227
ryujin::TimeIntegrator::step
Number step(StateVector &state_vector, Number t, Number t_final=std::numeric_limits< Number >::max())
Definition: time_integrator.template.h:244
ryujin::TimeIntegrator::step_strang_erk_33_cn
Number step_strang_erk_33_cn(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:658
ryujin::IDViolationStrategy::warn
@ warn
ryujin::IDViolationStrategy::raise_exception
@ raise_exception
ryujin::TimeSteppingScheme::erk_22
@ erk_22
ryujin::TimeSteppingScheme::erk_33
@ erk_33
ryujin::TimeSteppingScheme::imex_11
@ imex_11
ryujin::TimeSteppingScheme::erk_43
@ erk_43
ryujin::TimeSteppingScheme::erk_11
@ erk_11
ryujin::TimeSteppingScheme::ssprk_33
@ ssprk_33
ryujin::TimeSteppingScheme::strang_erk_43_cn
@ strang_erk_43_cn
ryujin::TimeSteppingScheme::ssprk_22
@ ssprk_22
ryujin::TimeSteppingScheme::strang_erk_33_cn
@ strang_erk_33_cn
ryujin::TimeSteppingScheme::imex_22
@ imex_22
ryujin::TimeSteppingScheme::strang_ssprk_33_cn
@ strang_ssprk_33_cn
ryujin::TimeSteppingScheme::imex_33
@ imex_33
ryujin::TimeSteppingScheme::erk_54
@ erk_54
ryujin::CFLRecoveryStrategy::bang_bang_control
@ bang_bang_control
ryujin::CFLRecoveryStrategy::none
@ none
ryujin::CFLRecoveryStrategy::cruise_control
@ cruise_control
ryujin::Vectors::StateVector
std::tuple< MultiComponentVector< Number, problem_dim >, MultiComponentVector< Number, prec_dim >, BlockVector< Number > > StateVector
Definition: state_vector.h:51
ryujin::sadd
void sadd(StateVector &dst, const Number s, const Number b, const StateVector &src)
Definition: time_integrator.template.h:16
time_integrator.h
Generated by doxygen 1.9.4