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