ryujin 2.1.1 revision 0348cbb53a3e4b1da2a4c037e81f88f2d21ce219
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 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
61 if (ParabolicSystem::is_identity)
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>
75 void TimeIntegrator<Description, dim, Number>::prepare()
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_) {
84 case TimeSteppingScheme::ssprk_22:
85 temp_.resize(2);
86 efficiency_ = 1.;
87 break;
88 case TimeSteppingScheme::ssprk_33:
89 temp_.resize(2);
90 efficiency_ = 1.;
91 break;
92 case TimeSteppingScheme::erk_11:
93 temp_.resize(1);
94 efficiency_ = 1.;
95 break;
96 case TimeSteppingScheme::erk_22:
97 temp_.resize(2);
98 efficiency_ = 2.;
99 break;
100 case TimeSteppingScheme::erk_33:
101 temp_.resize(3);
102 efficiency_ = 3.;
103 break;
104 case TimeSteppingScheme::erk_43:
105 temp_.resize(4);
106 efficiency_ = 4.;
107 break;
108 case TimeSteppingScheme::erk_54:
109 temp_.resize(5);
110 efficiency_ = 5.;
111 break;
112 case TimeSteppingScheme::strang_ssprk_33_cn:
113 temp_.resize(3);
114 efficiency_ = 2.;
115 break;
116 case TimeSteppingScheme::strang_erk_33_cn:
117 temp_.resize(4);
118 efficiency_ = 6.;
119 break;
120 case TimeSteppingScheme::strang_erk_43_cn:
121 temp_.resize(4);
122 efficiency_ = 8.;
123 break;
124 case TimeSteppingScheme::imex_11:
125 temp_.resize(2);
126 efficiency_ = 1.;
127 break;
128 case TimeSteppingScheme::imex_22:
129 temp_.resize(4);
130 efficiency_ = 2.;
131 break;
132 case TimeSteppingScheme::imex_33:
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_) {
158 case TimeSteppingScheme::ssprk_22:
159 [[fallthrough]];
160 case TimeSteppingScheme::ssprk_33:
161 [[fallthrough]];
162 case TimeSteppingScheme::erk_11:
163 [[fallthrough]];
164 case TimeSteppingScheme::erk_22:
165 [[fallthrough]];
166 case TimeSteppingScheme::erk_33:
167 [[fallthrough]];
168 case TimeSteppingScheme::erk_43:
169 [[fallthrough]];
170 case TimeSteppingScheme::erk_54: {
171 AssertThrow(
172 ParabolicSystem::is_identity,
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 }
179 case TimeSteppingScheme::imex_11:
180 [[fallthrough]];
181 case TimeSteppingScheme::imex_22:
182 [[fallthrough]];
183 case TimeSteppingScheme::imex_33:
184 [[fallthrough]];
185 case TimeSteppingScheme::strang_ssprk_33_cn:
186 [[fallthrough]];
187 case TimeSteppingScheme::strang_erk_33_cn:
188 [[fallthrough]];
189 case TimeSteppingScheme::strang_erk_43_cn: {
190 AssertThrow(
191 !ParabolicSystem::is_identity,
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>
207 Number TimeIntegrator<Description, dim, Number>::step(
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_) {
219 case TimeSteppingScheme::ssprk_22:
220 return step_ssprk_22(state_vector, t, tau_max);
221 case TimeSteppingScheme::ssprk_33:
222 return step_ssprk_33(state_vector, t, tau_max);
223 case TimeSteppingScheme::erk_11:
224 return step_erk_11(state_vector, t, tau_max);
225 case TimeSteppingScheme::erk_22:
226 return step_erk_22(state_vector, t, tau_max);
227 case TimeSteppingScheme::erk_33:
228 return step_erk_33(state_vector, t, tau_max);
229 case TimeSteppingScheme::erk_43:
230 return step_erk_43(state_vector, t, tau_max);
231 case TimeSteppingScheme::erk_54:
232 return step_erk_54(state_vector, t, tau_max);
233 case TimeSteppingScheme::strang_ssprk_33_cn:
234 return step_strang_ssprk_33_cn(state_vector, t, tau_max);
235 case TimeSteppingScheme::strang_erk_33_cn:
236 return step_strang_erk_33_cn(state_vector, t, tau_max);
237 case TimeSteppingScheme::strang_erk_43_cn:
238 return step_strang_erk_43_cn(state_vector, t, tau_max);
239 case TimeSteppingScheme::imex_11:
240 return step_imex_11(state_vector, t, tau_max);
241 case TimeSteppingScheme::imex_22:
242 return step_imex_22(state_vector, t, tau_max);
243 case TimeSteppingScheme::imex_33:
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_ =
252 IDViolationStrategy::raise_exception;
253 parabolic_module_->id_violation_strategy_ =
254 IDViolationStrategy::raise_exception;
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>
279 Number TimeIntegrator<Description, dim, Number>::step_ssprk_22(
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>
302 Number TimeIntegrator<Description, dim, Number>::step_ssprk_33(
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>
332 Number TimeIntegrator<Description, dim, Number>::step_erk_11(
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);
343
344 state_vector.swap(temp_[0]);
345 return tau;
346 }
347
348
349 template <typename Description, int dim, typename Number>
350 Number TimeIntegrator<Description, dim, Number>::step_erk_22(
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>
373 Number TimeIntegrator<Description, dim, Number>::step_erk_33(
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);
400
401 state_vector.swap(temp_[2]);
402 return 3. * tau;
403 }
404
405
406 template <typename Description, int dim, typename Number>
407 Number TimeIntegrator<Description, dim, Number>::step_erk_43(
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>
446 Number TimeIntegrator<Description, dim, Number>::step_erk_54(
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>
516 Number TimeIntegrator<Description, dim, Number>::step_strang_ssprk_33_cn(
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>
564 Number TimeIntegrator<Description, dim, Number>::step_strang_erk_33_cn(
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>
620 Number TimeIntegrator<Description, dim, Number>::step_strang_erk_43_cn(
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>
683 Number TimeIntegrator<Description, dim, Number>::step_imex_11(
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>
704 Number TimeIntegrator<Description, dim, Number>::step_imex_22(
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>
737 Number TimeIntegrator<Description, dim, Number>::step_imex_33(
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 */
ryujin::TimeIntegrator::step_imex_33
Number step_imex_33(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:737
ryujin::TimeIntegrator::TimeIntegrator
TimeIntegrator(const MPIEnsemble &mpi_ensemble, const OfflineData< dim, Number > &offline_data, const HyperbolicModule< Description, dim, Number > &hyperbolic_module, const ParabolicModule< Description, dim, Number > &parabolic_module, const std::string &subsection="/TimeIntegrator")
Definition: time_integrator.template.h:29
ryujin::TimeIntegrator::step_strang_ssprk_33_cn
Number step_strang_ssprk_33_cn(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:516
ryujin::TimeIntegrator::step_erk_22
Number step_erk_22(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:350
ryujin::TimeIntegrator::step_imex_11
Number step_imex_11(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:683
ryujin::TimeIntegrator::step_strang_erk_43_cn
Number step_strang_erk_43_cn(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:620
ryujin::TimeIntegrator::step_erk_54
Number step_erk_54(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:446
ryujin::TimeIntegrator::step_erk_33
Number step_erk_33(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:373
ryujin::TimeIntegrator::step_erk_43
Number step_erk_43(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:407
ryujin::TimeIntegrator::step_imex_22
Number step_imex_22(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:704
ryujin::TimeIntegrator::StateVector
typename View::StateVector StateVector
Definition: time_integrator.h:257
ryujin::TimeIntegrator::step_erk_11
Number step_erk_11(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:332
ryujin::TimeIntegrator::step_ssprk_22
Number step_ssprk_22(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:279
ryujin::TimeIntegrator::step_ssprk_33
Number step_ssprk_33(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:302
ryujin::TimeIntegrator::step
Number step(StateVector &state_vector, Number t, Number t_final=std::numeric_limits< Number >::max())
Definition: time_integrator.template.h:207
ryujin::TimeIntegrator::step_strang_erk_33_cn
Number step_strang_erk_33_cn(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:564
ryujin::IDViolationStrategy::warn
@ warn
ryujin::IDViolationStrategy::raise_exception
@ raise_exception
ryujin::TimeSteppingScheme::erk_22
@ erk_22
ryujin::TimeSteppingScheme::erk_33
@ erk_33
ryujin::TimeSteppingScheme::imex_11
@ imex_11
ryujin::TimeSteppingScheme::erk_43
@ erk_43
ryujin::TimeSteppingScheme::erk_11
@ erk_11
ryujin::TimeSteppingScheme::ssprk_33
@ ssprk_33
ryujin::TimeSteppingScheme::strang_erk_43_cn
@ strang_erk_43_cn
ryujin::TimeSteppingScheme::ssprk_22
@ ssprk_22
ryujin::TimeSteppingScheme::strang_erk_33_cn
@ strang_erk_33_cn
ryujin::TimeSteppingScheme::imex_22
@ imex_22
ryujin::TimeSteppingScheme::strang_ssprk_33_cn
@ strang_ssprk_33_cn
ryujin::TimeSteppingScheme::imex_33
@ imex_33
ryujin::TimeSteppingScheme::erk_54
@ erk_54
ryujin::CFLRecoveryStrategy::bang_bang_control
@ bang_bang_control
ryujin::CFLRecoveryStrategy::none
@ none
ryujin::Vectors::StateVector
std::tuple< MultiComponentVector< Number, problem_dim >, MultiComponentVector< Number, prec_dim >, BlockVector< Number > > StateVector
Definition: state_vector.h:51
ryujin::sadd
void sadd(StateVector &dst, const Number s, const Number b, const StateVector &src)
Definition: time_integrator.template.h:20
time_integrator.h
Generated by doxygen 1.9.4