ryujin 2.1.1 revision 6dc06e5864abd5d99e5d7ab641dbe621936411d9
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 template <typename StateVector, typename Number>
15 void
16 sadd(StateVector &dst, const Number s, const Number b, const StateVector &src)
17 {
18 auto &dst_U = std::get<0>(dst);
19 auto &src_U = std::get<0>(src);
20 dst_U.sadd(s, b, src_U);
21
22 auto &dst_V = std::get<2>(dst);
23 auto &src_V = std::get<2>(src);
24 dst_V.sadd(s, b, src_V);
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_->crank_nicolson_step(temp_[0], t, temp_[2], 2.0 * tau);
543
544 /* Second SSPRK 3 step with final result in temp_[0]: */
545
546 hyperbolic_module_->prepare_state_vector( temp_[2], t + 1.0 * tau);
547 hyperbolic_module_->template step<0>( temp_[2], {}, {}, temp_[0], tau);
548
549 hyperbolic_module_->prepare_state_vector(temp_[0], t + 2.0 * tau);
550 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
551 sadd(temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), temp_[2]);
552
553 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.5 * tau);
554 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
555 sadd(temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), temp_[2]);
556
557 state_vector.swap(temp_[0]);
558 return 2. * tau;
559 }
560
561
562 template <typename Description, int dim, typename Number>
563 Number TimeIntegrator<Description, dim, Number>::step_strang_erk_33_cn(
564 StateVector &state_vector, Number t, Number tau_max)
565 {
566 // FIXME: refactor to avoid code duplication with step_erk_33
567
568#ifdef DEBUG_OUTPUT
569 std::cout << "TimeIntegrator<dim, Number>::step_strang_erk_33_cn()"
570 << std::endl;
571#endif
572
573 /* First explicit ERK(3,3,1) step with final result in temp_[2]: */
574
575 hyperbolic_module_->prepare_state_vector(state_vector, t);
576 Number tau = hyperbolic_module_->template step<0>(
577 state_vector, {}, {}, temp_[0], Number(0.), tau_max / 6.);
578
579 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
580 hyperbolic_module_->template step<1>(
581 temp_[0], {{ state_vector}}, {{Number(-1.)}}, temp_[1], tau);
582
583 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
584 hyperbolic_module_->template step<2>(temp_[1],
585 {{ state_vector, temp_[0]}},
586 {{Number(0.75), Number(-2.)}},
587 temp_[2],
588 tau);
589
590 /* Implicit Crank-Nicolson step with final result in temp_[3]: */
591
592 parabolic_module_->crank_nicolson_step(temp_[2], t, temp_[3], 6.0 * tau);
593
594 /* Second explicit ERK(3,3,1) 3 step with final result in temp_[2]: */
595
596 hyperbolic_module_->prepare_state_vector(temp_[3], t + 3.0 * tau);
597 hyperbolic_module_->template step<0>(
598 temp_[3], {}, {}, temp_[0], tau);
599
600 hyperbolic_module_->prepare_state_vector(temp_[0], t + 4.0 * tau);
601 hyperbolic_module_->template step<1>(
602 temp_[0], {{ temp_[3]}}, {{Number(-1.)}}, temp_[1], tau);
603
604 hyperbolic_module_->prepare_state_vector(temp_[1], t + 5.0 * tau);
605 hyperbolic_module_->template step<2>(temp_[1],
606 {{ temp_[3], temp_[0]}},
607 {{Number(0.75), Number(-2.)}},
608 temp_[2],
609 tau);
610
611 state_vector.swap(temp_[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, Number tau_max)
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 temp_[3]: */
628
629 hyperbolic_module_->prepare_state_vector(state_vector, t);
630 Number tau = hyperbolic_module_->template step<0>(
631 state_vector, {}, {}, temp_[0], Number(0.), tau_max / 8.);
632
633 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
634 hyperbolic_module_->template step<1>(
635 temp_[0], {{ state_vector}}, {{Number(-1.)}}, temp_[1], tau);
636
637 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
638 hyperbolic_module_->template step<1>(
639 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
640
641 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
642 hyperbolic_module_->template step<2>(temp_[2],
643 {{temp_[0], temp_[1]}},
644 {{Number(5. / 3.), Number(-10. / 3.)}},
645 temp_[3],
646 tau);
647
648 /* Implicit Crank-Nicolson step with final result in temp_[2]: */
649
650 parabolic_module_->crank_nicolson_step(temp_[3], t, temp_[2], 8.0 * tau);
651
652 /* Second explicit ERK(4,3,1) step with final result in temp_[3]: */
653
654 hyperbolic_module_->prepare_state_vector(temp_[2], t + 4.0 * tau);
655 hyperbolic_module_->template step<0>(
656 temp_[2], {}, {}, temp_[0], tau);
657
658 hyperbolic_module_->prepare_state_vector(temp_[0], t + 5.0 * tau);
659 hyperbolic_module_->template step<1>(
660 temp_[0], {{ temp_[2]}}, {{Number(-1.)}}, temp_[1], tau);
661
662 hyperbolic_module_->prepare_state_vector(temp_[1], t + 6.0 * tau);
663 hyperbolic_module_->template step<1>(
664 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
665
666 hyperbolic_module_->prepare_state_vector(temp_[2], t + 7.0 * tau);
667 hyperbolic_module_->template step<2>(temp_[2],
668 {{temp_[0], temp_[1]}},
669 {{Number(5. / 3.), Number(-10. / 3.)}},
670 temp_[3],
671 tau);
672
673 state_vector.swap(temp_[3]);
674 return 8. * tau;
675 }
676
677 template <typename Description, int dim, typename Number>
678 Number TimeIntegrator<Description, dim, Number>::step_imex_11(
679 StateVector &state_vector, Number t, Number tau_max)
680 {
681#ifdef DEBUG_OUTPUT
682 std::cout << "TimeIntegrator<dim, Number>::step_imex_11()" << std::endl;
683#endif
684
685 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
686 hyperbolic_module_->prepare_state_vector(state_vector, t);
687 Number tau = hyperbolic_module_->template step<0>(
688 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
689
690 /* Implicit step 1: T1 <- {T0, 1} at time t -> t + tau */
691 parabolic_module_->template backward_euler_step<0>(
692 temp_[0], t, {}, {}, temp_[1], 1.0 * tau);
693
694 state_vector.swap(temp_[1]);
695 return tau;
696 }
697
698 template <typename Description, int dim, typename Number>
699 Number TimeIntegrator<Description, dim, Number>::step_imex_22(
700 StateVector &state_vector, Number t, Number tau_max)
701 {
702#ifdef DEBUG_OUTPUT
703 std::cout << "TimeIntegrator<dim, Number>::step_imex_22()" << std::endl;
704#endif
705
706 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
707 hyperbolic_module_->prepare_state_vector(state_vector, t);
708 Number tau = hyperbolic_module_->template step<0>(
709 state_vector, {}, {}, temp_[0], Number(0.), tau_max / 2.);
710
711 /* Implicit step 1: T1 <- {T0, 1} at time t -> t + tau */
712 parabolic_module_->template backward_euler_step<0>(
713 temp_[0], t, {}, {}, temp_[1], tau);
714
715 /* Explicit step 2: T2 <- {T1, 2} and {U_old, -1} at t + tau -> t + 2 tau */
716 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.0 * tau);
717 hyperbolic_module_->template step<1>(
718 temp_[1], {{state_vector}}, {{Number(-1.)}}, temp_[2], tau);
719
720 /* Implicit step 2: T3 <- {T2, 0} and {U_old, 1} at t + tau -> t + 2 tau */
721 parabolic_module_->template backward_euler_step<1>(temp_[2],
722 t + 1.0 * tau,
723 {{state_vector}},
724 {{Number(1.)}},
725 temp_[3],
726 tau);
727
728 state_vector.swap(temp_[3]);
729 return 2. * tau;
730 }
731
732 template <typename Description, int dim, typename Number>
733 Number TimeIntegrator<Description, dim, Number>::step_imex_33(
734 StateVector &state_vector, Number t, Number tau_max)
735 {
736#ifdef DEBUG_OUTPUT
737 std::cout << "TimeIntegrator<dim, Number>::step_imex_33()" << std::endl;
738#endif
739
740 /* IMEX(3, 3; 1), see @cite ErnGuermond2023, Sec. 4.3. */
741
742 const Number gamma = 0.5 + 0.5 * std::numbers::inv_sqrt3;
743
744 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
745 hyperbolic_module_->prepare_state_vector(state_vector, t);
746 Number tau = hyperbolic_module_->template step<0>(
747 state_vector, {}, {}, temp_[0], Number(0.), tau_max / 3.);
748
749 /* Implicit step 1: T1 <- {U_old, 1 - 3*gamma} at time t -> t + tau */
750 parabolic_module_->template backward_euler_step<1>(
751 temp_[0],
752 t,
753 {{state_vector}},
754 {{Number(1. - 3. * gamma)}},
755 temp_[1],
756 tau);
757
758 /* Explicit step 2: T2 <- {U_old, -1} and {T1, 2} at time t -> t + 2 tau */
759 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.0 * tau);
760 hyperbolic_module_->template step<1>(
761 temp_[1], {{state_vector}}, {{Number(-1.)}}, temp_[2], tau);
762
763 /*
764 * Implicit step 2:
765 * T3 <- {U_old, 6*gamma-1} and {T1, 2-9*gamma} at t -> t + * 2 tau
766 */
767 parabolic_module_->template backward_euler_step<2>(
768 temp_[2],
769 t + tau,
770 {{state_vector, temp_[1]}},
771 {{Number(6. * gamma - 1.), Number(2. - 9 * gamma)}},
772 temp_[3],
773 tau);
774
775 /* Explicit step 3: T4 <- {U_old, 3 / 4} and {T1, -2} at t -> t + 3 tau */
776 hyperbolic_module_->prepare_state_vector(temp_[3], t + 2. * tau);
777 hyperbolic_module_->template step<2>(temp_[3],
778 {{state_vector, temp_[1]}},
779 {{Number(0.75), Number(-2.)}},
780 temp_[4],
781 tau);
782
783 /* Implicit step 3: */
784 parabolic_module_->template backward_euler_step<3>(
785 temp_[4],
786 t + 2. * tau,
787 {{state_vector, temp_[1], temp_[3]}},
788 {{Number(0.75 - 3. * gamma),
789 Number(6. * gamma - 2.),
790 Number(9. / 4. - 3. * gamma)}},
791 temp_[5],
792 tau);
793
794 state_vector.swap(temp_[5]);
795 return 3. * tau;
796 }
797
798} /* namespace ryujin */
ryujin::TimeIntegrator::step_imex_33
Number step_imex_33(StateVector &state_vector, Number t, Number tau_max)
Definition: time_integrator.template.h:733
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:678
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:617
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:699
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:563
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:16
time_integrator.h
Generated by doxygen 1.9.4