ryujin 2.1.1 revision e1e0c9cf84519540dcc04ed5bb49dc8a8c052779
Loading...
Searching...
No Matches
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 - 2026 by the ryujin authors
4//
5
6#pragma once
7
8#include "time_integrator.h"
9
10#include <utility>
11
12namespace ryujin
13{
14 using namespace dealii;
15
16 template <typename HD, typename PD, int dim, typename Number>
17 TimeIntegrator<HD, PD, dim, Number>::TimeIntegrator(
18 const MPIEnsemble &mpi_ensemble,
19 const OfflineData<dim, Number> &offline_data,
20 const HyperbolicModule<HD, dim, Number> &hyperbolic_module,
21 const ParabolicModule<HD, PD, dim, Number> &parabolic_module,
22 const std::string &subsection /*= "TimeIntegrator"*/)
23 : ParameterAcceptor(subsection)
24 , mpi_ensemble_(mpi_ensemble)
25 , offline_data_(&offline_data)
26 , hyperbolic_module_(&hyperbolic_module)
27 , parabolic_module_(&parabolic_module)
28 {
29 cfl_min_ = Number(0.45);
30 add_parameter(
31 "cfl min",
32 cfl_min_,
33 "Minimal admissible relative CFL constant. How this parameter is used "
34 "depends on the chosen CFL recovery strategy");
35
36 cfl_max_ = Number(0.90);
37 add_parameter(
38 "cfl max",
39 cfl_max_,
40 "Maximal admissible relative CFL constant. How this parameter is used "
41 "depends on the chosen CFL recovery strategy");
42
43 cfl_recovery_strategy_ = CFLRecoveryStrategy::cruise_control;
44 add_parameter("cfl recovery strategy",
45 cfl_recovery_strategy_,
46 "CFL/invariant domain violation recovery strategy: none, "
47 "bang bang control, cruise control");
48
49 acceptable_tau_max_ratio_ = Number(2.0);
50 add_parameter("acceptable tau_max ratio",
51 acceptable_tau_max_ratio_,
52 "Maximal acceptable discrepancy between computed tau_max of "
53 "a (sub)step and enforced time-step size tau. If the ratio "
54 "is violated then a restart will be singnalled.");
55
56 tau_max_ = std::numeric_limits<Number>::max();
57 add_parameter("tau_max", tau_max_, "Largest time step size allowed.");
58
59 if (ParabolicSystem::is_identity)
60 time_stepping_scheme_ = TimeSteppingScheme::erk_33;
61 else
62 time_stepping_scheme_ = TimeSteppingScheme::strang_erk_33_cn;
63 add_parameter("time stepping scheme",
64 time_stepping_scheme_,
65 "Time stepping scheme: ssprk 22, ssprk 33, erk 11, erk 22, "
66 "erk 33, erk 43, erk "
67 "54, strang ssprk 33 cn, strang erk 33 cn, strang erk 43 cn, "
68 "imex 11, imex 22, imex 33");
69 }
70
71
72 template <typename HD, typename PD, int dim, typename Number>
73 void TimeIntegrator<HD, PD, dim, Number>::prepare()
74 {
75#ifdef DEBUG_OUTPUT
76 std::cout << "TimeIntegrator<dim, Number>::prepare()" << std::endl;
77#endif
78
79 /* Resize temporary storage to appropriate sizes: */
80
81 switch (time_stepping_scheme_) {
82 case TimeSteppingScheme::ssprk_22:
83 temp_.resize(2);
84 efficiency_ = 1.;
85 break;
86 case TimeSteppingScheme::ssprk_33:
87 temp_.resize(2);
88 efficiency_ = 1.;
89 break;
90 case TimeSteppingScheme::erk_11:
91 temp_.resize(1);
92 efficiency_ = 1.;
93 break;
94 case TimeSteppingScheme::erk_22:
95 temp_.resize(2);
96 efficiency_ = 2.;
97 break;
98 case TimeSteppingScheme::erk_33:
99 temp_.resize(3);
100 efficiency_ = 3.;
101 break;
102 case TimeSteppingScheme::erk_43:
103 temp_.resize(4);
104 efficiency_ = 4.;
105 break;
106 case TimeSteppingScheme::erk_54:
107 temp_.resize(5);
108 efficiency_ = 5.;
109 break;
110 case TimeSteppingScheme::strang_ssprk_33_cn:
111 temp_.resize(3);
112 efficiency_ = 2.;
113 break;
114 case TimeSteppingScheme::strang_erk_33_cn:
115 temp_.resize(4);
116 efficiency_ = 6.;
117 break;
118 case TimeSteppingScheme::strang_erk_43_cn:
119 temp_.resize(4);
120 efficiency_ = 8.;
121 break;
122 case TimeSteppingScheme::imex_11:
123 temp_.resize(2);
124 efficiency_ = 1.;
125 break;
126 case TimeSteppingScheme::imex_22:
127 temp_.resize(4);
128 efficiency_ = 2.;
129 break;
130 case TimeSteppingScheme::imex_33:
131 temp_.resize(6);
132 efficiency_ = 3.;
133 break;
134 }
135
136 /* Initialize temporary vectors: */
137
138 for (auto &it : temp_) {
139 hyperbolic_module_->reinit_state_vector(it);
140 parabolic_module_->reinit_state_vector(it);
141 }
142
143 /* Reset CFL to starting value, set maximal acceptable tau_max ratio: */
144
145 AssertThrow(cfl_min_ > 0., ExcMessage("cfl min must be a positive value"));
146 AssertThrow(cfl_max_ >= cfl_min_,
147 ExcMessage("cfl max must be greater than or equal to cfl min"));
148
149 AssertThrow(
150 acceptable_tau_max_ratio_ >= 1.0,
151 ExcMessage(
152 "acceptable tau_max ratio must be greater than or equal to 1."));
153
154 hyperbolic_module_->set_cfl(cfl_max_);
155 hyperbolic_module_->set_acceptable_tau_max_ratio(acceptable_tau_max_ratio_);
156
157 const auto check_whether_timestepping_makes_sense = [&]() {
158 /*
159 * Make sure the user selects an appropriate time-stepping scheme.
160 */
161
162 switch (time_stepping_scheme_) {
163 case TimeSteppingScheme::ssprk_22:
164 [[fallthrough]];
165 case TimeSteppingScheme::ssprk_33:
166 [[fallthrough]];
167 case TimeSteppingScheme::erk_11:
168 [[fallthrough]];
169 case TimeSteppingScheme::erk_22:
170 [[fallthrough]];
171 case TimeSteppingScheme::erk_33:
172 [[fallthrough]];
173 case TimeSteppingScheme::erk_43:
174 [[fallthrough]];
175 case TimeSteppingScheme::erk_54: {
176 AssertThrow(
177 ParabolicSystem::is_identity,
178 dealii::ExcMessage(
179 "The selected equation consists of a hyperbolic and nontrivial "
180 "parabolic subsystem and requires an IMEX timestepping "
181 "scheme such as »strang erk 33 cn«."));
182 break;
183 }
184 case TimeSteppingScheme::imex_11:
185 [[fallthrough]];
186 case TimeSteppingScheme::imex_22:
187 [[fallthrough]];
188 case TimeSteppingScheme::imex_33:
189 [[fallthrough]];
190 case TimeSteppingScheme::strang_ssprk_33_cn:
191 [[fallthrough]];
192 case TimeSteppingScheme::strang_erk_33_cn:
193 [[fallthrough]];
194 case TimeSteppingScheme::strang_erk_43_cn: {
195 AssertThrow(
196 !ParabolicSystem::is_identity,
197 dealii::ExcMessage(
198 "The selected equation has a trivial parabolic subsystem and "
199 "should not be run with an IMEX timestepping scheme."));
200 break;
201 }
202 }
203 };
204
205 check_whether_timestepping_makes_sense();
206 this->parse_parameters_call_back.connect(
207 check_whether_timestepping_makes_sense);
208 }
209
210
211 /*
212 * -------------------------------------------------------------------------
213 * Prepare state vector:
214 * -------------------------------------------------------------------------
215 */
216
217
218 template <typename HD, typename PD, int dim, typename Number>
219 void TimeIntegrator<HD, PD, dim, Number>::prepare_state_vector(
220 StateVector<Number> &state_vector, Number t) const
221 {
222 if (!ParabolicSystem::is_identity)
223 parabolic_module_->prepare_state_vector(state_vector, t);
224 hyperbolic_module_->prepare_state_vector(state_vector, t);
225 }
226
227
228 /*
229 * -------------------------------------------------------------------------
230 * High level step function implementing various CFLRecoveryStrategy
231 * -------------------------------------------------------------------------
232 */
233
234
235 template <typename HD, typename PD, int dim, typename Number>
236 Number TimeIntegrator<HD, PD, dim, Number>::step(
237 StateVector<Number> &state_vector,
238 Number t,
239 Number t_final /*=std::numeric_limits<Number>::max()*/)
240 {
241 Number tau_max =
242 std::min(tau_max_, t_final - t); /* enforces t <= t_final */
243
244#ifdef DEBUG_OUTPUT
245 std::cout << "TimeIntegrator<dim, Number>::step()" << std::endl;
246 std::cout << " enforcing tau_max <= " << tau_max << std::endl;
247#endif
248
249 const auto single_step = [&]() {
250 switch (time_stepping_scheme_) {
251 case TimeSteppingScheme::ssprk_22:
252 return step_ssprk_22(state_vector, t, tau_max);
253 case TimeSteppingScheme::ssprk_33:
254 return step_ssprk_33(state_vector, t, tau_max);
255 case TimeSteppingScheme::erk_11:
256 return step_erk_11(state_vector, t, tau_max);
257 case TimeSteppingScheme::erk_22:
258 return step_erk_22(state_vector, t, tau_max);
259 case TimeSteppingScheme::erk_33:
260 return step_erk_33(state_vector, t, tau_max);
261 case TimeSteppingScheme::erk_43:
262 return step_erk_43(state_vector, t, tau_max);
263 case TimeSteppingScheme::erk_54:
264 return step_erk_54(state_vector, t, tau_max);
265 case TimeSteppingScheme::strang_ssprk_33_cn:
266 return step_strang_ssprk_33_cn(state_vector, t, tau_max);
267 case TimeSteppingScheme::strang_erk_33_cn:
268 return step_strang_erk_33_cn(state_vector, t, tau_max);
269 case TimeSteppingScheme::strang_erk_43_cn:
270 return step_strang_erk_43_cn(state_vector, t, tau_max);
271 case TimeSteppingScheme::imex_11:
272 return step_imex_11(state_vector, t, tau_max);
273 case TimeSteppingScheme::imex_22:
274 return step_imex_22(state_vector, t, tau_max);
275 case TimeSteppingScheme::imex_33:
276 return step_imex_33(state_vector, t, tau_max);
277 default:
278 __builtin_unreachable();
279 }
280 };
281
282 if (cfl_recovery_strategy_ != CFLRecoveryStrategy::none) {
283 hyperbolic_module_->set_id_violation_strategy(
284 IDViolationStrategy::raise_exception);
285 parabolic_module_->set_id_violation_strategy(
286 IDViolationStrategy::raise_exception);
287 hyperbolic_module_->set_cfl(cfl_max_);
288 }
289
290 try {
291 return single_step();
292
293 } catch (const Restart &restart) {
294
295 AssertThrow(cfl_recovery_strategy_ != CFLRecoveryStrategy::none,
296 dealii::ExcInternalError());
297
298 hyperbolic_module_->set_id_violation_strategy(IDViolationStrategy::warn);
299 parabolic_module_->set_id_violation_strategy(IDViolationStrategy::warn);
300
301 if (cfl_recovery_strategy_ == CFLRecoveryStrategy::bang_bang_control) {
302 /* Retry with cfl_min instead of cfl_max: */
303#ifdef DEBUG_OUTPUT
304 std::cout
305 << " restart with bang bang control: setting cfl to cfl_min"
306 << std::endl;
307#endif
308 hyperbolic_module_->set_cfl(cfl_min_);
309 }
310
311 if (cfl_recovery_strategy_ == CFLRecoveryStrategy::cruise_control) {
312 /* Retry with the suggested tau_max: */
313#ifdef DEBUG_OUTPUT
314 std::cout
315 << " restart with cruise control: using suggested_tau_max"
316 << std::endl;
317#endif
318 //
319 // Multiply the suggested tau_max value with the efficiency.
320 //
321 // We have to account for the fact that the e Restart exception is
322 // thrown within a substep of the hyperbolic or parabolic module.
323 // This implies that the suggested_tau_max is computed for that
324 // particular substep and not for the full combined method (where
325 // tau_max can be larger). We thus multiply tau_max with the
326 // efficiency factor.
327 //
328 tau_max =
329 std::min(tau_max, efficiency_ * Number(restart.suggested_tau_max));
330 }
331
332 return single_step();
333 }
334 }
335
336
337 /*
338 * -------------------------------------------------------------------------
339 * Concrete implementation of ERK / IMEX time stepping strategies.
340 * -------------------------------------------------------------------------
341 */
342
343
344 template <typename HD, typename PD, int dim, typename Number>
345 Number TimeIntegrator<HD, PD, dim, Number>::step_ssprk_22(
346 StateVector<Number> &state_vector, Number t, Number tau_max)
347 {
348 /* SSP-RK2, see @cite Shu1988, Eq. 2.15. */
349
350#ifdef DEBUG_OUTPUT
351 std::cout << "TimeIntegrator<dim, Number>::step_ssprk_22()" << std::endl;
352#endif
353
354 Assert(efficiency_ == 1., dealii::ExcInternalError());
355
356 /* Step 1: T0 = U_old + tau * L(U_old) at t -> t + tau */
357 Number tau = hyperbolic_module_->template step<0>(
358 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
359
360 /* Step 2: T1 = T0 + tau L(T0) at time t + tau -> t + 2*tau */
361 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
362 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
363
364 /* Step 2: convex combination: T1 = 1/2 U_old + 1/2 T1 at time t + tau */
365 sadd<selected_memory_space_t>(
366 temp_[1], Number(1.0 / 2.0), Number(1.0 / 2.0), state_vector);
367
368 std::swap(state_vector, temp_[1]);
369 return tau;
370 }
371
372
373 template <typename HD, typename PD, int dim, typename Number>
374 Number TimeIntegrator<HD, PD, dim, Number>::step_ssprk_33(
375 StateVector<Number> &state_vector, Number t, Number tau_max)
376 {
377 /* SSP-RK3, see @cite Shu1988, Eq. 2.18. */
378
379#ifdef DEBUG_OUTPUT
380 std::cout << "TimeIntegrator<dim, Number>::step_ssprk_33()" << std::endl;
381#endif
382
383 Assert(efficiency_ == 1., dealii::ExcInternalError());
384
385 /* Step 1: T0 = U_old + tau * L(U_old) at time t -> t + tau */
386 Number tau = hyperbolic_module_->template step<0>(
387 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
388
389 /* Step 2: T1 = T0 + tau L(T0) at time t + tau -> t + 2*tau */
390 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
391 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
392
393 /* Step 2: convex combination T1 = 3/4 U_old + 1/4 T1 at time t + 0.5*tau */
394 sadd<selected_memory_space_t>(
395 temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), state_vector);
396
397 /* Step 3: T0 = T1 + tau L(T1) at time t + 0.5*tau -> t + 1.5*tau */
398 hyperbolic_module_->prepare_state_vector(temp_[1], t + 0.5 * tau);
399 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
400
401 /* Step 3: convex combination: T0 = 1/3 U_old + 2/3 T0 at time t + tau */
402 sadd<selected_memory_space_t>(
403 temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), state_vector);
404
405 std::swap(state_vector, temp_[0]);
406 return tau;
407 }
408
409
410 template <typename HD, typename PD, int dim, typename Number>
411 Number TimeIntegrator<HD, PD, dim, Number>::step_erk_11(
412 StateVector<Number> &state_vector, Number /*t*/, Number tau_max)
413 {
414#ifdef DEBUG_OUTPUT
415 std::cout << "TimeIntegrator<dim, Number>::step_erk_11()" << std::endl;
416#endif
417
418 Assert(efficiency_ == 1., dealii::ExcInternalError());
419
420 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
421 Number tau = hyperbolic_module_->template step<0>(
422 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
423
424 std::swap(state_vector, temp_[0]);
425 return tau;
426 }
427
428
429 template <typename HD, typename PD, int dim, typename Number>
430 Number TimeIntegrator<HD, PD, dim, Number>::step_erk_22(
431 StateVector<Number> &state_vector, Number t, Number tau_max)
432 {
433#ifdef DEBUG_OUTPUT
434 std::cout << "TimeIntegrator<dim, Number>::step_erk_22()" << std::endl;
435#endif
436
437 Assert(efficiency_ == 2., dealii::ExcInternalError());
438
439 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
440 Number tau = hyperbolic_module_->template step<0>(
441 state_vector, {}, {}, temp_[0], Number(.0), tau_max / efficiency_);
442
443 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + tau -> t + 2*tau */
444 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
445 hyperbolic_module_->template step<1>(
446 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
447
448 std::swap(state_vector, temp_[1]);
449 return efficiency_ * tau;
450 }
451
452
453 template <typename HD, typename PD, int dim, typename Number>
454 Number TimeIntegrator<HD, PD, dim, Number>::step_erk_33(
455 StateVector<Number> &state_vector, Number t, Number tau_max)
456 {
457#ifdef DEBUG_OUTPUT
458 std::cout << "TimeIntegrator<dim, Number>::step_erk_33()" << std::endl;
459#endif
460
461 Assert(efficiency_ == 3., dealii::ExcInternalError());
462
463 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
464 Number tau = hyperbolic_module_->template step<0>(
465 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
466
467 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + 1*tau -> t + 2*tau */
468 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
469 hyperbolic_module_->template step<1>(
470 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
471
472 /*
473 * Step 3: T2 <- {T1, 9/4} and {T0, -2} and {U_old, 3/4}
474 * at time t + 2*tau -> t + 3*tau
475 */
476 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
477 hyperbolic_module_->template step<2>(temp_[1],
478 {{state_vector, temp_[0]}},
479 {{Number(0.75), Number(-2.)}},
480 temp_[2],
481 tau);
482
483 std::swap(state_vector, temp_[2]);
484 return efficiency_ * tau;
485 }
486
487
488 template <typename HD, typename PD, int dim, typename Number>
489 Number TimeIntegrator<HD, PD, dim, Number>::step_erk_43(
490 StateVector<Number> &state_vector, Number t, Number tau_max)
491 {
492#ifdef DEBUG_OUTPUT
493 std::cout << "TimeIntegrator<dim, Number>::step_erk_43()" << std::endl;
494#endif
495
496 Assert(efficiency_ == 4., dealii::ExcInternalError());
497
498 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
499 Number tau = hyperbolic_module_->template step<0>(
500 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
501
502 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + 1*tau -> t + 2*tau */
503 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
504 hyperbolic_module_->template step<1>(
505 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
506
507 /* Step 3: T2 <- {T1, 2} and {T0, -1} at time t + 2*tau -> t + 3*tau */
508 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
509 hyperbolic_module_->template step<1>(
510 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
511
512 /*
513 * Step 4: T3 <- {T2, 8/3} and {T1,-10/3} and {T0, 5/3}
514 * at time t + 3*tau -> t + 4*tau
515 */
516 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
517 hyperbolic_module_->template step<2>(temp_[2],
518 {{temp_[0], temp_[1]}},
519 {{Number(5. / 3.), Number(-10. / 3.)}},
520 temp_[3],
521 tau);
522
523 std::swap(state_vector, temp_[3]);
524 return efficiency_ * tau;
525 }
526
527
528 template <typename HD, typename PD, int dim, typename Number>
529 Number TimeIntegrator<HD, PD, dim, Number>::step_erk_54(
530 StateVector<Number> &state_vector, Number t, Number tau_max)
531 {
532#ifdef DEBUG_OUTPUT
533 std::cout << "TimeIntegrator<dim, Number>::step_erk_54()" << std::endl;
534#endif
535
536 Assert(efficiency_ == 5., dealii::ExcInternalError());
537
538 constexpr Number c = 0.2; /* equidistant c_i */
539 constexpr Number a_21 = +0.2;
540 constexpr Number a_31 = +0.26075582269554909;
541 constexpr Number a_32 = +0.13924417730445096;
542 constexpr Number a_41 = -0.25856517872570289;
543 constexpr Number a_42 = +0.91136274166280729;
544 constexpr Number a_43 = -0.05279756293710430;
545 constexpr Number a_51 = +0.21623276431503774;
546 constexpr Number a_52 = +0.51534223099602405;
547 constexpr Number a_53 = -0.81662794199265554;
548 constexpr Number a_54 = +0.88505294668159373;
549 constexpr Number a_61 = -0.10511678454691901; /* aka b_1 */
550 constexpr Number a_62 = +0.87880047152100838; /* aka b_2 */
551 constexpr Number a_63 = -0.58903404061484477; /* aka b_3 */
552 constexpr Number a_64 = +0.46213380485434047; /* aka b_4 */
553 constexpr Number a_65 [[maybe_unused]] = +0.35321654878641495; /* aka b_5 */
554
555 /* Step 1: at time t -> t + 1*tau */
556 Number tau = hyperbolic_module_->template step<0>(
557 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
558
559 /* Step 2: at time t + 1*tau -> t + 2*tau */
560 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
561 hyperbolic_module_->template step<1>(
562 temp_[0], {{state_vector}}, {{(a_31 - a_21) / c}}, temp_[1], tau);
563
564 /* Step 3: at time t + 2*tau -> t + 3*tau */
565 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
566 hyperbolic_module_->template step<2>(
567 temp_[1],
568 {{state_vector, temp_[0]}},
569 {{(a_41 - a_31) / c, (a_42 - a_32) / c}},
570 temp_[2],
571 tau);
572
573 /* Step 4: at time t + 3*tau -> t + 4*tau */
574 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
575 hyperbolic_module_->template step<3>(
576 temp_[2],
577 {{state_vector, temp_[0], temp_[1]}},
578 {{(a_51 - a_41) / c, (a_52 - a_42) / c, (a_53 - a_43) / c}},
579 temp_[3],
580 tau);
581
582 /* Step 5: at time t + 4*tau -> t + 5*tau */
583 hyperbolic_module_->prepare_state_vector(temp_[3], t + 4.0 * tau);
584 hyperbolic_module_->template step<4>(
585 temp_[3],
586 {{state_vector, temp_[0], temp_[1], temp_[2]}},
587 {{(a_61 - a_51) / c,
588 (a_62 - a_52) / c,
589 (a_63 - a_53) / c,
590 (a_64 - a_54) / c}},
591 temp_[4],
592 tau);
593
594 std::swap(state_vector, temp_[4]);
595 return efficiency_ * tau;
596 }
597
598
599 template <typename HD, typename PD, int dim, typename Number>
600 Number TimeIntegrator<HD, PD, dim, Number>::step_strang_ssprk_33_cn(
601 StateVector<Number> &state_vector, Number t, Number tau_max)
602 {
603 // FIXME: avoid code duplication with step_ssprk_33
604
605#ifdef DEBUG_OUTPUT
606 std::cout << "TimeIntegrator<dim, Number>::step_strang_ssprk_33_cn()"
607 << std::endl;
608#endif
609
610 Assert(efficiency_ == 2., dealii::ExcInternalError());
611
612 /* First explicit SSPRK 3 step with final result in temp_[0]: */
613
614 Number tau = hyperbolic_module_->template step<0>(
615 state_vector, {}, {}, temp_[0], Number(0.0), tau_max / efficiency_);
616
617 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
618 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
619 sadd<selected_memory_space_t>(
620 temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), state_vector);
621
622 hyperbolic_module_->prepare_state_vector(temp_[1], t + 0.5 * tau);
623 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
624 sadd<selected_memory_space_t>(
625 temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), state_vector);
626
627 /* Implicit Crank-Nicolson step with final result in temp_[2]: */
628
629 try {
630 parabolic_module_->crank_nicolson_step(temp_[0], t, temp_[2], 2.0 * tau);
631 } catch (Restart &restart) {
632 /* Adjust suggested_tau_max. We multiply with efficiency_ again later */
633 restart.suggested_tau_max /= efficiency_;
634 throw;
635 }
636
637 /* Second SSPRK 3 step with final result in temp_[0]: */
638
639 hyperbolic_module_->prepare_state_vector( temp_[2], t + 1.0 * tau);
640 hyperbolic_module_->template step<0>( temp_[2], {}, {}, temp_[0], tau);
641
642 hyperbolic_module_->prepare_state_vector(temp_[0], t + 2.0 * tau);
643 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
644 sadd<selected_memory_space_t>(
645 temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), temp_[2]);
646
647 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.5 * tau);
648 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
649 sadd<selected_memory_space_t>(
650 temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), temp_[2]);
651
652 std::swap(state_vector, temp_[0]);
653 return efficiency_ * tau;
654 }
655
656
657 template <typename HD, typename PD, int dim, typename Number>
658 Number TimeIntegrator<HD, PD, dim, Number>::step_strang_erk_33_cn(
659 StateVector<Number> &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 std::swap(state_vector, temp_[2]);
714 return efficiency_ * tau;
715 }
716
717
718 template <typename HD, typename PD, int dim, typename Number>
719 Number TimeIntegrator<HD, PD, dim, Number>::step_strang_erk_43_cn(
720 StateVector<Number> &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 std::swap(state_vector, temp_[3]);
783 return efficiency_ * tau;
784 }
785
786
787 template <typename HD, typename PD, int dim, typename Number>
788 Number TimeIntegrator<HD, PD, dim, Number>::step_imex_11(
789 StateVector<Number> &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 std::swap(state_vector, temp_[1]);
806 return tau;
807 }
808
809
810 template <typename HD, typename PD, int dim, typename Number>
811 Number TimeIntegrator<HD, PD, dim, Number>::step_imex_22(
812 StateVector<Number> &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 std::swap(state_vector, temp_[3]);
842 return efficiency_ * tau;
843 }
844
845
846 template <typename HD, typename PD, int dim, typename Number>
847 Number TimeIntegrator<HD, PD, dim, Number>::step_imex_33(
848 StateVector<Number> &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 std::swap(state_vector, temp_[5]);
910 return efficiency_ * tau;
911 }
912
913} /* namespace ryujin */
ryujin::TimeIntegrator::prepare_state_vector
void prepare_state_vector(StateVector< Number > &state_vector, Number t) const
Definition time_integrator.template.h:219
ryujin::TimeIntegrator::TimeIntegrator
TimeIntegrator(const MPIEnsemble &mpi_ensemble, const OfflineData< dim, Number > &offline_data, const HyperbolicModule< HyperbolicDescription, dim, Number > &hyperbolic_module, const ParabolicModule< HyperbolicDescription, ParabolicDescription, dim, Number > &parabolic_module, const std::string &subsection="/TimeIntegrator")
Definition time_integrator.template.h:17
ryujin::TimeIntegrator::step
Number step(StateVector< Number > &state_vector, Number t, Number t_final=std::numeric_limits< Number >::max())
Definition time_integrator.template.h:236
ryujin::IDViolationStrategy::warn
@ warn
ryujin::IDViolationStrategy::raise_exception
@ raise_exception
ryujin::ParabolicModule
typename ParabolicDescription::template ParabolicModule< HyperbolicDescription, dim, Number > ParabolicModule
Definition parabolic_module.h:21
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
time_integrator.h
Generated by doxygen 1.9.8