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