ExaDG
Loading...
Searching...
No Matches
explicit_runge_kutta.h
1/* ______________________________________________________________________
2 *
3 * ExaDG - High-Order Discontinuous Galerkin for the Exa-Scale
4 *
5 * Copyright (C) 2021 by the ExaDG authors
6 *
7 * This program is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program. If not, see <https://www.gnu.org/licenses/>.
19 * ______________________________________________________________________
20 */
21
22#ifndef INCLUDE_CONVECTION_DIFFUSION_EXPLICIT_RUNGE_KUTTA_H_
23#define INCLUDE_CONVECTION_DIFFUSION_EXPLICIT_RUNGE_KUTTA_H_
24
25namespace ExaDG
26{
27template<typename Operator, typename VectorType>
28class ExplicitTimeIntegrator
29{
30public:
31 ExplicitTimeIntegrator(std::shared_ptr<Operator> operator_in) : underlying_operator(operator_in)
32 {
33 }
34
35 virtual ~ExplicitTimeIntegrator()
36 {
37 }
38
39 virtual void
40 solve_timestep(VectorType & dst, VectorType & src, double const time, double const time_step) = 0;
41
42 virtual unsigned int
43 get_order() const = 0;
44
45protected:
46 std::shared_ptr<Operator> underlying_operator;
47};
48
49/*
50 * Classical, explicit Runge-Kutta time integration schemes of order 1-4
51 */
52template<typename Operator, typename VectorType>
53class ExplicitRungeKuttaTimeIntegrator : public ExplicitTimeIntegrator<Operator, VectorType>
54{
55public:
56 // Constructor
57 ExplicitRungeKuttaTimeIntegrator(unsigned int order_time_integrator,
58 std::shared_ptr<Operator> const operator_in)
59 : ExplicitTimeIntegrator<Operator, VectorType>(operator_in), order(order_time_integrator)
60 {
61 // initialize vectors
62 if(order >= 2)
63 this->underlying_operator->initialize_dof_vector(vec_rhs);
64 if(order >= 3)
65 this->underlying_operator->initialize_dof_vector(vec_temp);
66 }
67
68 void
69 solve_timestep(VectorType & dst,
70 VectorType & src,
71 double const time,
72 double const time_step) final
73 {
74 if(order == 1) // explicit Euler method
75 {
76 /*
77 * Butcher table
78 *
79 * 0 |
80 * --------
81 * | 1
82 *
83 */
84
85 this->underlying_operator->evaluate(dst, src, time);
86 dst.sadd(time_step, 1., src);
87 }
88 else if(order == 2) // Runge-Kutta method of order 2
89 {
90 /*
91 * Butcher table
92 *
93 * 0 |
94 * 1/2 | 1/2
95 * ---------------
96 * | 0 1
97 *
98 */
99
100 // stage 1
101 this->underlying_operator->evaluate(vec_rhs, src, time);
102
103 // stage 2
104 vec_rhs.sadd(time_step / 2., 1., src);
105 this->underlying_operator->evaluate(dst, vec_rhs, time + time_step / 2.);
106 dst.sadd(time_step, 1., src);
107 }
108 else if(order == 3) // Heun's method of order 3
109 {
110 /*
111 * Butcher table
112 *
113 * 0 |
114 * 1/3 | 1/3
115 * 2/3 | 0 2/3
116 * -------------------
117 * | 1/4 0 3/4
118 *
119 */
120
121 dst = src;
122
123 // stage 1
124 this->underlying_operator->evaluate(vec_temp, src, time);
125 dst.add(1. * time_step / 4., vec_temp);
126
127 // stage 2
128 vec_rhs.equ(1., src);
129 vec_rhs.add(time_step / 3., vec_temp);
130 this->underlying_operator->evaluate(vec_temp, vec_rhs, time + time_step / 3.);
131
132 // stage 3
133 vec_rhs.equ(1., src);
134 vec_rhs.add(2. * time_step / 3., vec_temp);
135 this->underlying_operator->evaluate(vec_temp, vec_rhs, time + 2. * time_step / 3.);
136 dst.add(3. * time_step / 4., vec_temp);
137 }
138 else if(order == 4) // classical 4th order Runge-Kutta method
139 {
140 /*
141 * Butcher table
142 *
143 * 0 |
144 * 1/2 | 1/2
145 * 1/2 | 0 1/2
146 * 1 | 0 0 1
147 * -----------------------
148 * | 1/6 2/6 2/6 1/6
149 *
150 */
151
152 dst = src;
153
154 // stage 1
155 this->underlying_operator->evaluate(vec_temp, src, time);
156 dst.add(time_step / 6., vec_temp);
157
158 // stage 2
159 vec_rhs.equ(1., src);
160 vec_rhs.add(time_step / 2., vec_temp);
161 this->underlying_operator->evaluate(vec_temp, vec_rhs, time + time_step / 2.);
162 dst.add(time_step / 3., vec_temp);
163
164 // stage 3
165 vec_rhs.equ(1., src);
166 vec_rhs.add(time_step / 2., vec_temp);
167 this->underlying_operator->evaluate(vec_temp, vec_rhs, time + time_step / 2.);
168 dst.add(time_step / 3., vec_temp);
169
170 // stage 4
171 vec_rhs.equ(1., src);
172 vec_rhs.add(time_step, vec_temp);
173 this->underlying_operator->evaluate(vec_temp, vec_rhs, time + time_step);
174 dst.add(time_step / 6., vec_temp);
175 }
176 else
177 {
178 AssertThrow(order <= 4,
179 dealii::ExcMessage(
180 "Explicit Runge-Kutta method only implemented for order <= 4!"));
181 }
182 }
183
184 unsigned int
185 get_order() const final
186 {
187 return order;
188 }
189
190private:
191 unsigned int order;
192
193 VectorType vec_rhs, vec_temp;
194};
195
196
197
198/****************************************************************************************
199 * *
200 * Low-storage, explicit Runge-Kutta methods according to *
201 * *
202 * Kennedy et al. (2000, "Low-storage, explicit Runge-Kutta schemes for the *
203 * compressible Navier-Stokes equations") *
204 * *
205 * classification of methods: *
206 * *
207 * RKq(p)s[rR]X where ... *
208 * *
209 * ... q is the order of the main scheme *
210 * ... p is the order of the embedded scheme *
211 * ... s is the number of stages *
212 * ... r is the number of registers used *
213 * ... X denotes the optimization strategy/criterion, e.g, *
214 * *
215 * X = C: linear-stability-error compromise *
216 * X = S: maximum linear stability *
217 * X = F: FSAL - first same as last *
218 * X = M: minimum truncation error *
219 * X = N: maximum nonlinear stability *
220 * X = P: minimum phase error *
221 * *
222 * For applications where the time step is selected according to stability *
223 * considerations (time step size close to CFL limit, temporal error small enough), *
224 * schemes of type C or type S should be used. *
225 * *
226 ****************************************************************************************/
227
228/*
229 * Low storage Runge-Kutta method of order 3 with 4 stages and 2 registers according to
230 * Kennedy et al. (2000), where this method is denoted as RK3(2)4[2R+]C,
231 * see Table 1 on page 189 for the coefficients.
232 */
233template<typename Operator, typename VectorType>
234class LowStorageRK3Stage4Reg2C : public ExplicitTimeIntegrator<Operator, VectorType>
235{
236public:
237 LowStorageRK3Stage4Reg2C(std::shared_ptr<Operator> const operator_in)
238 : ExplicitTimeIntegrator<Operator, VectorType>(operator_in)
239 {
240 }
241
242 void
243 solve_timestep(VectorType & vec_np,
244 VectorType & vec_n,
245 double const time,
246 double const time_step) final
247 {
248 if(not vec_tmp1.partitioners_are_globally_compatible(*vec_n.get_partitioner()))
249 {
250 vec_tmp1.reinit(vec_np);
251 }
252
253 /*
254 * Butcher table
255 *
256 * c1 |
257 * c2 | a21
258 * c3 | b1 a32
259 * c4 | b1 b2 a43
260 * ---------------------
261 * | b1 b2 b3 b4
262 *
263 */
264
265 double const a21 = 11847461282814. / 36547543011857.;
266 double const a32 = 3943225443063. / 7078155732230.;
267 double const a43 = -346793006927. / 4029903576067.;
268
269 double const b1 = 1017324711453. / 9774461848756.;
270 double const b2 = 8237718856693. / 13685301971492.;
271 double const b3 = 57731312506979. / 19404895981398.;
272 double const b4 = -101169746363290. / 37734290219643.;
273
274 double const c1 = 0.;
275 double const c2 = a21;
276 double const c3 = b1 + a32;
277 double const c4 = b1 + b2 + a43;
278
279 // stage 1
280 this->underlying_operator->evaluate(vec_tmp1, vec_n /* u_1 */, time + c1 * time_step);
281 vec_n.add(a21 * time_step, vec_tmp1); /* = u_2 */
282 vec_np = vec_n;
283 vec_np.add((b1 - a21) * time_step, vec_tmp1); /* = u_p */
284
285 // stage 2
286 this->underlying_operator->evaluate(vec_tmp1, vec_n /* u_2 */, time + c2 * time_step);
287 vec_np.add(a32 * time_step, vec_tmp1); /* = u_3 */
288 vec_n = vec_np;
289 vec_n.add((b2 - a32) * time_step, vec_tmp1); /* = u_p */
290
291 // stage 3
292 this->underlying_operator->evaluate(vec_tmp1, vec_np /* u_3 */, time + c3 * time_step);
293 vec_n.add(a43 * time_step, vec_tmp1); /* = u_4 */
294 vec_np = vec_n;
295 vec_np.add((b3 - a43) * time_step, vec_tmp1); /* = u_p */
296
297 // stage 4
298 this->underlying_operator->evaluate(vec_tmp1, vec_n /* u_3 */, time + c4 * time_step);
299 vec_np.add(b4 * time_step, vec_tmp1); /* = u_p */
300 }
301
302 unsigned int
303 get_order() const final
304 {
305 return 3;
306 }
307
308private:
309 VectorType vec_tmp1;
310};
311
312
313/*
314 * Low storage Runge-Kutta method of order 4 with 5 stages and 2 registers according to
315 * Kennedy et al. (2000), where this method is denoted as RK4(3)5[2R+]C,
316 * see Table 1 on page 189 for the coefficients.
317 */
318template<typename Operator, typename VectorType>
319class LowStorageRK4Stage5Reg2C : public ExplicitTimeIntegrator<Operator, VectorType>
320{
321public:
322 LowStorageRK4Stage5Reg2C(std::shared_ptr<Operator> const operator_in)
323 : ExplicitTimeIntegrator<Operator, VectorType>(operator_in)
324 {
325 }
326
327 void
328 solve_timestep(VectorType & vec_np,
329 VectorType & vec_n,
330 double const time,
331 double const time_step) final
332 {
333 if(not vec_tmp1.partitioners_are_globally_compatible(*vec_n.get_partitioner()))
334 {
335 vec_tmp1.reinit(vec_np);
336 }
337
338 /*
339 * Butcher table
340 *
341 * c1 |
342 * c2 | a21
343 * c3 | b1 a32
344 * c4 | b1 b2 a43
345 * c5 | b1 b2 b3 a54
346 * ------------------------
347 * | b1 b2 b3 b4 b5
348 *
349 */
350
351 double const a21 = 970286171893. / 4311952581923.;
352 double const a32 = 6584761158862. / 12103376702013.;
353 double const a43 = 2251764453980. / 15575788980749.;
354 double const a54 = 26877169314380. / 34165994151039.;
355
356 double const b1 = 1153189308089. / 22510343858157.;
357 double const b2 = 1772645290293. / 4653164025191.;
358 double const b3 = -1672844663538. / 4480602732383.;
359 double const b4 = 2114624349019. / 3568978502595.;
360 double const b5 = 5198255086312. / 14908931495163.;
361
362 double const c1 = 0.;
363 double const c2 = a21;
364 double const c3 = b1 + a32;
365 double const c4 = b1 + b2 + a43;
366 double const c5 = b1 + b2 + b3 + a54;
367
368 // stage 1
369 this->underlying_operator->evaluate(vec_tmp1, vec_n /* u_1 */, time + c1 * time_step);
370 vec_n.add(a21 * time_step, vec_tmp1); /* = u_2 */
371 vec_np = vec_n;
372 vec_np.add((b1 - a21) * time_step, vec_tmp1); /* = u_p */
373
374 // stage 2
375 this->underlying_operator->evaluate(vec_tmp1, vec_n /* u_2 */, time + c2 * time_step);
376 vec_np.add(a32 * time_step, vec_tmp1); /* = u_3 */
377 vec_n = vec_np;
378 vec_n.add((b2 - a32) * time_step, vec_tmp1); /* = u_p */
379
380 // stage 3
381 this->underlying_operator->evaluate(vec_tmp1, vec_np /* u_3 */, time + c3 * time_step);
382 vec_n.add(a43 * time_step, vec_tmp1); /* = u_4 */
383 vec_np = vec_n;
384 vec_np.add((b3 - a43) * time_step, vec_tmp1); /* = u_p */
385
386 // stage 4
387 this->underlying_operator->evaluate(vec_tmp1, vec_n /* u_3 */, time + c4 * time_step);
388 vec_np.add(a54 * time_step, vec_tmp1); /* = u_5 */
389 vec_n = vec_np;
390 vec_n.add((b4 - a54) * time_step, vec_tmp1); /* = u_p */
391
392 // stage 5
393 this->underlying_operator->evaluate(vec_tmp1, vec_np /* u_4 */, time + c5 * time_step);
394 vec_np = vec_n;
395 vec_np.add(b5 * time_step, vec_tmp1);
396 }
397
398 unsigned int
399 get_order() const final
400 {
401 return 4;
402 }
403
404private:
405 VectorType vec_tmp1;
406};
407
408/*
409 * Low storage Runge-Kutta method of order 4 with 5 stages and 3 registers according to
410 * Kennedy et al. (2000), where this method is denoted as RK4(3)5[3R+]C,
411 * see Table 2 on page 190 for the coefficients.
412 */
413template<typename Operator, typename VectorType>
414class LowStorageRK4Stage5Reg3C : public ExplicitTimeIntegrator<Operator, VectorType>
415{
416public:
417 LowStorageRK4Stage5Reg3C(std::shared_ptr<Operator> const operator_in)
418 : ExplicitTimeIntegrator<Operator, VectorType>(operator_in)
419 {
420 }
421
422 void
423 solve_timestep(VectorType & vec_np,
424 VectorType & vec_n,
425 double const time,
426 double const time_step) final
427 {
428 if(not vec_tmp1.partitioners_are_globally_compatible(*vec_n.get_partitioner()))
429 {
430 vec_tmp1.reinit(vec_np, true);
431 vec_tmp2.reinit(vec_np, true);
432 }
433
434 /*
435 * Butcher table
436 *
437 * c1 |
438 * c2 | a21
439 * c3 | a31 a32
440 * c4 | b1 a42 a43
441 * c5 | b1 b2 a53 a54
442 * ------------------------
443 * | b1 b2 b3 b4 b5
444 *
445 */
446 double const a21 = 2365592473904. / 8146167614645.;
447 double const a32 = 4278267785271. / 6823155464066.;
448 double const a43 = 2789585899612. / 8986505720531.;
449 double const a54 = 15310836689591. / 24358012670437.;
450
451 double const a31 = -722262345248. / 10870640012513.;
452 double const a42 = 1365858020701. / 8494387045469.;
453 double const a53 = 3819021186. / 2763618202291.;
454
455 double const b1 = 846876320697. / 6523801458457.;
456 double const b2 = 3032295699695. / 12397907741132.;
457 double const b3 = 612618101729. / 6534652265123.;
458 double const b4 = 1155491934595. / 2954287928812.;
459 double const b5 = 707644755468. / 5028292464395.;
460
461 double const c1 = 0.;
462 double const c2 = a21;
463 double const c3 = a31 + a32;
464 double const c4 = b1 + a42 + a43;
465 double const c5 = b1 + b2 + a53 + a54;
466
467 // stage 1
468 this->underlying_operator->evaluate(vec_np /* F_1 */, vec_n /* u_1 */, time + c1 * time_step);
469 vec_n.add(a21 * time_step, vec_np /* F_1 */); /* = u_2 */
470 vec_tmp1 = vec_n /* u_2 */;
471 vec_tmp1.add((b1 - a21) * time_step, vec_np /* F_1 */); /* = u_p */
472
473 // stage 2
474 this->underlying_operator->evaluate(vec_tmp2 /* F_2 */, vec_n /* u_2 */, time + c2 * time_step);
475 vec_tmp1.add(a32 * time_step,
476 vec_tmp2 /* F_2 */,
477 (a31 - b1) * time_step,
478 vec_np /* F_1 */); /* = u_3 */
479 vec_np.sadd((b1 - a31) * time_step, 1, vec_tmp1 /* u_3 */);
480 vec_np.add((b2 - a32) * time_step, vec_tmp2 /* F_2 */); /* u_p */
481
482 // stage 3
483 this->underlying_operator->evaluate(vec_n /* F_3 */, vec_tmp1 /* u_3 */, time + c3 * time_step);
484 vec_np.add(a43 * time_step,
485 vec_n /* F_3 */,
486 (a42 - b2) * time_step,
487 vec_tmp2 /* F_2 */); /* = u_4 */
488 vec_tmp2.sadd((b2 - a42) * time_step, 1, vec_np /* u_4 */);
489 vec_tmp2.add((b3 - a43) * time_step, vec_n /* F_3 */); /* u_p */
490
491 // stage 4
492 this->underlying_operator->evaluate(vec_tmp1 /* F_4 */,
493 vec_np /* u_4 */,
494 time + c4 * time_step);
495 vec_tmp2.add(a54 * time_step,
496 vec_tmp1 /* F_4 */,
497 (a53 - b3) * time_step,
498 vec_n /* F_3 */); /* = u_5 */
499 vec_n.sadd((b3 - a53) * time_step, 1, vec_tmp2 /* u_5 */);
500 vec_n.add((b4 - a54) * time_step, vec_tmp1 /* F_4 */); /* = u_p */
501
502 // stage 5
503 this->underlying_operator->evaluate(vec_tmp1 /* F_5 */,
504 vec_tmp2 /* u_5 */,
505 time + c5 * time_step);
506 vec_np = vec_n;
507 vec_np.add(b5 * time_step, vec_tmp1 /* F_5 */); /* = u_p */
508 }
509
510 unsigned int
511 get_order() const final
512 {
513 return 4;
514 }
515
516private:
517 VectorType vec_tmp1, vec_tmp2;
518};
519
520/*
521 * Low storage Runge-Kutta method of order 5 with 9 stages and 2 registers according to
522 * Kennedy et al. (2000), where this method is denoted as RK5(4)9[2R+]S,
523 * see Table 1 on page 189 for the coefficients.
524 */
525template<typename Operator, typename VectorType>
526class LowStorageRK5Stage9Reg2S : public ExplicitTimeIntegrator<Operator, VectorType>
527{
528public:
529 LowStorageRK5Stage9Reg2S(std::shared_ptr<Operator> const operator_in)
530 : ExplicitTimeIntegrator<Operator, VectorType>(operator_in)
531 {
532 }
533
534 void
535 solve_timestep(VectorType & vec_np,
536 VectorType & vec_n,
537 double const time,
538 double const time_step) final
539 {
540 if(not vec_tmp1.partitioners_are_globally_compatible(*vec_n.get_partitioner()))
541 {
542 vec_tmp1.reinit(vec_np);
543 }
544
545 /*
546 * Butcher table
547 *
548 * c1 |
549 * c2 | a21
550 * c3 | b1 a32
551 * c4 | b1 b2 a43
552 * c5 | b1 b2 b3 a54
553 * c6 | b1 b2 b3 b4 a65
554 * c7 | b1 b2 b3 b4 b5 a76
555 * c8 | b1 b2 b3 b4 b5 b6 a87
556 * c9 | b1 b2 b3 b4 b5 b6 b7 a98
557 * -----------------------------------------
558 * | b1 b2 b3 b4 b5 b6 b7 b8 b9
559 *
560 */
561
562 double const a21 = 1107026461565. / 5417078080134.;
563 double const a32 = 38141181049399. / 41724347789894.;
564 double const a43 = 493273079041. / 11940823631197.;
565 double const a54 = 1851571280403. / 6147804934346.;
566 double const a65 = 11782306865191. / 62590030070788.;
567 double const a76 = 9452544825720. / 13648368537481.;
568 double const a87 = 4435885630781. / 26285702406235.;
569 double const a98 = 2357909744247. / 11371140753790.;
570
571 double const b1 = 2274579626619. / 23610510767302.;
572 double const b2 = 693987741272. / 12394497460941.;
573 double const b3 = -347131529483. / 15096185902911.;
574 double const b4 = 1144057200723. / 32081666971178.;
575 double const b5 = 1562491064753. / 11797114684756.;
576 double const b6 = 13113619727965. / 44346030145118.;
577 double const b7 = 393957816125. / 7825732611452.;
578 double const b8 = 720647959663. / 6565743875477.;
579 double const b9 = 3559252274877. / 14424734981077.;
580
581 double const c1 = 0.;
582 double const c2 = a21;
583 double const c3 = b1 + a32;
584 double const c4 = b1 + b2 + a43;
585 double const c5 = b1 + b2 + b3 + a54;
586 double const c6 = b1 + b2 + b3 + b4 + a65;
587 double const c7 = b1 + b2 + b3 + b4 + b5 + a76;
588 double const c8 = b1 + b2 + b3 + b4 + b5 + b6 + a87;
589 double const c9 = b1 + b2 + b3 + b4 + b5 + b6 + b7 + a98;
590
591 // stage 1
592 this->underlying_operator->evaluate(vec_tmp1, vec_n /* u_1 */, time + c1 * time_step);
593 vec_n.add(a21 * time_step, vec_tmp1); /* = u_2 */
594 vec_np = vec_n;
595 vec_np.add((b1 - a21) * time_step, vec_tmp1); /* = u_p */
596
597 // stage 2
598 this->underlying_operator->evaluate(vec_tmp1, vec_n /* u_2 */, time + c2 * time_step);
599 vec_np.add(a32 * time_step, vec_tmp1); /* = u_3 */
600 vec_n = vec_np;
601 vec_n.add((b2 - a32) * time_step, vec_tmp1); /* = u_p */
602
603 // stage 3
604 this->underlying_operator->evaluate(vec_tmp1, vec_np /* u_3 */, time + c3 * time_step);
605 vec_n.add(a43 * time_step, vec_tmp1); /* = u_4 */
606 vec_np = vec_n;
607 vec_np.add((b3 - a43) * time_step, vec_tmp1); /* = u_p */
608
609 // stage 4
610 this->underlying_operator->evaluate(vec_tmp1, vec_n /* u_4 */, time + c4 * time_step);
611 vec_np.add(a54 * time_step, vec_tmp1); /* = u_5 */
612 vec_n = vec_np;
613 vec_n.add((b4 - a54) * time_step, vec_tmp1); /* = u_p */
614
615 // stage 5
616 this->underlying_operator->evaluate(vec_tmp1, vec_np /* u_5 */, time + c5 * time_step);
617 vec_n.add(a65 * time_step, vec_tmp1); /* = u_6 */
618 vec_np = vec_n;
619 vec_np.add((b5 - a65) * time_step, vec_tmp1); /* = u_p */
620
621 // stage 6
622 this->underlying_operator->evaluate(vec_tmp1, vec_n /* u_6 */, time + c6 * time_step);
623 vec_np.add(a76 * time_step, vec_tmp1); /* = u_7 */
624 vec_n = vec_np;
625 vec_n.add((b6 - a76) * time_step, vec_tmp1); /* = u_p */
626
627 // stage 7
628 this->underlying_operator->evaluate(vec_tmp1, vec_np /* u_7 */, time + c7 * time_step);
629 vec_n.add(a87 * time_step, vec_tmp1); /* = u_8 */
630 vec_np = vec_n;
631 vec_np.add((b7 - a87) * time_step, vec_tmp1); /* = u_p */
632
633 // stage 8
634 this->underlying_operator->evaluate(vec_tmp1, vec_n /* u_8 */, time + c8 * time_step);
635 vec_np.add(a98 * time_step, vec_tmp1); /* = u_9 */
636 vec_n = vec_np;
637 vec_n.add((b8 - a98) * time_step, vec_tmp1); /* = u_p */
638
639 // stage 9
640 this->underlying_operator->evaluate(vec_tmp1, vec_np /* u_9 */, time + c9 * time_step);
641 vec_np = vec_n;
642 vec_np.add(b9 * time_step, vec_tmp1);
643 }
644
645 unsigned int
646 get_order() const final
647 {
648 return 5;
649 }
650
651private:
652 VectorType vec_tmp1;
653};
654
655
656/*
657 * Explicit Runge-Kutta of Toulorge & Desmet (2011) in low-storage format (2N scheme)
658 * of order q with additional stages s>q in order to optimize the stability region of
659 * the time integration scheme and to optimize costs = stages / CFL_crit:
660 *
661 * We consider two methods of order 3 and 4 with 7 and 8 stages, respectively.
662 * The Runge-Kutta coefficients for these schemes can be found in Toulorge & Desmet
663 * in Tables A.15 for the RKC73 scheme and Table A.21 for the RKC84 scheme.
664 * In our nomenclature, these time integration schemes are deonted as
665 * ExplRK3Stage7Reg2 and ExplRK4Stage8Reg2, respectively.
666 */
667template<typename Operator, typename VectorType>
668class LowStorageRKTD : public ExplicitTimeIntegrator<Operator, VectorType>
669{
670public:
671 LowStorageRKTD(std::shared_ptr<Operator> const operator_in,
672 unsigned int const order_in,
673 unsigned int const stages_in)
674 : ExplicitTimeIntegrator<Operator, VectorType>(operator_in), order(order_in), stages(stages_in)
675 {
676 }
677
678 void
679 solve_timestep(VectorType & vec_np,
680 VectorType & vec_n,
681 double const time,
682 double const time_step) final
683 {
684 if(not vec_tmp.partitioners_are_globally_compatible(*vec_n.get_partitioner()))
685 {
686 vec_tmp.reinit(vec_np);
687 }
688
689 std::vector<double> A;
690 A.resize(stages);
691 std::vector<double> B;
692 B.resize(stages);
693 std::vector<double> c;
694 c.resize(stages);
695
696 if(order == 3 and stages == 7)
697 {
698 A[0] = 0.;
699 A[1] = -0.8083163874983830;
700 A[2] = -1.503407858773331;
701 A[3] = -1.053064525050744;
702 A[4] = -1.463149119280508;
703 A[5] = -0.6592881281087830;
704 A[6] = -1.667891931891068;
705
706 B[0] = 0.01197052673097840;
707 B[1] = 0.8886897793820711;
708 B[2] = 0.4578382089261419;
709 B[3] = 0.5790045253338471;
710 B[4] = 0.3160214638138484;
711 B[5] = 0.2483525368264122;
712 B[6] = 0.06771230959408840;
713
714 c[0] = 0.;
715 c[1] = 0.01197052673097840;
716 c[2] = 0.1823177940361990;
717 c[3] = 0.5082168062551849;
718 c[4] = 0.6532031220148590;
719 c[5] = 0.8534401385678250;
720 c[6] = 0.9980466084623790;
721 }
722 else if(order == 4 and stages == 8)
723 {
724 A[0] = 0.;
725 A[1] = -0.7212962482279240;
726 A[2] = -0.01077336571612980;
727 A[3] = -0.5162584698930970;
728 A[4] = -1.730100286632201;
729 A[5] = -5.200129304403076;
730 A[6] = 0.7837058945416420;
731 A[7] = -0.5445836094332190;
732
733 B[0] = 0.2165936736758085;
734 B[1] = 0.1773950826411583;
735 B[2] = 0.01802538611623290;
736 B[3] = 0.08473476372541490;
737 B[4] = 0.8129106974622483;
738 B[5] = 1.903416030422760;
739 B[6] = 0.1314841743399048;
740 B[7] = 0.2082583170674149;
741
742 c[0] = 0.;
743 c[1] = 0.2165936736758085;
744 c[2] = 0.2660343487538170;
745 c[3] = 0.2840056122522720;
746 c[4] = 0.3251266843788570;
747 c[5] = 0.4555149599187530;
748 c[6] = 0.7713219317101170;
749 c[7] = 0.9199028964538660;
750 }
751 else
752 {
753 AssertThrow(false, dealii::ExcMessage("Not implemented."));
754 }
755
756 // loop over all stages
757 for(unsigned int s = 0; s < stages; s++)
758 {
759 this->underlying_operator->evaluate(vec_np, vec_n, time + c[s] * time_step);
760 vec_tmp.sadd(A[s], time_step, vec_np);
761 vec_n.add(B[s], vec_tmp);
762 }
763
764 vec_np = vec_n;
765 }
766
767 unsigned int
768 get_order() const final
769 {
770 return order;
771 }
772
773private:
774 VectorType vec_tmp;
775 unsigned int order;
776 unsigned int stages;
777};
778
779} // namespace ExaDG
780
781#endif /* INCLUDE_CONVECTION_DIFFUSION_EXPLICIT_RUNGE_KUTTA_H_ */
ExaDG
Definition driver.cpp:33
Generated by doxygen 1.13.2