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