ryujin 2.1.1 revision b6c23cccd75d6a01937bdcd5d1dbe2aae9b59c42
Loading...
Searching...
No Matches
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 - 2025 by the ryujin authors
4//
5
6#pragma once
7
8#include "time_integrator.h"
9
10#include <utility>
11
12namespace ryujin
13{
14 using namespace dealii;
15
16 template <typename StateVector, typename Number>
17 void
18 sadd(StateVector &dst, const Number s, const Number b, const StateVector &src)
19 {
20 /* Perform sadd() on the default memory space: */
21
22 using MemorySpace = selected_memory_space_t;
23
24 const auto dst_U_view = std::get<0>(dst).template view<MemorySpace>();
25 const auto src_U_view =
26 std::get<0>(std::as_const(src)).template view<MemorySpace>();
27 dst_U_view.zero_out_ghost_values();
28 dst_U_view.sadd(s, b, src_U_view);
29
30 /*
31 * FIXME: At the moment the parabolic state resides only on the host memory
32 * space. Update the following to view<MemorySpace>() when we have
33 * refactored the parabolic modules.
34 */
35
36 auto &dst_V = std::get<2>(dst);
37 const auto &src_V = std::get<2>(src);
38 AssertDimension(dst_V.size(), src_V.size());
39
40 for (std::size_t k = 0; k < dst_V.size(); ++k) {
41 const auto dst_view = dst_V[k].view();
42 const auto src_view = src_V[k].view();
43 dst_view.zero_out_ghost_values();
44 dst_view.sadd(s, b, src_view);
45 }
46 }
47
48
49 template <typename Description, int dim, typename Number>
50 TimeIntegrator<Description, dim, Number>::TimeIntegrator(
51 const MPIEnsemble &mpi_ensemble,
52 const OfflineData<dim, Number> &offline_data,
53 const HyperbolicModule<Description, dim, Number> &hyperbolic_module,
54 const ParabolicModule<Description, dim, Number> &parabolic_module,
55 const std::string &subsection /*= "TimeIntegrator"*/)
56 : ParameterAcceptor(subsection)
57 , mpi_ensemble_(mpi_ensemble)
58 , offline_data_(&offline_data)
59 , hyperbolic_module_(&hyperbolic_module)
60 , parabolic_module_(&parabolic_module)
61 {
62 cfl_min_ = Number(0.45);
63 add_parameter(
64 "cfl min",
65 cfl_min_,
66 "Minimal admissible relative CFL constant. How this parameter is used "
67 "depends on the chosen CFL recovery strategy");
68
69 cfl_max_ = Number(0.90);
70 add_parameter(
71 "cfl max",
72 cfl_max_,
73 "Maximal admissible relative CFL constant. How this parameter is used "
74 "depends on the chosen CFL recovery strategy");
75
76 cfl_recovery_strategy_ = CFLRecoveryStrategy::cruise_control;
77 add_parameter("cfl recovery strategy",
78 cfl_recovery_strategy_,
79 "CFL/invariant domain violation recovery strategy: none, "
80 "bang bang control, cruise control");
81
82 acceptable_tau_max_ratio_ = Number(2.0);
83 add_parameter("acceptable tau_max ratio",
84 acceptable_tau_max_ratio_,
85 "Maximal acceptable discrepancy between computed tau_max of "
86 "a (sub)step and enforced time-step size tau. If the ratio "
87 "is violated then a restart will be singnalled.");
88
89 tau_max_ = std::numeric_limits<Number>::max();
90 add_parameter("tau_max", tau_max_, "Largest time step size allowed.");
91
92 if (ParabolicSystem::is_identity)
93 time_stepping_scheme_ = TimeSteppingScheme::erk_33;
94 else
95 time_stepping_scheme_ = TimeSteppingScheme::strang_erk_33_cn;
96 add_parameter("time stepping scheme",
97 time_stepping_scheme_,
98 "Time stepping scheme: ssprk 22, ssprk 33, erk 11, erk 22, "
99 "erk 33, erk 43, erk "
100 "54, strang ssprk 33 cn, strang erk 33 cn, strang erk 43 cn, "
101 "imex 11, imex 22, imex 33");
102 }
103
104
105 template <typename Description, int dim, typename Number>
106 void TimeIntegrator<Description, dim, Number>::prepare()
107 {
108#ifdef DEBUG_OUTPUT
109 std::cout << "TimeIntegrator<dim, Number>::prepare()" << std::endl;
110#endif
111
112 /* Resize temporary storage to appropriate sizes: */
113
114 switch (time_stepping_scheme_) {
115 case TimeSteppingScheme::ssprk_22:
116 temp_.resize(2);
117 efficiency_ = 1.;
118 break;
119 case TimeSteppingScheme::ssprk_33:
120 temp_.resize(2);
121 efficiency_ = 1.;
122 break;
123 case TimeSteppingScheme::erk_11:
124 temp_.resize(1);
125 efficiency_ = 1.;
126 break;
127 case TimeSteppingScheme::erk_22:
128 temp_.resize(2);
129 efficiency_ = 2.;
130 break;
131 case TimeSteppingScheme::erk_33:
132 temp_.resize(3);
133 efficiency_ = 3.;
134 break;
135 case TimeSteppingScheme::erk_43:
136 temp_.resize(4);
137 efficiency_ = 4.;
138 break;
139 case TimeSteppingScheme::erk_54:
140 temp_.resize(5);
141 efficiency_ = 5.;
142 break;
143 case TimeSteppingScheme::strang_ssprk_33_cn:
144 temp_.resize(3);
145 efficiency_ = 2.;
146 break;
147 case TimeSteppingScheme::strang_erk_33_cn:
148 temp_.resize(4);
149 efficiency_ = 6.;
150 break;
151 case TimeSteppingScheme::strang_erk_43_cn:
152 temp_.resize(4);
153 efficiency_ = 8.;
154 break;
155 case TimeSteppingScheme::imex_11:
156 temp_.resize(2);
157 efficiency_ = 1.;
158 break;
159 case TimeSteppingScheme::imex_22:
160 temp_.resize(4);
161 efficiency_ = 2.;
162 break;
163 case TimeSteppingScheme::imex_33:
164 temp_.resize(6);
165 efficiency_ = 3.;
166 break;
167 }
168
169 /* Initialize temporary vectors: */
170
171 for (auto &it : temp_) {
172 hyperbolic_module_->reinit_state_vector(it);
173 parabolic_module_->reinit_state_vector(it);
174 }
175
176 /* Reset CFL to starting value, set maximal acceptable tau_max ratio: */
177
178 AssertThrow(cfl_min_ > 0., ExcMessage("cfl min must be a positive value"));
179 AssertThrow(cfl_max_ >= cfl_min_,
180 ExcMessage("cfl max must be greater than or equal to cfl min"));
181
182 AssertThrow(
183 acceptable_tau_max_ratio_ >= 1.0,
184 ExcMessage(
185 "acceptable tau_max ratio must be greater than or equal to 1."));
186
187 hyperbolic_module_->set_cfl(cfl_max_);
188 hyperbolic_module_->set_acceptable_tau_max_ratio(acceptable_tau_max_ratio_);
189
190 const auto check_whether_timestepping_makes_sense = [&]() {
191 /*
192 * Make sure the user selects an appropriate time-stepping scheme.
193 */
194
195 switch (time_stepping_scheme_) {
196 case TimeSteppingScheme::ssprk_22:
197 [[fallthrough]];
198 case TimeSteppingScheme::ssprk_33:
199 [[fallthrough]];
200 case TimeSteppingScheme::erk_11:
201 [[fallthrough]];
202 case TimeSteppingScheme::erk_22:
203 [[fallthrough]];
204 case TimeSteppingScheme::erk_33:
205 [[fallthrough]];
206 case TimeSteppingScheme::erk_43:
207 [[fallthrough]];
208 case TimeSteppingScheme::erk_54: {
209 AssertThrow(
210 ParabolicSystem::is_identity,
211 dealii::ExcMessage(
212 "The selected equation consists of a hyperbolic and nontrivial "
213 "parabolic subsystem and requires an IMEX timestepping "
214 "scheme such as »strang erk 33 cn«."));
215 break;
216 }
217 case TimeSteppingScheme::imex_11:
218 [[fallthrough]];
219 case TimeSteppingScheme::imex_22:
220 [[fallthrough]];
221 case TimeSteppingScheme::imex_33:
222 [[fallthrough]];
223 case TimeSteppingScheme::strang_ssprk_33_cn:
224 [[fallthrough]];
225 case TimeSteppingScheme::strang_erk_33_cn:
226 [[fallthrough]];
227 case TimeSteppingScheme::strang_erk_43_cn: {
228 AssertThrow(
229 !ParabolicSystem::is_identity,
230 dealii::ExcMessage(
231 "The selected equation has a trivial parabolic subsystem and "
232 "should not be run with an IMEX timestepping scheme."));
233 break;
234 }
235 }
236 };
237
238 check_whether_timestepping_makes_sense();
239 this->parse_parameters_call_back.connect(
240 check_whether_timestepping_makes_sense);
241 }
242
243
244 /*
245 * -------------------------------------------------------------------------
246 * Prepare state vector:
247 * -------------------------------------------------------------------------
248 */
249
250
251 template <typename Description, int dim, typename Number>
252 void TimeIntegrator<Description, dim, Number>::prepare_state_vector(
253 StateVector &state_vector, Number t) const
254 {
255 if (!ParabolicSystem::is_identity)
256 parabolic_module_->prepare_state_vector(state_vector, t);
257 hyperbolic_module_->prepare_state_vector(state_vector, t);
258 }
259
260
261 /*
262 * -------------------------------------------------------------------------
263 * High level step function implementing various CFLRecoveryStrategy
264 * -------------------------------------------------------------------------
265 */
266
267
268 template <typename Description, int dim, typename Number>
269 Number TimeIntegrator<Description, dim, Number>::step(
270 StateVector &state_vector,
271 Number t,
272 Number t_final /*=std::numeric_limits<Number>::max()*/)
273 {
274 Number tau_max =
275 std::min(tau_max_, t_final - t); /* enforces t <= t_final */
276
277#ifdef DEBUG_OUTPUT
278 std::cout << "TimeIntegrator<dim, Number>::step()" << std::endl;
279 std::cout << " enforcing tau_max <= " << tau_max << std::endl;
280#endif
281
282 const auto single_step = [&]() {
283 switch (time_stepping_scheme_) {
284 case TimeSteppingScheme::ssprk_22:
285 return step_ssprk_22(state_vector, t, tau_max);
286 case TimeSteppingScheme::ssprk_33:
287 return step_ssprk_33(state_vector, t, tau_max);
288 case TimeSteppingScheme::erk_11:
289 return step_erk_11(state_vector, t, tau_max);
290 case TimeSteppingScheme::erk_22:
291 return step_erk_22(state_vector, t, tau_max);
292 case TimeSteppingScheme::erk_33:
293 return step_erk_33(state_vector, t, tau_max);
294 case TimeSteppingScheme::erk_43:
295 return step_erk_43(state_vector, t, tau_max);
296 case TimeSteppingScheme::erk_54:
297 return step_erk_54(state_vector, t, tau_max);
298 case TimeSteppingScheme::strang_ssprk_33_cn:
299 return step_strang_ssprk_33_cn(state_vector, t, tau_max);
300 case TimeSteppingScheme::strang_erk_33_cn:
301 return step_strang_erk_33_cn(state_vector, t, tau_max);
302 case TimeSteppingScheme::strang_erk_43_cn:
303 return step_strang_erk_43_cn(state_vector, t, tau_max);
304 case TimeSteppingScheme::imex_11:
305 return step_imex_11(state_vector, t, tau_max);
306 case TimeSteppingScheme::imex_22:
307 return step_imex_22(state_vector, t, tau_max);
308 case TimeSteppingScheme::imex_33:
309 return step_imex_33(state_vector, t, tau_max);
310 default:
311 __builtin_unreachable();
312 }
313 };
314
315 if (cfl_recovery_strategy_ != CFLRecoveryStrategy::none) {
316 hyperbolic_module_->set_id_violation_strategy(
317 IDViolationStrategy::raise_exception);
318 parabolic_module_->set_id_violation_strategy(
319 IDViolationStrategy::raise_exception);
320 hyperbolic_module_->set_cfl(cfl_max_);
321 }
322
323 try {
324 return single_step();
325
326 } catch (const Restart &restart) {
327
328 AssertThrow(cfl_recovery_strategy_ != CFLRecoveryStrategy::none,
329 dealii::ExcInternalError());
330
331 hyperbolic_module_->set_id_violation_strategy(IDViolationStrategy::warn);
332 parabolic_module_->set_id_violation_strategy(IDViolationStrategy::warn);
333
334 if (cfl_recovery_strategy_ == CFLRecoveryStrategy::bang_bang_control) {
335 /* Retry with cfl_min instead of cfl_max: */
336#ifdef DEBUG_OUTPUT
337 std::cout
338 << " restart with bang bang control: setting cfl to cfl_min"
339 << std::endl;
340#endif
341 hyperbolic_module_->set_cfl(cfl_min_);
342 }
343
344 if (cfl_recovery_strategy_ == CFLRecoveryStrategy::cruise_control) {
345 /* Retry with the suggested tau_max: */
346#ifdef DEBUG_OUTPUT
347 std::cout
348 << " restart with cruise control: using suggested_tau_max"
349 << std::endl;
350#endif
351 //
352 // Multiply the suggested tau_max value with the efficiency.
353 //
354 // We have to account for the fact that the e Restart exception is
355 // thrown within a substep of the hyperbolic or parabolic module.
356 // This implies that the suggested_tau_max is computed for that
357 // particular substep and not for the full combined method (where
358 // tau_max can be larger). We thus multiply tau_max with the
359 // efficiency factor.
360 //
361 tau_max =
362 std::min(tau_max, efficiency_ * Number(restart.suggested_tau_max));
363 }
364
365 return single_step();
366 }
367 }
368
369
370 /*
371 * -------------------------------------------------------------------------
372 * Concrete implementation of ERK / IMEX time stepping strategies.
373 * -------------------------------------------------------------------------
374 */
375
376
377 template <typename Description, int dim, typename Number>
378 Number TimeIntegrator<Description, dim, Number>::step_ssprk_22(
379 StateVector &state_vector, Number t, Number tau_max)
380 {
381 /* SSP-RK2, see @cite Shu1988, Eq. 2.15. */
382
383#ifdef DEBUG_OUTPUT
384 std::cout << "TimeIntegrator<dim, Number>::step_ssprk_22()" << std::endl;
385#endif
386
387 Assert(efficiency_ == 1., dealii::ExcInternalError());
388
389 /* Step 1: T0 = U_old + tau * L(U_old) at t -> t + tau */
390 Number tau = hyperbolic_module_->template step<0>(
391 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
392
393 /* Step 2: T1 = T0 + tau L(T0) at time t + tau -> t + 2*tau */
394 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
395 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
396
397 /* Step 2: convex combination: T1 = 1/2 U_old + 1/2 T1 at time t + tau */
398 sadd(temp_[1], Number(1.0 / 2.0), Number(1.0 / 2.0), state_vector);
399
400 state_vector.swap(temp_[1]);
401 return tau;
402 }
403
404
405 template <typename Description, int dim, typename Number>
406 Number TimeIntegrator<Description, dim, Number>::step_ssprk_33(
407 StateVector &state_vector, Number t, Number tau_max)
408 {
409 /* SSP-RK3, see @cite Shu1988, Eq. 2.18. */
410
411#ifdef DEBUG_OUTPUT
412 std::cout << "TimeIntegrator<dim, Number>::step_ssprk_33()" << std::endl;
413#endif
414
415 Assert(efficiency_ == 1., dealii::ExcInternalError());
416
417 /* Step 1: T0 = U_old + tau * L(U_old) at time t -> t + tau */
418 Number tau = hyperbolic_module_->template step<0>(
419 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
420
421 /* Step 2: T1 = T0 + tau L(T0) at time t + tau -> t + 2*tau */
422 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
423 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
424
425 /* Step 2: convex combination T1 = 3/4 U_old + 1/4 T1 at time t + 0.5*tau */
426 sadd(temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), state_vector);
427
428 /* Step 3: T0 = T1 + tau L(T1) at time t + 0.5*tau -> t + 1.5*tau */
429 hyperbolic_module_->prepare_state_vector(temp_[1], t + 0.5 * tau);
430 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
431
432 /* Step 3: convex combination: T0 = 1/3 U_old + 2/3 T0 at time t + tau */
433 sadd(temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), state_vector);
434
435 state_vector.swap(temp_[0]);
436 return tau;
437 }
438
439
440 template <typename Description, int dim, typename Number>
441 Number TimeIntegrator<Description, dim, Number>::step_erk_11(
442 StateVector &state_vector, Number /*t*/, Number tau_max)
443 {
444#ifdef DEBUG_OUTPUT
445 std::cout << "TimeIntegrator<dim, Number>::step_erk_11()" << std::endl;
446#endif
447
448 Assert(efficiency_ == 1., dealii::ExcInternalError());
449
450 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
451 Number tau = hyperbolic_module_->template step<0>(
452 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
453
454 state_vector.swap(temp_[0]);
455 return tau;
456 }
457
458
459 template <typename Description, int dim, typename Number>
460 Number TimeIntegrator<Description, dim, Number>::step_erk_22(
461 StateVector &state_vector, Number t, Number tau_max)
462 {
463#ifdef DEBUG_OUTPUT
464 std::cout << "TimeIntegrator<dim, Number>::step_erk_22()" << std::endl;
465#endif
466
467 Assert(efficiency_ == 2., dealii::ExcInternalError());
468
469 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
470 Number tau = hyperbolic_module_->template step<0>(
471 state_vector, {}, {}, temp_[0], Number(.0), tau_max / efficiency_);
472
473 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + tau -> t + 2*tau */
474 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
475 hyperbolic_module_->template step<1>(
476 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
477
478 state_vector.swap(temp_[1]);
479 return efficiency_ * tau;
480 }
481
482
483 template <typename Description, int dim, typename Number>
484 Number TimeIntegrator<Description, dim, Number>::step_erk_33(
485 StateVector &state_vector, Number t, Number tau_max)
486 {
487#ifdef DEBUG_OUTPUT
488 std::cout << "TimeIntegrator<dim, Number>::step_erk_33()" << std::endl;
489#endif
490
491 Assert(efficiency_ == 3., dealii::ExcInternalError());
492
493 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
494 Number tau = hyperbolic_module_->template step<0>(
495 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
496
497 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + 1*tau -> t + 2*tau */
498 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
499 hyperbolic_module_->template step<1>(
500 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
501
502 /*
503 * Step 3: T2 <- {T1, 9/4} and {T0, -2} and {U_old, 3/4}
504 * at time t + 2*tau -> t + 3*tau
505 */
506 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
507 hyperbolic_module_->template step<2>(temp_[1],
508 {{state_vector, temp_[0]}},
509 {{Number(0.75), Number(-2.)}},
510 temp_[2],
511 tau);
512
513 state_vector.swap(temp_[2]);
514 return efficiency_ * tau;
515 }
516
517
518 template <typename Description, int dim, typename Number>
519 Number TimeIntegrator<Description, dim, Number>::step_erk_43(
520 StateVector &state_vector, Number t, Number tau_max)
521 {
522#ifdef DEBUG_OUTPUT
523 std::cout << "TimeIntegrator<dim, Number>::step_erk_43()" << std::endl;
524#endif
525
526 Assert(efficiency_ == 4., dealii::ExcInternalError());
527
528 /* Step 1: T0 <- {U_old, 1} at time t -> t + tau */
529 Number tau = hyperbolic_module_->template step<0>(
530 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
531
532 /* Step 2: T1 <- {T0, 2} and {U_old, -1} at time t + 1*tau -> t + 2*tau */
533 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
534 hyperbolic_module_->template step<1>(
535 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
536
537 /* Step 3: T2 <- {T1, 2} and {T0, -1} at time t + 2*tau -> t + 3*tau */
538 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
539 hyperbolic_module_->template step<1>(
540 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
541
542 /*
543 * Step 4: T3 <- {T2, 8/3} and {T1,-10/3} and {T0, 5/3}
544 * at time t + 3*tau -> t + 4*tau
545 */
546 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
547 hyperbolic_module_->template step<2>(temp_[2],
548 {{temp_[0], temp_[1]}},
549 {{Number(5. / 3.), Number(-10. / 3.)}},
550 temp_[3],
551 tau);
552
553 state_vector.swap(temp_[3]);
554 return efficiency_ * tau;
555 }
556
557
558 template <typename Description, int dim, typename Number>
559 Number TimeIntegrator<Description, dim, Number>::step_erk_54(
560 StateVector &state_vector, Number t, Number tau_max)
561 {
562#ifdef DEBUG_OUTPUT
563 std::cout << "TimeIntegrator<dim, Number>::step_erk_54()" << std::endl;
564#endif
565
566 Assert(efficiency_ == 5., dealii::ExcInternalError());
567
568 constexpr Number c = 0.2; /* equidistant c_i */
569 constexpr Number a_21 = +0.2;
570 constexpr Number a_31 = +0.26075582269554909;
571 constexpr Number a_32 = +0.13924417730445096;
572 constexpr Number a_41 = -0.25856517872570289;
573 constexpr Number a_42 = +0.91136274166280729;
574 constexpr Number a_43 = -0.05279756293710430;
575 constexpr Number a_51 = +0.21623276431503774;
576 constexpr Number a_52 = +0.51534223099602405;
577 constexpr Number a_53 = -0.81662794199265554;
578 constexpr Number a_54 = +0.88505294668159373;
579 constexpr Number a_61 = -0.10511678454691901; /* aka b_1 */
580 constexpr Number a_62 = +0.87880047152100838; /* aka b_2 */
581 constexpr Number a_63 = -0.58903404061484477; /* aka b_3 */
582 constexpr Number a_64 = +0.46213380485434047; /* aka b_4 */
583 constexpr Number a_65 [[maybe_unused]] = +0.35321654878641495; /* aka b_5 */
584
585 /* Step 1: at time t -> t + 1*tau */
586 Number tau = hyperbolic_module_->template step<0>(
587 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
588
589 /* Step 2: at time t + 1*tau -> t + 2*tau */
590 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
591 hyperbolic_module_->template step<1>(
592 temp_[0], {{state_vector}}, {{(a_31 - a_21) / c}}, temp_[1], tau);
593
594 /* Step 3: at time t + 2*tau -> t + 3*tau */
595 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
596 hyperbolic_module_->template step<2>(
597 temp_[1],
598 {{state_vector, temp_[0]}},
599 {{(a_41 - a_31) / c, (a_42 - a_32) / c}},
600 temp_[2],
601 tau);
602
603 /* Step 4: at time t + 3*tau -> t + 4*tau */
604 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
605 hyperbolic_module_->template step<3>(
606 temp_[2],
607 {{state_vector, temp_[0], temp_[1]}},
608 {{(a_51 - a_41) / c, (a_52 - a_42) / c, (a_53 - a_43) / c}},
609 temp_[3],
610 tau);
611
612 /* Step 5: at time t + 4*tau -> t + 5*tau */
613 hyperbolic_module_->prepare_state_vector(temp_[3], t + 4.0 * tau);
614 hyperbolic_module_->template step<4>(
615 temp_[3],
616 {{state_vector, temp_[0], temp_[1], temp_[2]}},
617 {{(a_61 - a_51) / c,
618 (a_62 - a_52) / c,
619 (a_63 - a_53) / c,
620 (a_64 - a_54) / c}},
621 temp_[4],
622 tau);
623
624 state_vector.swap(temp_[4]);
625 return efficiency_ * tau;
626 }
627
628
629 template <typename Description, int dim, typename Number>
630 Number TimeIntegrator<Description, dim, Number>::step_strang_ssprk_33_cn(
631 StateVector &state_vector, Number t, Number tau_max)
632 {
633 // FIXME: avoid code duplication with step_ssprk_33
634
635#ifdef DEBUG_OUTPUT
636 std::cout << "TimeIntegrator<dim, Number>::step_strang_ssprk_33_cn()"
637 << std::endl;
638#endif
639
640 Assert(efficiency_ == 2., dealii::ExcInternalError());
641
642 /* First explicit SSPRK 3 step with final result in temp_[0]: */
643
644 Number tau = hyperbolic_module_->template step<0>(
645 state_vector, {}, {}, temp_[0], Number(0.0), tau_max / efficiency_);
646
647 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
648 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
649 sadd(temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), state_vector);
650
651 hyperbolic_module_->prepare_state_vector(temp_[1], t + 0.5 * tau);
652 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
653 sadd(temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), state_vector);
654
655 /* Implicit Crank-Nicolson step with final result in temp_[2]: */
656
657 try {
658 parabolic_module_->crank_nicolson_step(temp_[0], t, temp_[2], 2.0 * tau);
659 } catch (Restart &restart) {
660 /* Adjust suggested_tau_max. We multiply with efficiency_ again later */
661 restart.suggested_tau_max /= efficiency_;
662 throw;
663 }
664
665 /* Second SSPRK 3 step with final result in temp_[0]: */
666
667 hyperbolic_module_->prepare_state_vector( temp_[2], t + 1.0 * tau);
668 hyperbolic_module_->template step<0>( temp_[2], {}, {}, temp_[0], tau);
669
670 hyperbolic_module_->prepare_state_vector(temp_[0], t + 2.0 * tau);
671 hyperbolic_module_->template step<0>(temp_[0], {}, {}, temp_[1], tau);
672 sadd(temp_[1], Number(1.0 / 4.0), Number(3.0 / 4.0), temp_[2]);
673
674 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.5 * tau);
675 hyperbolic_module_->template step<0>(temp_[1], {}, {}, temp_[0], tau);
676 sadd(temp_[0], Number(2.0 / 3.0), Number(1.0 / 3.0), temp_[2]);
677
678 state_vector.swap(temp_[0]);
679 return efficiency_ * tau;
680 }
681
682
683 template <typename Description, int dim, typename Number>
684 Number TimeIntegrator<Description, dim, Number>::step_strang_erk_33_cn(
685 StateVector &state_vector, Number t, Number tau_max)
686 {
687 // FIXME: refactor to avoid code duplication with step_erk_33
688
689#ifdef DEBUG_OUTPUT
690 std::cout << "TimeIntegrator<dim, Number>::step_strang_erk_33_cn()"
691 << std::endl;
692#endif
693
694 Assert(efficiency_ == 6., dealii::ExcInternalError());
695
696 /* First explicit ERK(3,3,1) step with final result in temp_[2]: */
697
698 Number tau = hyperbolic_module_->template step<0>(
699 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
700
701 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
702 hyperbolic_module_->template step<1>(
703 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
704
705 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
706 hyperbolic_module_->template step<2>(temp_[1],
707 {{state_vector, temp_[0]}},
708 {{Number(0.75), Number(-2.)}},
709 temp_[2],
710 tau);
711
712 /* Implicit Crank-Nicolson step with final result in temp_[3]: */
713
714 try {
715 parabolic_module_->crank_nicolson_step(temp_[2], t, temp_[3], 6.0 * tau);
716 } catch (Restart &restart) {
717 /* Adjust suggested_tau_max. We multiply with efficiency_ again later */
718 restart.suggested_tau_max /= efficiency_;
719 throw;
720 }
721
722 /* Second explicit ERK(3,3,1) 3 step with final result in temp_[2]: */
723
724 hyperbolic_module_->prepare_state_vector(temp_[3], t + 3.0 * tau);
725 hyperbolic_module_->template step<0>(
726 temp_[3], {}, {}, temp_[0], tau);
727
728 hyperbolic_module_->prepare_state_vector(temp_[0], t + 4.0 * tau);
729 hyperbolic_module_->template step<1>(
730 temp_[0], {{ temp_[3]}}, {{Number(-1.)}}, temp_[1], tau);
731
732 hyperbolic_module_->prepare_state_vector(temp_[1], t + 5.0 * tau);
733 hyperbolic_module_->template step<2>(temp_[1],
734 {{ temp_[3], temp_[0]}},
735 {{Number(0.75), Number(-2.)}},
736 temp_[2],
737 tau);
738
739 state_vector.swap(temp_[2]);
740 return efficiency_ * tau;
741 }
742
743
744 template <typename Description, int dim, typename Number>
745 Number TimeIntegrator<Description, dim, Number>::step_strang_erk_43_cn(
746 StateVector &state_vector, Number t, Number tau_max)
747 {
748 // FIXME: refactor to avoid code duplication with step_erk_43
749
750#ifdef DEBUG_OUTPUT
751 std::cout << "TimeIntegrator<dim, Number>::step_strang_erk_43_cn()"
752 << std::endl;
753#endif
754
755 Assert(efficiency_ == 8., dealii::ExcInternalError());
756
757 /* First explicit ERK(4,3,1) step with final result in temp_[3]: */
758
759 Number tau = hyperbolic_module_->template step<0>(
760 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
761
762 hyperbolic_module_->prepare_state_vector(temp_[0], t + 1.0 * tau);
763 hyperbolic_module_->template step<1>(
764 temp_[0], {{state_vector}}, {{Number(-1.)}}, temp_[1], tau);
765
766 hyperbolic_module_->prepare_state_vector(temp_[1], t + 2.0 * tau);
767 hyperbolic_module_->template step<1>(
768 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
769
770 hyperbolic_module_->prepare_state_vector(temp_[2], t + 3.0 * tau);
771 hyperbolic_module_->template step<2>(temp_[2],
772 {{temp_[0], temp_[1]}},
773 {{Number(5. / 3.), Number(-10. / 3.)}},
774 temp_[3],
775 tau);
776
777 /* Implicit Crank-Nicolson step with final result in temp_[2]: */
778
779 try {
780 parabolic_module_->crank_nicolson_step(temp_[3], t, temp_[2], 8.0 * tau);
781 } catch (Restart &restart) {
782 /* Adjust suggested_tau_max. We multiply with efficiency_ again later */
783 restart.suggested_tau_max /= efficiency_;
784 throw;
785 }
786
787 /* Second explicit ERK(4,3,1) step with final result in temp_[3]: */
788
789 hyperbolic_module_->prepare_state_vector(temp_[2], t + 4.0 * tau);
790 hyperbolic_module_->template step<0>(
791 temp_[2], {}, {}, temp_[0], tau);
792
793 hyperbolic_module_->prepare_state_vector(temp_[0], t + 5.0 * tau);
794 hyperbolic_module_->template step<1>(
795 temp_[0], {{ temp_[2]}}, {{Number(-1.)}}, temp_[1], tau);
796
797 hyperbolic_module_->prepare_state_vector(temp_[1], t + 6.0 * tau);
798 hyperbolic_module_->template step<1>(
799 temp_[1], {{temp_[0]}}, {{Number(-1.)}}, temp_[2], tau);
800
801 hyperbolic_module_->prepare_state_vector(temp_[2], t + 7.0 * tau);
802 hyperbolic_module_->template step<2>(temp_[2],
803 {{temp_[0], temp_[1]}},
804 {{Number(5. / 3.), Number(-10. / 3.)}},
805 temp_[3],
806 tau);
807
808 state_vector.swap(temp_[3]);
809 return efficiency_ * tau;
810 }
811
812
813 template <typename Description, int dim, typename Number>
814 Number TimeIntegrator<Description, dim, Number>::step_imex_11(
815 StateVector &state_vector, Number t, Number tau_max)
816 {
817#ifdef DEBUG_OUTPUT
818 std::cout << "TimeIntegrator<dim, Number>::step_imex_11()" << std::endl;
819#endif
820
821 Assert(efficiency_ == 1., dealii::ExcInternalError());
822
823 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
824 Number tau = hyperbolic_module_->template step<0>(
825 state_vector, {}, {}, temp_[0], Number(0.), tau_max);
826
827 /* Implicit step 1: T1 <- {T0, 1} at time t -> t + tau */
828 parabolic_module_->template backward_euler_step<0>(
829 temp_[0], t, {}, {}, temp_[1], 1.0 * tau);
830
831 state_vector.swap(temp_[1]);
832 return tau;
833 }
834
835
836 template <typename Description, int dim, typename Number>
837 Number TimeIntegrator<Description, dim, Number>::step_imex_22(
838 StateVector &state_vector, Number t, Number tau_max)
839 {
840#ifdef DEBUG_OUTPUT
841 std::cout << "TimeIntegrator<dim, Number>::step_imex_22()" << std::endl;
842#endif
843
844 Assert(efficiency_ == 2., dealii::ExcInternalError());
845
846 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
847 Number tau = hyperbolic_module_->template step<0>(
848 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
849
850 /* Implicit step 1: T1 <- {T0, 1} at time t -> t + tau */
851 parabolic_module_->template backward_euler_step<0>(
852 temp_[0], t, {}, {}, temp_[1], tau);
853
854 /* Explicit step 2: T2 <- {T1, 2} and {U_old, -1} at t + tau -> t + 2 tau */
855 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.0 * tau);
856 hyperbolic_module_->template step<1>(
857 temp_[1], {{state_vector}}, {{Number(-1.)}}, temp_[2], tau);
858
859 /* Implicit step 2: T3 <- {T2, 0} and {U_old, 1} at t + tau -> t + 2 tau */
860 parabolic_module_->template backward_euler_step<1>(temp_[2],
861 t + 1.0 * tau,
862 {{state_vector}},
863 {{Number(1.)}},
864 temp_[3],
865 tau);
866
867 state_vector.swap(temp_[3]);
868 return efficiency_ * tau;
869 }
870
871
872 template <typename Description, int dim, typename Number>
873 Number TimeIntegrator<Description, dim, Number>::step_imex_33(
874 StateVector &state_vector, Number t, Number tau_max)
875 {
876#ifdef DEBUG_OUTPUT
877 std::cout << "TimeIntegrator<dim, Number>::step_imex_33()" << std::endl;
878#endif
879
880 Assert(efficiency_ == 3., dealii::ExcInternalError());
881
882 /* IMEX(3, 3; 1), see @cite ErnGuermond2023, Sec. 4.3. */
883
884 const Number gamma = Number(0.5) + std::sqrt(Number(3.0)) / Number(6.0);
885
886 /* Explicit step 1: T0 <- {U_old, 1} at time t -> t + tau */
887 Number tau = hyperbolic_module_->template step<0>(
888 state_vector, {}, {}, temp_[0], Number(0.), tau_max / efficiency_);
889
890 /* Implicit step 1: T1 <- {U_old, 1 - 3*gamma} at time t -> t + tau */
891 parabolic_module_->template backward_euler_step<1>(
892 temp_[0],
893 t,
894 {{state_vector}},
895 {{Number(1. - 3. * gamma)}},
896 temp_[1],
897 tau);
898
899 /* Explicit step 2: T2 <- {U_old, -1} and {T1, 2} at time t -> t + 2 tau */
900 hyperbolic_module_->prepare_state_vector(temp_[1], t + 1.0 * tau);
901 hyperbolic_module_->template step<1>(
902 temp_[1], {{state_vector}}, {{Number(-1.)}}, temp_[2], tau);
903
904 /*
905 * Implicit step 2:
906 * T3 <- {U_old, 6*gamma-1} and {T1, 2-9*gamma} at t -> t + * 2 tau
907 */
908 parabolic_module_->template backward_euler_step<2>(
909 temp_[2],
910 t + tau,
911 {{state_vector, temp_[1]}},
912 {{Number(6. * gamma - 1.), Number(2. - 9 * gamma)}},
913 temp_[3],
914 tau);
915
916 /* Explicit step 3: T4 <- {U_old, 3 / 4} and {T1, -2} at t -> t + 3 tau */
917 hyperbolic_module_->prepare_state_vector(temp_[3], t + 2. * tau);
918 hyperbolic_module_->template step<2>(temp_[3],
919 {{state_vector, temp_[1]}},
920 {{Number(0.75), Number(-2.)}},
921 temp_[4],
922 tau);
923
924 /* Implicit step 3: */
925 parabolic_module_->template backward_euler_step<3>(
926 temp_[4],
927 t + 2. * tau,
928 {{state_vector, temp_[1], temp_[3]}},
929 {{Number(0.75 - 3. * gamma),
930 Number(6. * gamma - 2.),
931 Number(9. / 4. - 3. * gamma)}},
932 temp_[5],
933 tau);
934
935 state_vector.swap(temp_[5]);
936 return efficiency_ * tau;
937 }
938
939} /* namespace ryujin */
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:50
ryujin::TimeIntegrator::StateVector
typename View::StateVector StateVector
Definition time_integrator.h:266
ryujin::TimeIntegrator::prepare_state_vector
void prepare_state_vector(StateVector &state_vector, Number t) const
Definition time_integrator.template.h:252
ryujin::TimeIntegrator::step
Number step(StateVector &state_vector, Number t, Number t_final=std::numeric_limits< Number >::max())
Definition time_integrator.template.h:269
ryujin::selected_memory_space_t
std::conditional_t< have_separate_memory_spaces, dealii::MemorySpace::Default, dealii::MemorySpace::Host > selected_memory_space_t
Definition gpu.h:65
ryujin::IDViolationStrategy::warn
@ warn
ryujin::IDViolationStrategy::raise_exception
@ raise_exception
ryujin::ParabolicModule
typename Description::template ParabolicModule< dim, Number > ParabolicModule
Definition parabolic_module.h:19
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::CFLRecoveryStrategy::cruise_control
@ cruise_control
ryujin::sadd
void sadd(StateVector &dst, const Number s, const Number b, const StateVector &src)
Definition time_integrator.template.h:18
time_integrator.h
Generated by doxygen 1.9.8