ryujin 2.1.1 revision 2062b43aae59b7ee510dec4ae666e29d56c4f0be
time_integrator.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: MIT
3// Copyright (C) 2020 - 2023 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 Description, int dim, typename Number>
16 const MPI_Comm &mpi_communicator,
17 std::map<std::string, dealii::Timer> &computing_timer,
18 const OfflineData<dim, Number> &offline_data,
19 const HyperbolicModule<Description, dim, Number> &hyperbolic_module,
20 const ParabolicModule<Description, dim, Number> &parabolic_module,
21 const std::string &subsection /*= "TimeIntegrator"*/)
22 : ParameterAcceptor(subsection)
23 , mpi_communicator_(mpi_communicator)
24 , computing_timer_(computing_timer)
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::bang_bang_control;
44 add_parameter("cfl recovery strategy",
45 cfl_recovery_strategy_,
46 "CFL/invariant domain violation recovery strategy: none, "
47 "bang bang control");
48
50 time_stepping_scheme_ = TimeSteppingScheme::erk_33;
51 else
52 time_stepping_scheme_ = TimeSteppingScheme::strang_erk_33_cn;
53 add_parameter("time stepping scheme",
54 time_stepping_scheme_,
55 "Time stepping scheme: ssprk 22, ssprk 33, erk 11, erk 22, "
56 "erk 33, erk 43, erk "
57 "54, strang ssprk 33 cn, strang erk 33 cn, strang erk 43 cn");
58 }
59
60
61 template <typename Description, int dim, typename Number>
63 {
64#ifdef DEBUG_OUTPUT
65 std::cout << "TimeIntegrator<dim, Number>::prepare()" << std::endl;
66#endif
67
68 /* Resize temporary storage to appropriate sizes: */
69
70 switch (time_stepping_scheme_) {
72 U_.resize(2);
73 precomputed_.resize(1);
74 efficiency_ = 1.;
75 break;
77 U_.resize(2);
78 precomputed_.resize(1);
79 efficiency_ = 1.;
80 break;
82 U_.resize(1);
83 precomputed_.resize(1);
84 efficiency_ = 1.;
85 break;
87 U_.resize(2);
88 precomputed_.resize(2);
89 efficiency_ = 2.;
90 break;
92 U_.resize(3);
93 precomputed_.resize(3);
94 efficiency_ = 3.;
95 break;
97 U_.resize(4);
98 precomputed_.resize(4);
99 efficiency_ = 4.;
100 break;
102 U_.resize(5);
103 precomputed_.resize(5);
104 efficiency_ = 5.;
105 break;
107 U_.resize(3);
108 precomputed_.resize(1);
109 efficiency_ = 2.;
110 break;
112 U_.resize(4);
113 precomputed_.resize(3);
114 efficiency_ = 6.;
115 break;
117 U_.resize(4); // FIXME
118 precomputed_.resize(4);
119 efficiency_ = 8.;
120 break;
121 }
122
123 /* Initialize temporary vectors and matrices: */
124
125 const auto &vector_partitioner = offline_data_->vector_partitioner();
126 for (auto &it : U_)
127 it.reinit(vector_partitioner);
128
129 const auto &scalar_partitioner = offline_data_->scalar_partitioner();
130 for (auto &it : precomputed_)
131 it.reinit_with_scalar_partitioner(scalar_partitioner);
132
133 /* Reset CFL to canonical starting value: */
134
135 AssertThrow(cfl_min_ > 0., ExcMessage("cfl min must be a positive value"));
136 AssertThrow(cfl_max_ >= cfl_min_,
137 ExcMessage("cfl max must be greater than or equal to cfl min"));
138
139 hyperbolic_module_->cfl(cfl_max_);
140
141 const auto check_whether_timestepping_makes_sense = [&]() {
142 /*
143 * Make sure the user selects an appropriate time-stepping scheme.
144 */
145
146 switch (time_stepping_scheme_) {
148 [[fallthrough]];
150 [[fallthrough]];
152 [[fallthrough]];
154 [[fallthrough]];
156 [[fallthrough]];
158 [[fallthrough]];
160 AssertThrow(
162 dealii::ExcMessage(
163 "The selected equation consists of a hyperbolic and nontrivial "
164 "parabolic subsystem and requires an IMEX timestepping "
165 "scheme such as »strang erk 33 cn«."));
166 break;
167 }
169 [[fallthrough]];
171 [[fallthrough]];
173 AssertThrow(
175 dealii::ExcMessage(
176 "The selected equation has a trivial parabolic subsystem and "
177 "should not be run with an IMEX timestepping scheme."));
178 break;
179 }
180 }
181 };
182
183 check_whether_timestepping_makes_sense();
184 this->parse_parameters_call_back.connect(
185 check_whether_timestepping_makes_sense);
186 }
187
188
189 template <typename Description, int dim, typename Number>
191 Number t)
192 {
193#ifdef DEBUG_OUTPUT
194 std::cout << "TimeIntegrator<dim, Number>::step()" << std::endl;
195#endif
196
197 const auto single_step = [&]() {
198 switch (time_stepping_scheme_) {
200 return step_ssprk_22(U, t);
202 return step_ssprk_33(U, t);
204 return step_erk_11(U, t);
206 return step_erk_22(U, t);
208 return step_erk_33(U, t);
210 return step_erk_43(U, t);
212 return step_erk_54(U, t);
214 return step_strang_ssprk_33_cn(U, t);
216 return step_strang_erk_33_cn(U, t);
218 return step_strang_erk_43_cn(U, t);
219 default:
220 __builtin_unreachable();
221 }
222 };
223
224 if (cfl_recovery_strategy_ == CFLRecoveryStrategy::bang_bang_control) {
225 hyperbolic_module_->id_violation_strategy_ =
227 parabolic_module_->id_violation_strategy_ =
229 hyperbolic_module_->cfl(cfl_max_);
230 }
231
232 try {
233 return single_step();
234
235 } catch (Restart) {
236
237 AssertThrow(cfl_recovery_strategy_ != CFLRecoveryStrategy::none,
238 dealii::ExcInternalError());
239
240 if (cfl_recovery_strategy_ == CFLRecoveryStrategy::bang_bang_control) {
241 hyperbolic_module_->id_violation_strategy_ = IDViolationStrategy::warn;
242 parabolic_module_->id_violation_strategy_ = IDViolationStrategy::warn;
243 hyperbolic_module_->cfl(cfl_min_);
244 return single_step();
245 }
246
247 __builtin_unreachable();
248 }
249 }
250
251 template <typename Description, int dim, typename Number>
253 Number t)
254 {
255 /* SSP-RK3, see @cite Shu1988, Eq. 2.15. */
256
257 /* Step 1: U1 = U_old + tau * L(U_old) at time t + tau */
258 Number tau = hyperbolic_module_->template step<0>(
259 U, {}, {}, {}, U_[0], precomputed_[0]);
260 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
261
262 /* Step 2: U2 = 1/2 U_old + 1/2 (U1 + tau L(U1)) at time t + tau */
263 hyperbolic_module_->template step<0>(
264 U_[0], {}, {}, {}, U_[1], precomputed_[0], tau);
265 U_[1].sadd(Number(1. / 2.), Number(1. / 2.), U);
266 hyperbolic_module_->apply_boundary_conditions(U_[1], t + tau);
267
268 U.swap(U_[1]);
269 return tau;
270 }
271
272
273 template <typename Description, int dim, typename Number>
275 Number t)
276 {
277 /* SSP-RK3, see @cite Shu1988, Eq. 2.18. */
278
279 /* Step 1: U1 = U_old + tau * L(U_old) at time t + tau */
280 Number tau = hyperbolic_module_->template step<0>(
281 U, {}, {}, {}, U_[0], precomputed_[0]);
282 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
283
284 /* Step 2: U2 = 3/4 U_old + 1/4 (U1 + tau L(U1)) at time t + 0.5 * tau */
285 hyperbolic_module_->template step<0>(
286 U_[0], {}, {}, {}, U_[1], precomputed_[0], tau);
287 U_[1].sadd(Number(1. / 4.), Number(3. / 4.), U);
288 hyperbolic_module_->apply_boundary_conditions(U_[1], t + 0.5 * tau);
289
290 /* Step 3: U3 = 1/3 U_old + 2/3 (U2 + tau L(U2)) at final time t + tau */
291 hyperbolic_module_->template step<0>(
292 U_[1], {}, {}, {}, U_[0], precomputed_[0], tau);
293 U_[0].sadd(Number(2. / 3.), Number(1. / 3.), U);
294 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
295
296 U.swap(U_[0]);
297 return tau;
298 }
300
301 template <typename Description, int dim, typename Number>
303 Number t)
304 {
305#ifdef DEBUG_OUTPUT
306 std::cout << "TimeIntegrator<dim, Number>::step_erk_11()" << std::endl;
307#endif
308
309 /* Step 1: U1 <- {U, 1} at time t + tau */
310 Number tau = hyperbolic_module_->template step<0>(
311 U, {}, {}, {}, U_[0], precomputed_[0]);
312 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
314 U.swap(U_[0]);
315 return tau;
316 }
317
318
319 template <typename Description, int dim, typename Number>
321 Number t)
322 {
323#ifdef DEBUG_OUTPUT
324 std::cout << "TimeIntegrator<dim, Number>::step_erk_22()" << std::endl;
325#endif
326
327 /* Step 1: U1 <- {U, 1} at time t + tau */
328 Number tau = hyperbolic_module_->template step<0>(
329 U, {}, {}, {}, U_[0], precomputed_[0]);
330 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
331
332 /* Step 2: U2 <- {U1, 2} and {U, -1} at time t + 2 tau */
333 hyperbolic_module_->template step<1>(U_[0],
334 {{U}},
335 {{precomputed_[0]}},
336 {{Number(-1.)}},
337 U_[1],
338 precomputed_[1],
339 tau);
340 hyperbolic_module_->apply_boundary_conditions(U_[1], t + 2. * tau);
341
342 U.swap(U_[1]);
343 return 2. * tau;
344 }
345
346
347 template <typename Description, int dim, typename Number>
349 Number t)
350 {
351#ifdef DEBUG_OUTPUT
352 std::cout << "TimeIntegrator<dim, Number>::step_erk_33()" << std::endl;
353#endif
354
355 /* Step 1: U1 <- {U, 1} at time t + tau */
356 Number tau = hyperbolic_module_->template step<0>(
357 U, {}, {}, {}, U_[0], precomputed_[0]);
358 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
359
360 /* Step 2: U2 <- {U1, 2} and {U, -1} at time t + 2 tau */
361 hyperbolic_module_->template step<1>(U_[0],
362 {{U}},
363 {{precomputed_[0]}},
364 {{Number(-1.)}},
365 U_[1],
366 precomputed_[1],
367 tau);
368 hyperbolic_module_->apply_boundary_conditions(U_[1], t + 2. * tau);
369
370 /* Step 3: U3 <- {U2, 9/4} and {U1, -2} and {U, 3/4} at time t + 3 tau */
371 hyperbolic_module_->template step<2>(U_[1],
372 {{U, U_[0]}},
373 {{precomputed_[0], precomputed_[1]}},
374 {{Number(0.75), Number(-2.)}},
375 U_[2],
376 precomputed_[2],
377 tau);
378 hyperbolic_module_->apply_boundary_conditions(U_[2], t + 3. * tau);
379
380 U.swap(U_[2]);
381 return 3. * tau;
382 }
383
384
385 template <typename Description, int dim, typename Number>
387 Number t)
388 {
389#ifdef DEBUG_OUTPUT
390 std::cout << "TimeIntegrator<dim, Number>::step_erk_43()" << std::endl;
391#endif
392
393 /* Step 1: U1 <- {U, 1} at time t + tau */
394 Number tau = hyperbolic_module_->template step<0>(
395 U, {}, {}, {}, U_[0], precomputed_[0]);
396 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
397
398 /* Step 2: U2 <- {U1, 2} and {U, -1} at time t + 2 tau */
399 hyperbolic_module_->template step<1>(U_[0],
400 {{U}},
401 {{precomputed_[0]}},
402 {{Number(-1.)}},
403 U_[1],
404 precomputed_[1],
405 tau);
406 hyperbolic_module_->apply_boundary_conditions(U_[1], t + 2. * tau);
407
408 /* Step 3: U3 <- {U2, 2} and {U1, -1} at time t + 3 tau */
409 hyperbolic_module_->template step<1>(U_[1],
410 {{U_[0]}},
411 {{precomputed_[1]}},
412 {{Number(-1.)}},
413 U_[2],
414 precomputed_[2],
415 tau);
416 hyperbolic_module_->apply_boundary_conditions(U_[2], t + 3. * tau);
417
418 /* Step 4: U4 <- {U3, 8/3} and {U2,-10/3} and {U1, 5/3} at time t + 4 tau */
419 hyperbolic_module_->template step<2>(U_[2],
420 {{U_[0], U_[1]}},
421 {{precomputed_[1], precomputed_[2]}},
422 {{Number(5. / 3.), Number(-10. / 3.)}},
423 U_[3],
424 precomputed_[3],
425 tau);
426 hyperbolic_module_->apply_boundary_conditions(U_[3], t + 4. * tau);
427
428 U.swap(U_[3]);
429 return 4. * tau;
430 }
431
432
433 template <typename Description, int dim, typename Number>
435 Number t)
436 {
437#ifdef DEBUG_OUTPUT
438 std::cout << "TimeIntegrator<dim, Number>::step_erk_54()" << std::endl;
439#endif
440
441 constexpr Number c = 0.2; /* equidistant c_i */
442 constexpr Number a_21 = +0.2;
443 constexpr Number a_31 = +0.26075582269554909;
444 constexpr Number a_32 = +0.13924417730445096;
445 constexpr Number a_41 = -0.25856517872570289;
446 constexpr Number a_42 = +0.91136274166280729;
447 constexpr Number a_43 = -0.05279756293710430;
448 constexpr Number a_51 = +0.21623276431503774;
449 constexpr Number a_52 = +0.51534223099602405;
450 constexpr Number a_53 = -0.81662794199265554;
451 constexpr Number a_54 = +0.88505294668159373;
452 constexpr Number a_61 = -0.10511678454691901; /* aka b_1 */
453 constexpr Number a_62 = +0.87880047152100838; /* aka b_2 */
454 constexpr Number a_63 = -0.58903404061484477; /* aka b_3 */
455 constexpr Number a_64 = +0.46213380485434047; /* aka b_4 */
456 // constexpr Number a_65 = +0.35321654878641495; /* aka b_5 */
457
458 /* Step 1: */
459 Number tau = hyperbolic_module_->template step<0>(
460 U, {}, {}, {}, U_[0], precomputed_[0]);
461 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
462
463 /* Step 2: */
464 hyperbolic_module_->template step<1>(U_[0],
465 {{U}},
466 {{precomputed_[0]}},
467 {{(a_31 - a_21) / c}},
468 U_[1],
469 precomputed_[1],
470 tau);
471 hyperbolic_module_->apply_boundary_conditions(U_[1], t + 2. * tau);
472
473 /* Step 3: */
474 hyperbolic_module_->template step<2>(
475 U_[1],
476 {{U, U_[0]}},
477 {{precomputed_[0], precomputed_[1]}},
478 {{(a_41 - a_31) / c, (a_42 - a_32) / c}},
479 U_[2],
480 precomputed_[2],
481 tau);
482 hyperbolic_module_->apply_boundary_conditions(U_[2], t + 3. * tau);
483
484 /* Step 4: */
485 hyperbolic_module_->template step<3>(
486 U_[2],
487 {{U, U_[0], U_[1]}},
488 {{precomputed_[0], precomputed_[1], precomputed_[2]}},
489 {{(a_51 - a_41) / c, (a_52 - a_42) / c, (a_53 - a_43) / c}},
490 U_[3],
491 precomputed_[3],
492 tau);
493 hyperbolic_module_->apply_boundary_conditions(U_[3], t + 4. * tau);
494
495 /* Step 5: */
496 hyperbolic_module_->template step<4>(
497 U_[3],
498 {{U, U_[0], U_[1], U_[2]}},
499 {{precomputed_[0], precomputed_[1], precomputed_[2], precomputed_[3]}},
500 {{(a_61 - a_51) / c,
501 (a_62 - a_52) / c,
502 (a_63 - a_53) / c,
503 (a_64 - a_54) / c}},
504 U_[4],
505 precomputed_[4],
506 tau);
507 hyperbolic_module_->apply_boundary_conditions(U_[4], t + 5. * tau);
508
509 U.swap(U_[4]);
510 return 5. * tau;
511 }
512
513
514 template <typename Description, int dim, typename Number>
516 vector_type &U, Number t)
517 {
518 // FIXME: avoid code duplication with step_ssprk_33
519
520#ifdef DEBUG_OUTPUT
521 std::cout << "TimeIntegrator<dim, Number>::step_strang_ssprk_33_cn()"
522 << std::endl;
523#endif
524
525 /* First explicit SSPRK 3 step with final result in U_[0]: */
526
527 Number tau = hyperbolic_module_->template step<0>(
528 /*input*/ U, {}, {}, {}, U_[0], precomputed_[0]);
529 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
530
531 hyperbolic_module_->template step<0>(
532 U_[0], {}, {}, {}, U_[1], precomputed_[0], tau);
533 U_[1].sadd(Number(1. / 4.), Number(3. / 4.), /*input*/ U);
534 hyperbolic_module_->apply_boundary_conditions(U_[1], t + 0.5 * tau);
535
536 hyperbolic_module_->template step<0>(
537 U_[1], {}, {}, {}, U_[0], precomputed_[0], tau);
538 U_[0].sadd(Number(2. / 3.), Number(1. / 3.), /*input*/ U);
539 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
540
541 /* Implicit Crank-Nicolson step with final result in U_[2]: */
542
543 parabolic_module_->crank_nicolson_step(U_[0], t, U_[2], 2.0 * tau);
544
545 /* Second SSPRK 3 step with final result in U_[0]: */
546
547 hyperbolic_module_->template step<0>(
548 /*intermediate*/ U_[2], {}, {}, {}, U_[0], precomputed_[0], tau);
549 hyperbolic_module_->apply_boundary_conditions(U_[0], t + 2.0 * tau);
550
551 hyperbolic_module_->template step<0>(
552 U_[0], {}, {}, {}, U_[1], precomputed_[0], tau);
553 U_[1].sadd(Number(1. / 4.), Number(3. / 4.), /*intermediate*/ U_[2]);
554 hyperbolic_module_->apply_boundary_conditions(U_[1], t + 1.5 * tau);
555
556 hyperbolic_module_->template step<0>(
557 U_[1], {}, {}, {}, U_[0], precomputed_[0], tau);
558 U_[0].sadd(Number(2. / 3.), Number(1. / 3.), /*intermediate*/ U_[2]);
559 hyperbolic_module_->apply_boundary_conditions(U_[0], t + 2.0 * tau);
560
561 U.swap(U_[0]);
562 return 2.0 * tau;
563 }
564
565
566 template <typename Description, int dim, typename Number>
568 vector_type &U, Number t)
569 {
570 // FIXME: refactor to avoid code duplication with step_erk_33
571
572#ifdef DEBUG_OUTPUT
573 std::cout << "TimeIntegrator<dim, Number>::step_strang_erk_33_cn()"
574 << std::endl;
575#endif
576
577 /* First explicit ERK(3,3,1) step with final result in U_[2]: */
578
579 Number tau = hyperbolic_module_->template step<0>(
580 /*input*/ U, {}, {}, {}, U_[0], precomputed_[0]);
581 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
582
583 hyperbolic_module_->template step<1>(U_[0],
584 {{/*input*/ U}},
585 {{precomputed_[0]}},
586 {{Number(-1.)}},
587 U_[1],
588 precomputed_[1],
589 tau);
590 hyperbolic_module_->apply_boundary_conditions(U_[1], t + 2. * tau);
591
592 hyperbolic_module_->template step<2>(U_[1],
593 {{/*input*/ U, U_[0]}},
594 {{precomputed_[0], precomputed_[1]}},
595 {{Number(0.75), Number(-2.)}},
596 U_[2],
597 precomputed_[2],
598 tau);
599 hyperbolic_module_->apply_boundary_conditions(U_[2], t + 3. * tau);
600
601 /* Implicit Crank-Nicolson step with final result in U_[3]: */
602
603 parabolic_module_->crank_nicolson_step(U_[2], t, U_[3], 6.0 * tau);
604
605 /* First explicit SSPRK 3 step with final result in U_[2]: */
606
607 hyperbolic_module_->template step<0>(
608 /*intermediate*/ U_[3], {}, {}, {}, U_[0], precomputed_[0], tau);
609 hyperbolic_module_->apply_boundary_conditions(U_[0], t + 4. * tau);
610
611 hyperbolic_module_->template step<1>(U_[0],
612 {{/*intermediate*/ U_[3]}},
613 {{precomputed_[0]}},
614 {{Number(-1.)}},
615 U_[1],
616 precomputed_[1],
617 tau);
618 hyperbolic_module_->apply_boundary_conditions(U_[1], t + 5. * tau);
619
620 hyperbolic_module_->template step<2>(U_[1],
621 {{/*intermediate*/ U_[3], U_[0]}},
622 {{precomputed_[0], precomputed_[1]}},
623 {{Number(0.75), Number(-2.)}},
624 U_[2],
625 precomputed_[2],
626 tau);
627 hyperbolic_module_->apply_boundary_conditions(U_[2], t + 6. * tau);
628
629 U.swap(U_[2]);
630 return 6. * tau;
631 }
632
633
634 template <typename Description, int dim, typename Number>
636 vector_type &U, Number t)
637 {
638 // FIXME: refactor to avoid code duplication with step_erk_43
639
640#ifdef DEBUG_OUTPUT
641 std::cout << "TimeIntegrator<dim, Number>::step_strang_erk_43_cn()"
642 << std::endl;
643#endif
644
645 /* First explicit ERK(4,3,1) step with final result in U_[3]: */
646
647 /* Step 1: U1 <- {U, 1} at time t + tau */
648 Number tau = hyperbolic_module_->template step<0>(
649 /*input*/ U, {}, {}, {}, U_[0], precomputed_[0]);
650 hyperbolic_module_->apply_boundary_conditions(U_[0], t + tau);
651
652 /* Step 2: U2 <- {U1, 2} and {U, -1} at time t + 2 tau */
653 hyperbolic_module_->template step<1>(U_[0],
654 {{/*input*/ U}},
655 {{precomputed_[0]}},
656 {{Number(-1.)}},
657 U_[1],
658 precomputed_[1],
659 tau);
660 hyperbolic_module_->apply_boundary_conditions(U_[1], t + 2. * tau);
661
662 /* Step 3: U3 <- {U2, 2} and {U1, -1} at time t + 3 tau */
663 hyperbolic_module_->template step<1>(U_[1],
664 {{U_[0]}},
665 {{precomputed_[1]}},
666 {{Number(-1.)}},
667 U_[2],
668 precomputed_[2],
669 tau);
670 hyperbolic_module_->apply_boundary_conditions(U_[2], t + 3. * tau);
671
672 /* Step 4: U4 <- {U3, 8/3} and {U2,-10/3} and {U1, 5/3} at time t + 4 tau */
673 hyperbolic_module_->template step<2>(U_[2],
674 {{U_[0], U_[1]}},
675 {{precomputed_[1], precomputed_[2]}},
676 {{Number(5. / 3.), Number(-10. / 3.)}},
677 U_[3],
678 precomputed_[3],
679 tau);
680 hyperbolic_module_->apply_boundary_conditions(U_[3], t + 4. * tau);
681
682 /* Implicit Crank-Nicolson step with final result in U_[2]: */
683
684 parabolic_module_->crank_nicolson_step(U_[3], t, U_[2], 8.0 * tau);
685
686 /* First explicit SSPRK 3 step with final result in U_[3]: */
687
688 /* Step 1: U1 <- {U, 1} at time t + tau */
689 hyperbolic_module_->template step<0>(
690 /*intermediate*/ U_[2], {}, {}, {}, U_[0], precomputed_[0], tau);
691 hyperbolic_module_->apply_boundary_conditions(U_[0], t + 5. * tau);
692
693 /* Step 2: U2 <- {U1, 2} and {U, -1} at time t + 2 tau */
694 hyperbolic_module_->template step<1>(U_[0],
695 {{/*intermediate*/ U_[2]}},
696 {{precomputed_[0]}},
697 {{Number(-1.)}},
698 U_[1],
699 precomputed_[1],
700 tau);
701 hyperbolic_module_->apply_boundary_conditions(U_[1], t + 6. * tau);
702
703 /* Step 3: U3 <- {U2, 2} and {U1, -1} at time t + 3 tau */
704 hyperbolic_module_->template step<1>(U_[1],
705 {{U_[0]}},
706 {{precomputed_[1]}},
707 {{Number(-1.)}},
708 U_[2],
709 precomputed_[2],
710 tau);
711 hyperbolic_module_->apply_boundary_conditions(U_[2], t + 7. * tau);
712
713 /* Step 4: U4 <- {U3, 8/3} and {U2,-10/3} and {U1, 5/3} at time t + 4 tau */
714 hyperbolic_module_->template step<2>(U_[2],
715 {{U_[0], U_[1]}},
716 {{precomputed_[1], precomputed_[2]}},
717 {{Number(5. / 3.), Number(-10. / 3.)}},
718 U_[3],
719 precomputed_[3],
720 tau);
721 hyperbolic_module_->apply_boundary_conditions(U_[3], t + 8. * tau);
722
723 U.swap(U_[3]);
724 return 8. * tau;
725 }
726
727} /* namespace ryujin */
Number step_erk_22(vector_type &U, Number t)
TimeIntegrator(const MPI_Comm &mpi_communicator, std::map< std::string, dealii::Timer > &computing_timer, 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_erk_11(vector_type &U, Number t)
Number step_erk_33(vector_type &U, Number t)
Number step(vector_type &U, Number t)
Number step_strang_erk_43_cn(vector_type &U, Number t)
Number step_strang_ssprk_33_cn(vector_type &U, Number t)
Number step_erk_54(vector_type &U, Number t)
Number step_ssprk_22(vector_type &U, Number t)
Number step_erk_43(vector_type &U, Number t)
Number step_ssprk_33(vector_type &U, Number t)
Number step_strang_erk_33_cn(vector_type &U, Number t)