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