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