ryujin 2.1.1 revision c38afcfdec63b34961924f05408760247496a1f0
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
18 template <typename StateVector, typename Number>
19 void
20 sadd(StateVector &dst, const Number s, const Number b, const StateVector &src)
21 {
22 auto &dst_U = std::get<0>(dst);
23 auto &src_U = std::get<0>(src);
24 dst_U.sadd(s, b, src_U);
25 }
26
27
28 template <typename Description, int dim, typename Number>
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::bang_bang_control;
56 add_parameter("cfl recovery strategy",
57 cfl_recovery_strategy_,
58 "CFL/invariant domain violation recovery strategy: none, "
59 "bang bang control");
60
62 time_stepping_scheme_ = TimeSteppingScheme::erk_33;
63 else
64 time_stepping_scheme_ = TimeSteppingScheme::strang_erk_33_cn;
65 add_parameter("time stepping scheme",
66 time_stepping_scheme_,
67 "Time stepping scheme: ssprk 22, ssprk 33, erk 11, erk 22, "
68 "erk 33, erk 43, erk "
69 "54, strang ssprk 33 cn, strang erk 33 cn, strang erk 43 cn, "
70 "imex 11, imex 22, imex 33");
71 }
72
73
74 template <typename Description, int dim, typename Number>
76 {
77#ifdef DEBUG_OUTPUT
78 std::cout << "TimeIntegrator<dim, Number>::prepare()" << std::endl;
79#endif
80
81 /* Resize temporary storage to appropriate sizes: */
82
83 switch (time_stepping_scheme_) {
85 temp_.resize(2);
86 efficiency_ = 1.;
87 break;
89 temp_.resize(2);
90 efficiency_ = 1.;
91 break;
93 temp_.resize(1);
94 efficiency_ = 1.;
95 break;
97 temp_.resize(2);
98 efficiency_ = 2.;
99 break;
101 temp_.resize(3);
102 efficiency_ = 3.;
103 break;
105 temp_.resize(4);
106 efficiency_ = 4.;
107 break;
109 temp_.resize(5);
110 efficiency_ = 5.;
111 break;
113 temp_.resize(3);
114 efficiency_ = 2.;
115 break;
117 temp_.resize(4);
118 efficiency_ = 6.;
119 break;
121 temp_.resize(4);
122 efficiency_ = 8.;
123 break;
125 temp_.resize(2);
126 efficiency_ = 1.;
127 break;
129 temp_.resize(4);
130 efficiency_ = 2.;
131 break;
133 temp_.resize(6);
134 efficiency_ = 3.;
135 break;
136 }
137
138 /* Initialize temporary vectors: */
139
140 for (auto &it : temp_) {
141 Vectors::reinit_state_vector<Description>(it, *offline_data_);
142 }
143
144 /* Reset CFL to canonical starting value: */
145
146 AssertThrow(cfl_min_ > 0., ExcMessage("cfl min must be a positive value"));
147 AssertThrow(cfl_max_ >= cfl_min_,
148 ExcMessage("cfl max must be greater than or equal to cfl min"));
149
150 hyperbolic_module_->cfl(cfl_max_);
151
152 const auto check_whether_timestepping_makes_sense = [&]() {
153 /*
154 * Make sure the user selects an appropriate time-stepping scheme.
155 */
156
157 switch (time_stepping_scheme_) {
159 [[fallthrough]];
161 [[fallthrough]];
163 [[fallthrough]];
165 [[fallthrough]];
167 [[fallthrough]];
169 [[fallthrough]];
171 AssertThrow(
173 dealii::ExcMessage(
174 "The selected equation consists of a hyperbolic and nontrivial "
175 "parabolic subsystem and requires an IMEX timestepping "
176 "scheme such as »strang erk 33 cn«."));
177 break;
178 }
180 [[fallthrough]];
182 [[fallthrough]];
184 [[fallthrough]];
186 [[fallthrough]];
188 [[fallthrough]];
190 AssertThrow(
192 dealii::ExcMessage(
193 "The selected equation has a trivial parabolic subsystem and "
194 "should not be run with an IMEX timestepping scheme."));
195 break;
196 }
197 }
198 };
199
200 check_whether_timestepping_makes_sense();
201 this->parse_parameters_call_back.connect(
202 check_whether_timestepping_makes_sense);
203 }
204
205
206 template <typename Description, int dim, typename Number>
208 StateVector &state_vector,
209 Number t,
210 Number t_final /*=std::numeric_limits<Number>::max()*/)
211 {
212#ifdef DEBUG_OUTPUT
213 std::cout << "TimeIntegrator<dim, Number>::step()" << std::endl;
214#endif
215 Number tau_max = t_final - t;
216
217 const auto single_step = [&]() {
218 switch (time_stepping_scheme_) {
220 return step_ssprk_22(state_vector, t, tau_max);
222 return step_ssprk_33(state_vector, t, tau_max);
224 return step_erk_11(state_vector, t, tau_max);
226 return step_erk_22(state_vector, t, tau_max);
228 return step_erk_33(state_vector, t, tau_max);
230 return step_erk_43(state_vector, t, tau_max);
232 return step_erk_54(state_vector, t, tau_max);
234 return step_strang_ssprk_33_cn(state_vector, t, tau_max);
236 return step_strang_erk_33_cn(state_vector, t, tau_max);
238 return step_strang_erk_43_cn(state_vector, t, tau_max);
240 return step_imex_11(state_vector, t, tau_max);
242 return step_imex_22(state_vector, t, tau_max);
244 return step_imex_33(state_vector, t, tau_max);
245 default:
246 __builtin_unreachable();
247 }
248 };
249
250 if (cfl_recovery_strategy_ == CFLRecoveryStrategy::bang_bang_control) {
251 hyperbolic_module_->id_violation_strategy_ =
253 parabolic_module_->id_violation_strategy_ =
255 hyperbolic_module_->cfl(cfl_max_);
256 }
257
258 try {
259 return single_step();
260
261 } catch (Restart) {
262
263 AssertThrow(cfl_recovery_strategy_ != CFLRecoveryStrategy::none,
264 dealii::ExcInternalError());
265
266 if (cfl_recovery_strategy_ == CFLRecoveryStrategy::bang_bang_control) {
267 hyperbolic_module_->id_violation_strategy_ = IDViolationStrategy::warn;
268 parabolic_module_->id_violation_strategy_ = IDViolationStrategy::warn;
269 hyperbolic_module_->cfl(cfl_min_);
270 return single_step();
271 }
272
273 __builtin_unreachable();
274 }
275 }
276
277
278 template <typename Description, int dim, typename Number>
280 StateVector &state_vector, Number t, Number tau_max)
281 {
282 /* SSP-RK2, see @cite Shu1988, Eq. 2.15. */
283
284 /* Step 1: T0 = U_old + tau * L(U_old) at t -> t + tau */
285 hyperbolic_module_->prepare_state_vector(state_vector, t);
286 Number tau = hyperbolic_module_->template step<0>(
287 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
288
289 /* Step 2: T1 = T0 + tau L(T0) at time t + tau -> t + 2*tau */
290 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
291 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
292
293 /* Step 2: convex combination: T1 = 1/2 U_old + 1/2 T1 at time t + tau */
294 sadd(temp_[1], Number(1.0 / 2.0), Number(1.0 / 2.0), state_vector);
295
296 state_vector.swap(temp_[1]);
297 return tau;
298 }
299
300
301 template <typename Description, int dim, typename Number>
303 StateVector &state_vector, Number t, Number tau_max)
304 {
305 /* SSP-RK3, see @cite Shu1988, Eq. 2.18. */
306
307 /* Step 1: T0 = U_old + tau * L(U_old) at time t -> t + tau */
308 hyperbolic_module_->prepare_state_vector(state_vector, t);
309 Number tau = hyperbolic_module_->template step<0>(
310 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
311
312 /* Step 2: T1 = T0 + tau L(T0) at time t + tau -> t + 2*tau */
313 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
314 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
315
316 /* Step 2: convex combination T1 = 3/4 U_old + 1/4 T1 at time t + 0.5*tau */
317 sadd(temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), state_vector);
318
319 /* Step 3: T0 = T1 + tau L(T1) at time t + 0.5*tau -> t + 1.5*tau */
320 hyperbolic_module_->prepare_state_vector(temp_[1], t + 0.5 * tau);
321 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
322
323 /* Step 3: convex combination: T0 = 1/3 U_old + 2/3 T0 at time t + tau */
324 sadd(temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), state_vector);
325
326 state_vector.swap(temp_[0]);
327 return tau;
328 }
329
330
331 template <typename Description, int dim, typename Number>
333 StateVector &state_vector, Number t, Number tau_max)
334 {
335#ifdef DEBUG_OUTPUT
336 std::cout << "TimeIntegrator<dim, Number>::step_erk_11()" << std::endl;
337#endif
338
339 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
340 hyperbolic_module_->prepare_state_vector(state_vector, t);
341 Number tau = hyperbolic_module_->template step<0>(
342 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
344 state_vector.swap(temp_[0]);
345 return tau;
346 }
347
348
349 template <typename Description, int dim, typename Number>
351 StateVector &state_vector, Number t, Number tau_max)
352 {
353#ifdef DEBUG_OUTPUT
354 std::cout << "TimeIntegrator<dim, Number>::step_erk_22()" << std::endl;
355#endif
356
357 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
358 hyperbolic_module_->prepare_state_vector(state_vector, t);
359 Number tau = hyperbolic_module_->template step<0>(
360 state_vector, {}, {}, temp_[0], Number(.0), tau_max / 2.);
361
362 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + tau -> t + 2*tau */
363 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
364 hyperbolic_module_->template step<1>(
365 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
366
367 state_vector.swap(temp_[1]);
368 return 2. * tau;
369 }
370
371
372 template <typename Description, int dim, typename Number>
374 StateVector &state_vector, Number t, Number tau_max)
375 {
376#ifdef DEBUG_OUTPUT
377 std::cout << "TimeIntegrator<dim, Number>::step_erk_33()" << std::endl;
378#endif
379
380 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
381 hyperbolic_module_->prepare_state_vector(state_vector, t);
382 Number tau = hyperbolic_module_->template step<0>(
383 state_vector, {}, {}, temp_[0], Number(0.), tau_max / 3.);
384
385 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + 1*tau -> t + 2*tau */
386 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
387 hyperbolic_module_->template step<1>(
388 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
389
390 /*
391 * Step 3: T2 <- {T1, 9/4} and {T0, -2} and {U_old, 3/4}
392 * at time t + 2*tau -> t + 3*tau
393 */
394 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
395 hyperbolic_module_->template step<2>(temp_[1],
396 {{state_vector, temp_[0]}},
397 {{Number(0.75), Number(-2.)}},
398 temp_[2],
399 tau);
401 state_vector.swap(temp_[2]);
402 return 3. * tau;
403 }
404
405
406 template <typename Description, int dim, typename Number>
408 StateVector &state_vector, Number t, Number tau_max)
409 {
410#ifdef DEBUG_OUTPUT
411 std::cout << "TimeIntegrator<dim, Number>::step_erk_43()" << std::endl;
412#endif
413
414 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
415 hyperbolic_module_->prepare_state_vector(state_vector, t);
416 Number tau = hyperbolic_module_->template step<0>(
417 state_vector, {}, {}, temp_[0], Number(0.), tau_max / 4.);
418
419 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + 1*tau -> t + 2*tau */
420 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
421 hyperbolic_module_->template step<1>(
422 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
423
424 /* Step 3: T2 <- {T1, 2} and {T0, -1} at time t + 2*tau -> t + 3*tau */
425 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
426 hyperbolic_module_->template step<1>(
427 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
428
429 /*
430 * Step 4: T3 <- {T2, 8/3} and {T1,-10/3} and {T0, 5/3}
431 * at time t + 3*tau -> t + 4*tau
432 */
433 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
434 hyperbolic_module_->template step<2>(temp_[2],
435 {{temp_[0], temp_[1]}},
436 {{Number(5. / 3.), Number(-10. / 3.)}},
437 temp_[3],
438 tau);
439
440 state_vector.swap(temp_[3]);
441 return 4. * tau;
442 }
443
444
445 template <typename Description, int dim, typename Number>
447 StateVector &state_vector, Number t, Number tau_max)
448 {
449#ifdef DEBUG_OUTPUT
450 std::cout << "TimeIntegrator<dim, Number>::step_erk_54()" << std::endl;
451#endif
452
453 constexpr Number c = 0.2; /* equidistant c_i */
454 constexpr Number a_21 = +0.2;
455 constexpr Number a_31 = +0.26075582269554909;
456 constexpr Number a_32 = +0.13924417730445096;
457 constexpr Number a_41 = -0.25856517872570289;
458 constexpr Number a_42 = +0.91136274166280729;
459 constexpr Number a_43 = -0.05279756293710430;
460 constexpr Number a_51 = +0.21623276431503774;
461 constexpr Number a_52 = +0.51534223099602405;
462 constexpr Number a_53 = -0.81662794199265554;
463 constexpr Number a_54 = +0.88505294668159373;
464 constexpr Number a_61 = -0.10511678454691901; /* aka b_1 */
465 constexpr Number a_62 = +0.87880047152100838; /* aka b_2 */
466 constexpr Number a_63 = -0.58903404061484477; /* aka b_3 */
467 constexpr Number a_64 = +0.46213380485434047; /* aka b_4 */
468 // constexpr Number a_65 = +0.35321654878641495; /* aka b_5 */
469
470 /* Step 1: at time t -> t + 1*tau */
471 hyperbolic_module_->prepare_state_vector(state_vector, t);
472 Number tau = hyperbolic_module_->template step<0>(
473 state_vector, {}, {}, temp_[0], Number(0.), tau_max / 5.);
474
475 /* Step 2: 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}}, {{(a_31 - a_21) / c}}, temp_[1], tau);
479
480 /* Step 3: at time t + 2*tau -> t + 3*tau */
481 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
482 hyperbolic_module_->template step<2>(
483 temp_[1],
484 {{state_vector, temp_[0]}},
485 {{(a_41 - a_31) / c, (a_42 - a_32) / c}},
486 temp_[2],
487 tau);
488
489 /* Step 4: at time t + 3*tau -> t + 4*tau */
490 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
491 hyperbolic_module_->template step<3>(
492 temp_[2],
493 {{state_vector, temp_[0], temp_[1]}},
494 {{(a_51 - a_41) / c, (a_52 - a_42) / c, (a_53 - a_43) / c}},
495 temp_[3],
496 tau);
497
498 /* Step 5: at time t + 4*tau -> t + 5*tau */
499 hyperbolic_module_->prepare_state_vector(temp_[3], t + 4.0 * tau);
500 hyperbolic_module_->template step<4>(
501 temp_[3],
502 {{state_vector, temp_[0], temp_[1], temp_[2]}},
503 {{(a_61 - a_51) / c,
504 (a_62 - a_52) / c,
505 (a_63 - a_53) / c,
506 (a_64 - a_54) / c}},
507 temp_[4],
508 tau);
509
510 state_vector.swap(temp_[4]);
511 return 5. * tau;
512 }
513
514
515 template <typename Description, int dim, typename Number>
517 StateVector &state_vector, Number t, Number tau_max)
518 {
519 // FIXME: avoid code duplication with step_ssprk_33
520
521#ifdef DEBUG_OUTPUT
522 std::cout << "TimeIntegrator<dim, Number>::step_strang_ssprk_33_cn()"
523 << std::endl;
524#endif
525
526 /* First explicit SSPRK 3 step with final result in temp_[0]: */
527
528 hyperbolic_module_->prepare_state_vector( state_vector, t);
529 Number tau = hyperbolic_module_->template step<0>(
530 state_vector, {}, {}, temp_[0], Number(0.0), tau_max / 2.);
531
532 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
533 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
534 sadd(temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), state_vector);
535
536 hyperbolic_module_->prepare_state_vector(temp_[1], t + 0.5 * tau);
537 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
538 sadd(temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), state_vector);
539
540 /* Implicit Crank-Nicolson step with final result in temp_[2]: */
541
542 parabolic_module_->template step<0>(temp_[0], t, {}, {}, temp_[2], tau);
543 sadd(temp_[2], Number(2.), Number(-1.), temp_[0]);
544
545 /* Second SSPRK 3 step with final result in temp_[0]: */
546
547 hyperbolic_module_->prepare_state_vector( temp_[2], t + 1.0 * tau);
548 hyperbolic_module_->template step<0>( temp_[2], {}, {}, temp_[0], tau);
549
550 hyperbolic_module_->prepare_state_vector(temp_[0], t + 2.0 * tau);
551 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
552 sadd(temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), temp_[2]);
553
554 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.5 * tau);
555 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
556 sadd(temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), temp_[2]);
557
558 state_vector.swap(temp_[0]);
559 return 2. * tau;
560 }
561
562
563 template <typename Description, int dim, typename Number>
565 StateVector &state_vector, Number t, Number tau_max)
566 {
567 // FIXME: refactor to avoid code duplication with step_erk_33
568
569#ifdef DEBUG_OUTPUT
570 std::cout << "TimeIntegrator<dim, Number>::step_strang_erk_33_cn()"
571 << std::endl;
572#endif
573
574 /* First explicit ERK(3,3,1) step with final result in temp_[2]: */
575
576 hyperbolic_module_->prepare_state_vector(state_vector, t);
577 Number tau = hyperbolic_module_->template step<0>(
578 state_vector, {}, {}, temp_[0], Number(0.), tau_max / 6.);
579
580 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
581 hyperbolic_module_->template step<1>(
582 temp_[0], {{ state_vector}}, {{Number(-1.)}}, temp_[1], tau);
583
584 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
585 hyperbolic_module_->template step<2>(temp_[1],
586 {{ state_vector, temp_[0]}},
587 {{Number(0.75), Number(-2.)}},
588 temp_[2],
589 tau);
590
591 /* Implicit Crank-Nicolson step with final result in temp_[3]: */
592
593 parabolic_module_->template step<0>(
594 temp_[2], t, {}, {}, temp_[3], 3.0 * tau);
595 sadd(temp_[3], Number(2.), Number(-1.), temp_[2]);
596
597 /* Second explicit ERK(3,3,1) 3 step with final result in temp_[2]: */
598
599 hyperbolic_module_->prepare_state_vector(temp_[3], t + 3.0 * tau);
600 hyperbolic_module_->template step<0>(
601 temp_[3], {}, {}, temp_[0], tau);
602
603 hyperbolic_module_->prepare_state_vector(temp_[0], t + 4.0 * tau);
604 hyperbolic_module_->template step<1>(
605 temp_[0], {{ temp_[3]}}, {{Number(-1.)}}, temp_[1], tau);
606
607 hyperbolic_module_->prepare_state_vector(temp_[1], t + 5.0 * tau);
608 hyperbolic_module_->template step<2>(temp_[1],
609 {{ temp_[3], temp_[0]}},
610 {{Number(0.75), Number(-2.)}},
611 temp_[2],
612 tau);
613
614 state_vector.swap(temp_[2]);
615 return 6. * tau;
616 }
617
618
619 template <typename Description, int dim, typename Number>
621 StateVector &state_vector, Number t, Number tau_max)
622 {
623 // FIXME: refactor to avoid code duplication with step_erk_43
624
625#ifdef DEBUG_OUTPUT
626 std::cout << "TimeIntegrator<dim, Number>::step_strang_erk_43_cn()"
627 << std::endl;
628#endif
629
630 /* First explicit ERK(4,3,1) step with final result in temp_[3]: */
631
632 hyperbolic_module_->prepare_state_vector(state_vector, t);
633 Number tau = hyperbolic_module_->template step<0>(
634 state_vector, {}, {}, temp_[0], Number(0.), tau_max / 8.);
635
636 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
637 hyperbolic_module_->template step<1>(
638 temp_[0], {{ state_vector}}, {{Number(-1.)}}, temp_[1], tau);
639
640 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
641 hyperbolic_module_->template step<1>(
642 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
643
644 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
645 hyperbolic_module_->template step<2>(temp_[2],
646 {{temp_[0], temp_[1]}},
647 {{Number(5. / 3.), Number(-10. / 3.)}},
648 temp_[3],
649 tau);
650
651 /* Implicit Crank-Nicolson step with final result in temp_[2]: */
652
653 parabolic_module_->template step<0>(
654 temp_[3], t, {}, {}, temp_[2], 4.0 * tau);
655 sadd(temp_[2], Number(2.), Number(-1.), temp_[3]);
656
657 /* Second explicit ERK(4,3,1) step with final result in temp_[3]: */
658
659 hyperbolic_module_->prepare_state_vector(temp_[2], t + 4.0 * tau);
660 hyperbolic_module_->template step<0>(
661 temp_[2], {}, {}, temp_[0], tau);
662
663 hyperbolic_module_->prepare_state_vector(temp_[0], t + 5.0 * tau);
664 hyperbolic_module_->template step<1>(
665 temp_[0], {{ temp_[2]}}, {{Number(-1.)}}, temp_[1], tau);
666
667 hyperbolic_module_->prepare_state_vector(temp_[1], t + 6.0 * tau);
668 hyperbolic_module_->template step<1>(
669 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
670
671 hyperbolic_module_->prepare_state_vector(temp_[2], t + 7.0 * tau);
672 hyperbolic_module_->template step<2>(temp_[2],
673 {{temp_[0], temp_[1]}},
674 {{Number(5. / 3.), Number(-10. / 3.)}},
675 temp_[3],
676 tau);
677
678 state_vector.swap(temp_[3]);
679 return 8. * tau;
680 }
681
682 template <typename Description, int dim, typename Number>
684 StateVector &state_vector, Number t, Number tau_max)
685 {
686#ifdef DEBUG_OUTPUT
687 std::cout << "TimeIntegrator<dim, Number>::step_imex_11()" << std::endl;
688#endif
689
690 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
691 hyperbolic_module_->prepare_state_vector(state_vector, t);
692 Number tau = hyperbolic_module_->template step<0>(
693 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
694
695 /* Implicit step 1: T1 <- {T0, 1} at time t -> t + tau */
696 parabolic_module_->template step<0>(
697 temp_[0], t, {}, {}, temp_[1], 1.0 * tau);
698
699 state_vector.swap(temp_[1]);
700 return tau;
701 }
702
703 template <typename Description, int dim, typename Number>
705 StateVector &state_vector, Number t, Number tau_max)
706 {
707#ifdef DEBUG_OUTPUT
708 std::cout << "TimeIntegrator<dim, Number>::step_imex_22()" << std::endl;
709#endif
710
711 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
712 hyperbolic_module_->prepare_state_vector(state_vector, t);
713 Number tau = hyperbolic_module_->template step<0>(
714 state_vector, {}, {}, temp_[0], Number(0.), tau_max / 2.);
715
716 /* Implicit step 1: T1 <- {T0, 1} at time t -> t + tau */
717 parabolic_module_->template step<0>(temp_[0], t, {}, {}, temp_[1], tau);
718
719 /* Explicit step 2: T2 <- {T1, 2} and {U_old, -1} at t + tau -> t + 2 tau */
720 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.0 * tau);
721 hyperbolic_module_->template step<1>(
722 temp_[1], {{state_vector}}, {{Number(-1.)}}, temp_[2], tau);
723
724 /* Implicit step 2: T3 <- {T2, 0} and {U_old, 1} at t + tau -> t + 2 tau */
725 parabolic_module_->template step<1>(temp_[2],
726 t + 1.0 * tau,
727 {{state_vector}},
728 {{Number(1.)}},
729 temp_[3],
730 tau);
731
732 state_vector.swap(temp_[3]);
733 return 2. * tau;
734 }
735
736 template <typename Description, int dim, typename Number>
738 StateVector &state_vector, Number t, Number tau_max)
739 {
740#ifdef DEBUG_OUTPUT
741 std::cout << "TimeIntegrator<dim, Number>::step_imex_33()" << std::endl;
742#endif
743
744 /* IMEX(3, 3; 1), see @cite ErnGuermond2023, Sec. 4.3. */
745
746 const Number gamma = 0.5 + 0.5 * (1. / std::sqrt(3.));
747
748 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
749 hyperbolic_module_->prepare_state_vector(state_vector, t);
750 Number tau = hyperbolic_module_->template step<0>(
751 state_vector, {}, {}, temp_[0], Number(0.), tau_max / 3.);
752
753 /* Implicit step 1: T1 <- {U_old, 1 - 3*gamma} at time t -> t + tau */
754 parabolic_module_->template step<1>(temp_[0],
755 t,
756 {{state_vector}},
757 {{Number(1. - 3. * gamma)}},
758 temp_[1],
759 tau);
760
761 /* Explicit step 2: T2 <- {U_old, -1} and {T1, 2} at time t -> t + 2 tau */
762 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.0 * tau);
763 hyperbolic_module_->template step<1>(
764 temp_[1], {{state_vector}}, {{Number(-1.)}}, temp_[2], tau);
765
766 /*
767 * Implicit step 2:
768 * T3 <- {U_old, 6*gamma-1} and {T1, 2-9*gamma} at t -> t + * 2 tau
769 */
770 parabolic_module_->template step<2>(
771 temp_[2],
772 t + tau,
773 {{state_vector, temp_[1]}},
774 {{Number(6. * gamma - 1.), Number(2. - 9 * gamma)}},
775 temp_[3],
776 tau);
777
778 /* Explicit step 3: T4 <- {U_old, 3 / 4} and {T1, -2} at t -> t + 3 tau */
779 hyperbolic_module_->prepare_state_vector(temp_[3], t + 2. * tau);
780 hyperbolic_module_->template step<2>(temp_[3],
781 {{state_vector, temp_[1]}},
782 {{Number(0.75), Number(-2.)}},
783 temp_[4],
784 tau);
785
786 /* Implicit step 3: */
787 parabolic_module_->template step<3>(temp_[4],
788 t + 2. * tau,
789 {{state_vector, temp_[1], temp_[3]}},
790 {{Number(0.75 - 3. * gamma),
791 Number(6. * gamma - 2.),
792 Number(9. / 4. - 3. * gamma)}},
793 temp_[5],
794 tau);
795
796 state_vector.swap(temp_[5]);
797 return 3. * tau;
798 }
799
800} /* namespace ryujin */
Number step_imex_33(StateVector &state_vector, Number t, Number tau_max)
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")
Number step_strang_ssprk_33_cn(StateVector &state_vector, Number t, Number tau_max)
Number step_erk_22(StateVector &state_vector, Number t, Number tau_max)
Number step_imex_11(StateVector &state_vector, Number t, Number tau_max)
Number step_strang_erk_43_cn(StateVector &state_vector, Number t, Number tau_max)
Number step_erk_54(StateVector &state_vector, Number t, Number tau_max)
Number step_erk_33(StateVector &state_vector, Number t, Number tau_max)
Number step_erk_43(StateVector &state_vector, Number t, Number tau_max)
Number step_imex_22(StateVector &state_vector, Number t, Number tau_max)
typename View::StateVector StateVector
Number step_erk_11(StateVector &state_vector, Number t, Number tau_max)
Number step_ssprk_22(StateVector &state_vector, Number t, Number tau_max)
Number step_ssprk_33(StateVector &state_vector, Number t, Number tau_max)
Number step(StateVector &state_vector, Number t, Number t_final=std::numeric_limits< Number >::max())
Number step_strang_erk_33_cn(StateVector &state_vector, Number t, Number tau_max)
std::tuple< MultiComponentVector< Number, problem_dim >, MultiComponentVector< Number, prec_dim >, BlockVector< Number > > StateVector
Definition: state_vector.h:51
void sadd(StateVector &dst, const Number s, const Number b, const StateVector &src)