ExaDG
Loading...
Searching...
No Matches
ssp_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_EXADG_TIME_INTEGRATION_SSP_RUNGE_KUTTA_H_
23#define INCLUDE_EXADG_TIME_INTEGRATION_SSP_RUNGE_KUTTA_H_
24
25// deal.II
26#include <deal.II/lac/full_matrix.h>
27
28// ExaDG
29#include <exadg/time_integration/explicit_runge_kutta.h>
30
31namespace ExaDG
32{
33/*
34 * Strong-Stability-Preserving Runge-Kutta Methods according to
35 *
36 * Kubatko et al., Optimal Strong-Stability-Preserving Runge-Kutta Time Discretizations
37 * for Discontinuous Galerkin Methods, J Sci Comput (2014) 60:313-344.
38 *
39 * The Runge-Kutta scheme is implemented in Shu-Osher form instead of the Butcher form.
40 *
41 */
42template<typename Operator, typename VectorType>
43class SSPRK : public ExplicitTimeIntegrator<Operator, VectorType>
44{
45public:
46 SSPRK(std::shared_ptr<Operator> const operator_in,
47 unsigned int const order_in,
48 unsigned int const stages_in)
49 : ExplicitTimeIntegrator<Operator, VectorType>(operator_in), order(order_in)
50 {
51 initialize_coeffs(stages_in);
52 }
53
54 void
55 solve_timestep(VectorType & vec_np,
56 VectorType & vec_n,
57 double const time,
58 double const time_step) final;
59
60 unsigned int
61 get_order() const final
62 {
63 return order;
64 }
65
66private:
67 void
68 initialize_coeffs(unsigned int const stages);
69
70 dealii::FullMatrix<double> A, B;
71 std::vector<double> c;
72 unsigned int const order;
73
74 std::vector<VectorType> u_vec, F_vec;
75};
76
77template<typename Operator, typename VectorType>
78void
79SSPRK<Operator, VectorType>::solve_timestep(VectorType & vec_np,
80 VectorType & vec_n,
81 double const time,
82 double const time_step)
83{
84 unsigned int const stages = A.m();
85
86 // Initialize vectors if necessary
87 if(u_vec.empty() or not(u_vec[0].partitioners_are_globally_compatible(*vec_n.get_partitioner())))
88 {
89 u_vec.resize(stages + 1);
90 for(unsigned int d = 0; d < u_vec.size(); ++d)
91 u_vec[d].reinit(vec_n);
92
93 F_vec.resize(stages);
94 for(unsigned int d = 0; d < F_vec.size(); ++d)
95 F_vec[d].reinit(vec_n, true);
96 }
97
98 u_vec[0] = vec_n;
99
100 for(unsigned int s = 1; s <= stages; ++s)
101 {
102 this->underlying_operator->evaluate(F_vec[s - 1], u_vec[s - 1], time + c[s - 1] * time_step);
103
104 u_vec[s] = 0.; // Do not forget to reset u_vec[s]!
105 for(unsigned int l = 0; l < s; ++l)
106 {
107 u_vec[s].add(A[s - 1][l], u_vec[l], B[s - 1][l] * time_step, F_vec[l]);
108 }
109 }
110
111 vec_np = u_vec[stages];
112}
113
114template<typename Operator, typename VectorType>
115void
116SSPRK<Operator, VectorType>::initialize_coeffs(unsigned int const stages)
117{
118 A.reinit(stages, stages);
119 B.reinit(stages, stages);
120
121 bool coefficients_are_initialized = false;
122
123 if(order == 3)
124 {
125 if(stages == 4)
126 {
127 A[0][0] = 1.;
128 A[0][1] = 0.;
129 A[0][2] = 0.;
130 A[0][3] = 0.;
131 A[1][0] = 0.522361915162541;
132 A[1][1] = 0.477638084837459;
133 A[1][2] = 0.;
134 A[1][3] = 0.;
135 A[2][0] = 0.368530939472566;
136 A[2][1] = 0.;
137 A[2][2] = 0.631469060527434;
138 A[2][3] = 0.;
139 A[3][0] = 0.334082932462285;
140 A[3][1] = 0.006966183666289;
141 A[3][2] = 0.;
142 A[3][3] = 0.658950883871426;
143
144 B[0][0] = 0.594057152884440;
145 B[0][1] = 0.;
146 B[0][2] = 0.;
147 B[0][3] = 0.;
148 B[1][0] = 0.;
149 B[1][1] = 0.283744320787718;
150 B[1][2] = 0.;
151 B[1][3] = 0.;
152 B[2][0] = 0.000000038023030;
153 B[2][1] = 0.;
154 B[2][2] = 0.375128712231540;
155 B[2][3] = 0.;
156 B[3][0] = 0.116941419604231;
157 B[3][1] = 0.004138311235266;
158 B[3][2] = 0.;
159 B[3][3] = 0.391454485963345;
160
161 coefficients_are_initialized = true;
162 }
163 else if(stages == 5)
164 {
165 A[0][0] = 1.;
166 A[0][1] = 0.;
167 A[0][2] = 0.;
168 A[0][3] = 0.;
169 A[0][4] = 0.;
170 A[1][0] = 0.495124140877703;
171 A[1][1] = 0.504875859122297;
172 A[1][2] = 0.;
173 A[1][3] = 0.;
174 A[1][4] = 0.;
175 A[2][0] = 0.105701991897526;
176 A[2][1] = 0.;
177 A[2][2] = 0.894298008102474;
178 A[2][3] = 0.;
179 A[2][4] = 0.;
180 A[3][0] = 0.411551205755676;
181 A[3][1] = 0.011170516177380;
182 A[3][2] = 0.;
183 A[3][3] = 0.577278278066944;
184 A[3][4] = 0.;
185 A[4][0] = 0.186911123548222;
186 A[4][1] = 0.013354480555382;
187 A[4][2] = 0.012758264566319;
188 A[4][3] = 0.;
189 A[4][4] = 0.786976131330077;
190
191 B[0][0] = 0.418883109982196;
192 B[0][1] = 0.;
193 B[0][2] = 0.;
194 B[0][3] = 0.;
195 B[0][4] = 0.;
196 B[1][0] = 0.;
197 B[1][1] = 0.211483970024081;
198 B[1][2] = 0.;
199 B[1][3] = 0.;
200 B[1][4] = 0.;
201 B[2][0] = 0.000000000612488;
202 B[2][1] = 0.;
203 B[2][2] = 0.374606330884848;
204 B[2][3] = 0.;
205 B[2][4] = 0.;
206 B[3][0] = 0.046744815663888;
207 B[3][1] = 0.004679140556487;
208 B[3][2] = 0.;
209 B[3][3] = 0.241812120441849;
210 B[3][4] = 0.;
211 B[4][0] = 0.071938257223857;
212 B[4][1] = 0.005593966347235;
213 B[4][2] = 0.005344221539515;
214 B[4][3] = 0.;
215 B[4][4] = 0.329651009373300;
216
217 coefficients_are_initialized = true;
218 }
219 else if(stages == 6)
220 {
221 A[0][0] = 1.;
222 A[0][1] = 0.;
223 A[0][2] = 0.;
224 A[0][3] = 0.;
225 A[0][4] = 0.;
226 A[0][5] = 0.;
227 A[1][0] = 0.271376652410776;
228 A[1][1] = 0.728623347589224;
229 A[1][2] = 0.;
230 A[1][3] = 0.;
231 A[1][4] = 0.;
232 A[1][5] = 0.;
233 A[2][0] = 0.003607665467954;
234 A[2][1] = 0.;
235 A[2][2] = 0.996392334532046;
236 A[2][3] = 0.;
237 A[2][4] = 0.;
238 A[2][5] = 0.;
239 A[3][0] = 0.295174024904477;
240 A[3][1] = 0.104490494022953;
241 A[3][2] = 0.;
242 A[3][3] = 0.600335481072570;
243 A[3][4] = 0.;
244 A[3][5] = 0.;
245 A[4][0] = 0.300088895805571;
246 A[4][1] = 0.000000004174982;
247 A[4][2] = 0.000038417983374;
248 A[4][3] = 0.;
249 A[4][4] = 0.699872682036073;
250 A[4][5] = 0.;
251 A[5][0] = 0.057902281374384;
252 A[5][1] = 0.003951957060919;
253 A[5][2] = 0.179481122980769;
254 A[5][3] = 0.126656280556504;
255 A[5][4] = 0.;
256 A[5][5] = 0.632008358027424;
257
258 B[0][0] = 0.325620674236780;
259 B[0][1] = 0.;
260 B[0][2] = 0.;
261 B[0][3] = 0.;
262 B[0][4] = 0.;
263 B[0][5] = 0.;
264 B[1][0] = 0.;
265 B[1][1] = 0.237254825706663;
266 B[1][2] = 0.;
267 B[1][3] = 0.;
268 B[1][4] = 0.;
269 B[1][5] = 0.;
270 B[2][0] = 0.000014278868889;
271 B[2][1] = 0.;
272 B[2][2] = 0.324445943775684;
273 B[2][3] = 0.;
274 B[2][4] = 0.;
275 B[2][5] = 0.;
276 B[3][0] = 0.000000008816565;
277 B[3][1] = 0.034024265115088;
278 B[3][2] = 0.;
279 B[3][3] = 0.195481644115112;
280 B[3][4] = 0.;
281 B[3][5] = 0.;
282 B[4][0] = 0.0;
283 B[4][1] = 0.000000001359460;
284 B[4][2] = 0.000012509689649;
285 B[4][3] = 0.;
286 B[4][4] = 0.227893014604489;
287 B[4][5] = 0.;
288 B[5][0] = 0.033480821651945;
289 B[5][1] = 0.001286838922731;
290 B[5][2] = 0.058442764277772;
291 B[5][3] = 0.041241903471131;
292 B[5][4] = 0.;
293 B[5][5] = 0.205794987664170;
294
295 coefficients_are_initialized = true;
296 }
297 else if(stages == 7)
298 {
299 A[0][0] = 1.;
300 A[0][1] = 0.;
301 A[0][2] = 0.;
302 A[0][3] = 0.;
303 A[0][4] = 0.;
304 A[0][5] = 0.;
305 A[0][6] = 0.;
306 A[1][0] = 0.412429019730110;
307 A[1][1] = 0.587570980269890;
308 A[1][2] = 0.;
309 A[1][3] = 0.;
310 A[1][4] = 0.;
311 A[1][5] = 0.;
312 A[1][6] = 0.;
313 A[2][0] = 0.005800594241485;
314 A[2][1] = 0.0;
315 A[2][2] = 0.994199405758515;
316 A[2][3] = 0.;
317 A[2][4] = 0.;
318 A[2][5] = 0.;
319 A[2][6] = 0.;
320 A[3][0] = 0.162485678538202;
321 A[3][1] = 0.000000000270334;
322 A[3][2] = 0.;
323 A[3][3] = 0.837514321191464;
324 A[3][4] = 0.;
325 A[3][5] = 0.;
326 A[3][6] = 0.;
327 A[4][0] = 0.205239611567914;
328 A[4][1] = 0.000000000554433;
329 A[4][2] = 0.001461982584386;
330 A[4][3] = 0.;
331 A[4][4] = 0.793298405293266;
332 A[4][5] = 0.;
333 A[4][6] = 0.;
334 A[5][0] = 0.246951813330533;
335 A[5][1] = 0.000686077138452;
336 A[5][2] = 0.098274672761128;
337 A[5][3] = 0.125080337194733;
338 A[5][4] = 0.;
339 A[5][5] = 0.529007099575153;
340 A[5][6] = 0.;
341 A[6][0] = 0.003515397992512;
342 A[6][1] = 0.002051029751004;
343 A[6][2] = 0.037621575915744;
344 A[6][3] = 0.113733937331291;
345 A[6][4] = 0.000552268540167;
346 A[6][5] = 0.;
347 A[6][6] = 0.842525790469282;
348
349 B[0][0] = 0.267322588523961;
350 B[0][1] = 0.;
351 B[0][2] = 0.;
352 B[0][3] = 0.;
353 B[0][4] = 0.;
354 B[0][5] = 0.;
355 B[0][6] = 0.;
356 B[1][0] = 0.;
357 B[1][1] = 0.157070995387308;
358 B[1][2] = 0.;
359 B[1][3] = 0.;
360 B[1][4] = 0.;
361 B[1][5] = 0.;
362 B[1][6] = 0.;
363 B[2][0] = 0.019051847781300;
364 B[2][1] = 0.;
365 B[2][2] = 0.265771958656350;
366 B[2][3] = 0.;
367 B[2][4] = 0.;
368 B[2][5] = 0.;
369 B[2][6] = 0.;
370 B[3][0] = 0.014327744686556;
371 B[3][1] = 0.000000000072266;
372 B[3][2] = 0.;
373 B[3][3] = 0.223886496266790;
374 B[3][4] = 0.;
375 B[3][5] = 0.;
376 B[3][6] = 0.;
377 B[4][0] = 0.030979976588062;
378 B[4][1] = 0.000000000148213;
379 B[4][2] = 0.000390820968835;
380 B[4][3] = 0.;
381 B[4][4] = 0.212066583174926;
382 B[4][5] = 0.;
383 B[4][6] = 0.;
384 B[5][0] = 0.004054481853252;
385 B[5][1] = 0.000183403916578;
386 B[5][2] = 0.026271039908850;
387 B[5][3] = 0.033436799512346;
388 B[5][4] = 0.;
389 B[5][5] = 0.141415547205983;
390 B[5][6] = 0.;
391 B[6][0] = 0.021050441338920;
392 B[6][1] = 0.000548286582178;
393 B[6][2] = 0.010057097058147;
394 B[6][3] = 0.030403650530423;
395 B[6][4] = 0.000147633855718;
396 B[6][5] = 0.;
397 B[6][6] = 0.225226175206445;
398
399 coefficients_are_initialized = true;
400 }
401 else if(stages == 8)
402 {
403 A[0][0] = 1.;
404 A[0][1] = 0.;
405 A[0][2] = 0.;
406 A[0][3] = 0.;
407 A[0][4] = 0.;
408 A[0][5] = 0.;
409 A[0][6] = 0.;
410 A[0][7] = 0.;
411 A[1][0] = 0.108675201424538;
412 A[1][1] = 0.891324798575462;
413 A[1][2] = 0.;
414 A[1][3] = 0.;
415 A[1][4] = 0.;
416 A[1][5] = 0.;
417 A[1][6] = 0.;
418 A[1][7] = 0.;
419 A[2][0] = 0.008159777689219;
420 A[2][1] = 0.;
421 A[2][2] = 0.991840222310781;
422 A[2][3] = 0.;
423 A[2][4] = 0.;
424 A[2][5] = 0.;
425 A[2][6] = 0.;
426 A[2][7] = 0.;
427 A[3][0] = 0.000075204616622;
428 A[3][1] = 0.000017473454611;
429 A[3][2] = 0.;
430 A[3][3] = 0.999907321928768;
431 A[3][4] = 0.;
432 A[3][5] = 0.;
433 A[3][6] = 0.;
434 A[3][7] = 0.;
435 A[4][0] = 0.275083494553101;
436 A[4][1] = 0.013251614514063;
437 A[4][2] = 0.333930523474093;
438 A[4][3] = 0.;
439 A[4][4] = 0.377734367458743;
440 A[4][5] = 0.;
441 A[4][6] = 0.;
442 A[4][7] = 0.;
443 A[5][0] = 0.172210423641858;
444 A[5][1] = 0.067723791902171;
445 A[5][2] = 0.031061316699451;
446 A[5][3] = 0.018868041432255;
447 A[5][4] = 0.;
448 A[5][5] = 0.710136426324266;
449 A[5][6] = 0.;
450 A[5][7] = 0.;
451 A[6][0] = 0.155954681117895;
452 A[6][1] = 0.;
453 A[6][2] = 0.000000000164948;
454 A[6][3] = 0.000000000009768;
455 A[6][4] = 0.000000000001983;
456 A[6][5] = 0.;
457 A[6][6] = 0.844045318705406;
458 A[6][7] = 0.;
459 A[7][0] = 0.021413729448041;
460 A[7][1] = 0.000000000008708;
461 A[7][2] = 0.000000000028479;
462 A[7][3] = 0.072111559681400;
463 A[7][4] = 0.109489249417096;
464 A[7][5] = 0.046882143587611;
465 A[7][6] = 0.;
466 A[7][7] = 0.750103317828665;
467
468 B[0][0] = 0.227519284497891;
469 B[0][1] = 0.;
470 B[0][2] = 0.;
471 B[0][3] = 0.;
472 B[0][4] = 0.;
473 B[0][5] = 0.;
474 B[0][6] = 0.;
475 B[0][7] = 0.;
476 B[1][0] = 0.;
477 B[1][1] = 0.202793580427116;
478 B[1][2] = 0.;
479 B[1][3] = 0.;
480 B[1][4] = 0.;
481 B[1][5] = 0.;
482 B[1][6] = 0.;
483 B[1][7] = 0.;
484 B[2][0] = 0.026257801089696;
485 B[2][1] = 0.;
486 B[2][2] = 0.225662777716378;
487 B[2][3] = 0.;
488 B[2][4] = 0.;
489 B[2][5] = 0.;
490 B[2][6] = 0.;
491 B[2][7] = 0.;
492 B[3][0] = 0.000000000064040;
493 B[3][1] = 0.000003975547891;
494 B[3][2] = 0.;
495 B[3][3] = 0.227498198449436;
496 B[3][4] = 0.;
497 B[3][5] = 0.;
498 B[3][6] = 0.;
499 B[3][7] = 0.;
500 B[4][0] = 0.006914605853513;
501 B[4][1] = 0.003014997852682;
502 B[4][2] = 0.075975633772832;
503 B[4][3] = 0.;
504 B[4][4] = 0.085941853014477;
505 B[4][5] = 0.;
506 B[4][6] = 0.;
507 B[4][7] = 0.;
508 B[5][0] = 0.000733937709508;
509 B[5][1] = 0.015408468677066;
510 B[5][2] = 0.007067048551022;
511 B[5][3] = 0.004292843286543;
512 B[5][4] = 0.;
513 B[5][5] = 0.161569731613186;
514 B[5][6] = 0.;
515 B[5][7] = 0.;
516 B[6][0] = 0.000000000000089;
517 B[6][1] = 0.;
518 B[6][2] = 0.000000000037529;
519 B[6][3] = 0.000000000002222;
520 B[6][4] = 0.000000000000451;
521 B[6][5] = 0.;
522 B[6][6] = 0.192036586995649;
523 B[6][7] = 0.;
524 B[7][0] = 0.024945755721405;
525 B[7][1] = 0.000000000001981;
526 B[7][2] = 0.000000000006480;
527 B[7][3] = 0.016406770462739;
528 B[7][4] = 0.024910915687589;
529 B[7][5] = 0.010666591764781;
530 B[7][6] = 0.;
531 B[7][7] = 0.170662970171872;
532
533 coefficients_are_initialized = true;
534 }
535 else
536 {
537 AssertThrow(false, dealii::ExcNotImplemented());
538 }
539 }
540 else if(order == 4)
541 {
542 if(stages == 5)
543 {
544 A[0][0] = 1.;
545 A[0][1] = 0.;
546 A[0][2] = 0.;
547 A[0][3] = 0.;
548 A[0][4] = 0.;
549 A[1][0] = 0.261216512493821;
550 A[1][1] = 0.738783487506179;
551 A[1][2] = 0.;
552 A[1][3] = 0.;
553 A[1][4] = 0.;
554 A[2][0] = 0.623613752757655;
555 A[2][1] = 0.;
556 A[2][2] = 0.376386247242345;
557 A[2][3] = 0.;
558 A[2][4] = 0.;
559 A[3][0] = 0.444745181201454;
560 A[3][1] = 0.120932584902288;
561 A[3][2] = 0.;
562 A[3][3] = 0.434322233896258;
563 A[3][4] = 0.;
564 A[4][0] = 0.213357715199957;
565 A[4][1] = 0.209928473023448;
566 A[4][2] = 0.063353148180384;
567 A[4][3] = 0.;
568 A[4][4] = 0.513360663596212;
569
570 B[0][0] = 0.605491839566400;
571 B[0][1] = 0.;
572 B[0][2] = 0.;
573 B[0][3] = 0.;
574 B[0][4] = 0.;
575 B[1][0] = 0.;
576 B[1][1] = 0.447327372891397;
577 B[1][2] = 0.;
578 B[1][3] = 0.;
579 B[1][4] = 0.;
580 B[2][0] = 0.000000844149769;
581 B[2][1] = 0.;
582 B[2][2] = 0.227898801230261;
583 B[2][3] = 0.;
584 B[2][4] = 0.;
585 B[3][0] = 0.002856233144485;
586 B[3][1] = 0.073223693296006;
587 B[3][2] = 0.;
588 B[3][3] = 0.262978568366434;
589 B[3][4] = 0.;
590 B[4][0] = 0.002362549760441;
591 B[4][1] = 0.127109977308333;
592 B[4][2] = 0.038359814234063;
593 B[4][3] = 0.;
594 B[4][4] = 0.310835692561898;
595
596 coefficients_are_initialized = true;
597 }
598 else if(stages == 6)
599 {
600 A[0][0] = 1.;
601 A[0][1] = 0.;
602 A[0][2] = 0.;
603 A[0][3] = 0.;
604 A[0][4] = 0.;
605 A[0][5] = 0.;
606 A[1][0] = 0.441581886978406;
607 A[1][1] = 0.558418113021594;
608 A[1][2] = 0.;
609 A[1][3] = 0.;
610 A[1][4] = 0.;
611 A[1][5] = 0.;
612 A[2][0] = 0.496140382330059;
613 A[2][1] = 0.;
614 A[2][2] = 0.503859617669941;
615 A[2][3] = 0.;
616 A[2][4] = 0.;
617 A[2][5] = 0.;
618 A[3][0] = 0.392013998230666;
619 A[3][1] = 0.001687525300458;
620 A[3][2] = -0.;
621 A[3][3] = 0.606298476468875;
622 A[3][4] = 0.;
623 A[3][5] = 0.;
624 A[4][0] = 0.016884674246355;
625 A[4][1] = 0.000000050328214;
626 A[4][2] = 0.000018549175549;
627 A[4][3] = 0.;
628 A[4][4] = 0.983096726249882;
629 A[4][5] = 0.;
630 A[5][0] = 0.128599802059752;
631 A[5][1] = 0.150433518466544;
632 A[5][2] = 0.179199506866483;
633 A[5][3] = 0.173584325551242;
634 A[5][4] = 0.;
635 A[5][5] = 0.368182847055979;
636
637 B[0][0] = 0.448860018455995;
638 B[0][1] = 0.;
639 B[0][2] = 0.;
640 B[0][3] = 0.;
641 B[0][4] = 0.;
642 B[0][5] = 0.;
643 B[1][0] = 0.;
644 B[1][1] = 0.250651564517035;
645 B[1][2] = 0.;
646 B[1][3] = 0.;
647 B[1][4] = 0.;
648 B[1][5] = 0.;
649 B[2][0] = 0.004050697317371;
650 B[2][1] = 0.;
651 B[2][2] = 0.226162437286560;
652 B[2][3] = 0.;
653 B[2][4] = 0.;
654 B[2][5] = 0.;
655 B[3][0] = 0.000000073512372;
656 B[3][1] = 0.000757462637509;
657 B[3][2] = -0.;
658 B[3][3] = 0.272143145337661;
659 B[3][4] = 0.;
660 B[3][5] = 0.;
661 B[4][0] = 0.000592927398846;
662 B[4][1] = 0.000000022590323;
663 B[4][2] = 0.000008325983279;
664 B[4][3] = 0.;
665 B[4][4] = 0.441272814688551;
666 B[4][5] = 0.;
667 B[5][0] = 0.000000009191468;
668 B[5][1] = 0.067523591875293;
669 B[5][2] = 0.080435493959395;
670 B[5][3] = 0.077915063570602;
671 B[5][4] = 0.;
672 B[5][5] = 0.165262559524728;
673
674 coefficients_are_initialized = true;
675 }
676 else if(stages == 7)
677 {
678 A[0][0] = 1.;
679 A[0][1] = 0.;
680 A[0][2] = 0.;
681 A[0][3] = 0.;
682 A[0][4] = 0.;
683 A[0][5] = 0.;
684 A[0][6] = 0.;
685 A[1][0] = 0.277584603405600;
686 A[1][1] = 0.722415396594400;
687 A[1][2] = 0.;
688 A[1][3] = 0.;
689 A[1][4] = 0.;
690 A[1][5] = 0.;
691 A[1][6] = 0.;
692 A[2][0] = 0.528403304637363;
693 A[2][1] = 0.018109310473034;
694 A[2][2] = 0.453487384889603;
695 A[2][3] = 0.;
696 A[2][4] = 0.;
697 A[2][5] = 0.;
698 A[2][6] = 0.;
699 A[3][0] = 0.363822566916605;
700 A[3][1] = 0.025636760093079;
701 A[3][2] = 0.000072932527637;
702 A[3][3] = 0.610467740462679;
703 A[3][4] = 0.;
704 A[3][5] = 0.;
705 A[3][6] = 0.;
706 A[4][0] = 0.080433061177282;
707 A[4][1] = 0.000000001538366;
708 A[4][2] = 0.000000000000020;
709 A[4][3] = 0.000000000036824;
710 A[4][4] = 0.919566937247508;
711 A[4][5] = 0.;
712 A[4][6] = 0.;
713 A[5][0] = 0.305416318145737;
714 A[5][1] = 0.017282647045059;
715 A[5][2] = 0.214348299745317;
716 A[5][3] = 0.001174022148498;
717 A[5][4] = 0.003799138070873;
718 A[5][5] = 0.457979574844515;
719 A[5][6] = 0.;
720 A[6][0] = 0.112741543203136;
721 A[6][1] = 0.042888410429255;
722 A[6][2] = 0.185108001868376;
723 A[6][3] = 0.000003952121250;
724 A[6][4] = 0.230275526732661;
725 A[6][5] = 0.110240916986851;
726 A[6][6] = 0.318741648658470;
727
728 B[0][0] = 0.236998129331275;
729 B[0][1] = 0.;
730 B[0][2] = 0.;
731 B[0][3] = 0.;
732 B[0][4] = 0.;
733 B[0][5] = 0.;
734 B[0][6] = 0.;
735 B[1][0] = 0.001205136607466;
736 B[1][1] = 0.310012922173259;
737 B[1][2] = 0.;
738 B[1][3] = 0.;
739 B[1][4] = 0.;
740 B[1][5] = 0.;
741 B[1][6] = 0.;
742 B[2][0] = 0.000000000029361;
743 B[2][1] = 0.007771318668946;
744 B[2][2] = 0.194606801046999;
745 B[2][3] = 0.;
746 B[2][4] = 0.;
747 B[2][5] = 0.;
748 B[2][6] = 0.;
749 B[3][0] = 0.001612059039346;
750 B[3][1] = 0.011001602331536;
751 B[3][2] = 0.000031297818569;
752 B[3][3] = 0.261972390131100;
753 B[3][4] = 0.;
754 B[3][5] = 0.;
755 B[3][6] = 0.;
756 B[4][0] = 0.000000000027723;
757 B[4][1] = 0.000000000660165;
758 B[4][2] = 0.000000000000009;
759 B[4][3] = 0.000000000015802;
760 B[4][4] = 0.394617327778342;
761 B[4][5] = 0.;
762 B[4][6] = 0.;
763 B[5][0] = 0.115125889382648;
764 B[5][1] = 0.007416569384575;
765 B[5][2] = 0.091984117559200;
766 B[5][3] = 0.000503812679890;
767 B[5][4] = 0.001630338861330;
768 B[5][5] = 0.196534551952426;
769 B[5][6] = 0.;
770 B[6][0] = 0.000102167855778;
771 B[6][1] = 0.018404869978158;
772 B[6][2] = 0.079436115076445;
773 B[6][3] = 0.000001695989127;
774 B[6][4] = 0.098819030275264;
775 B[6][5] = 0.047308112450629;
776 B[6][6] = 0.136782840433305;
777
778 coefficients_are_initialized = true;
779 }
780 else if(stages == 8)
781 {
782 A[0][0] = 1.;
783 A[0][1] = 0.;
784 A[0][2] = 0.;
785 A[0][3] = 0.;
786 A[0][4] = 0.;
787 A[0][5] = 0.;
788 A[0][6] = 0.;
789 A[0][7] = 0.;
790 A[1][0] = 0.538569155333175;
791 A[1][1] = 0.461430844666825;
792 A[1][2] = 0.;
793 A[1][3] = 0.;
794 A[1][4] = 0.;
795 A[1][5] = 0.;
796 A[1][6] = 0.;
797 A[1][7] = 0.;
798 A[2][0] = 0.004485387460763;
799 A[2][1] = 0.;
800 A[2][2] = 0.995514612539237;
801 A[2][3] = 0.;
802 A[2][4] = 0.;
803 A[2][5] = 0.;
804 A[2][6] = 0.;
805 A[2][7] = 0.;
806 A[3][0] = 0.164495299288580;
807 A[3][1] = 0.016875060685979;
808 A[3][2] = 0.;
809 A[3][3] = 0.818629640025440;
810 A[3][4] = 0.;
811 A[3][5] = 0.;
812 A[3][6] = 0.;
813 A[3][7] = 0.;
814 A[4][0] = 0.426933682982668;
815 A[4][1] = 0.157047028197878;
816 A[4][2] = 0.023164224070770;
817 A[4][3] = 0.;
818 A[4][4] = 0.392855064748685;
819 A[4][5] = 0.;
820 A[4][6] = 0.;
821 A[4][7] = 0.;
822 A[5][0] = 0.082083400476958;
823 A[5][1] = 0.000000039091042;
824 A[5][2] = 0.033974171137350;
825 A[5][3] = 0.005505195713107;
826 A[5][4] = 0.;
827 A[5][5] = 0.878437193581543;
828 A[5][6] = 0.;
829 A[5][7] = 0.;
830 A[6][0] = 0.006736365648625;
831 A[6][1] = 0.010581829625529;
832 A[6][2] = 0.009353386191951;
833 A[6][3] = 0.101886062556838;
834 A[6][4] = 0.000023428364930;
835 A[6][5] = 0.;
836 A[6][6] = 0.871418927612128;
837 A[6][7] = 0.;
838 A[7][0] = 0.071115287415749;
839 A[7][1] = 0.018677648343953;
840 A[7][2] = 0.007902408660034;
841 A[7][3] = 0.319384027162348;
842 A[7][4] = 0.007121989995845;
843 A[7][5] = 0.001631615692736;
844 A[7][6] = 0.;
845 A[7][7] = 0.574167022729334;
846
847 B[0][0] = 0.282318339066479;
848 B[0][1] = 0.;
849 B[0][2] = 0.;
850 B[0][3] = 0.;
851 B[0][4] = 0.;
852 B[0][5] = 0.;
853 B[0][6] = 0.;
854 B[0][7] = 0.;
855 B[1][0] = 0.;
856 B[1][1] = 0.130270389660380;
857 B[1][2] = 0.;
858 B[1][3] = 0.;
859 B[1][4] = 0.;
860 B[1][5] = 0.;
861 B[1][6] = 0.;
862 B[1][7] = 0.;
863 B[2][0] = 0.003963092203460;
864 B[2][1] = 0.;
865 B[2][2] = 0.281052031928487;
866 B[2][3] = 0.;
867 B[2][4] = 0.;
868 B[2][5] = 0.;
869 B[2][6] = 0.;
870 B[2][7] = 0.;
871 B[3][0] = 0.000038019518678;
872 B[3][1] = 0.004764139104512;
873 B[3][2] = 0.;
874 B[3][3] = 0.231114160282572;
875 B[3][4] = 0.;
876 B[3][5] = 0.;
877 B[3][6] = 0.;
878 B[3][7] = 0.;
879 B[4][0] = 0.000019921336144;
880 B[4][1] = 0.044337256156151;
881 B[4][2] = 0.006539685265423;
882 B[4][3] = 0.;
883 B[4][4] = 0.110910189373703;
884 B[4][5] = 0.;
885 B[4][6] = 0.;
886 B[4][7] = 0.;
887 B[5][0] = 0.000000034006679;
888 B[5][1] = 0.000000011036118;
889 B[5][2] = 0.009591531566657;
890 B[5][3] = 0.001554217709960;
891 B[5][4] = 0.;
892 B[5][5] = 0.247998929466160;
893 B[5][6] = 0.;
894 B[5][7] = 0.;
895 B[6][0] = 0.013159891155054;
896 B[6][1] = 0.002987444564164;
897 B[6][2] = 0.002640632454359;
898 B[6][3] = 0.028764303955070;
899 B[6][4] = 0.000006614257074;
900 B[6][5] = 0.;
901 B[6][6] = 0.246017544274548;
902 B[6][7] = 0.;
903 B[7][0] = 0.000000010647874;
904 B[7][1] = 0.005273042658132;
905 B[7][2] = 0.002230994887525;
906 B[7][3] = 0.090167968072837;
907 B[7][4] = 0.002010668386475;
908 B[7][5] = 0.000460635032368;
909 B[7][6] = 0.;
910 B[7][7] = 0.162097880203691;
911
912 coefficients_are_initialized = true;
913 }
914 else
915 {
916 AssertThrow(false, dealii::ExcNotImplemented());
917 }
918 }
919 else
920 {
921 AssertThrow(false, dealii::ExcNotImplemented());
922 }
923
924 AssertThrow(coefficients_are_initialized, dealii::ExcNotImplemented());
925
926 // calculate coefficients c_i for non-autonomous systems, i.e.,
927 // systems with explicit time-dependency
928 dealii::FullMatrix<double> crk;
929 crk.reinit(stages, stages);
930
931 for(unsigned int k = 0; k < stages; ++k)
932 {
933 for(unsigned int i = 0; i < stages; ++i)
934 {
935 double csum = 0.;
936 for(unsigned int l = k + 1; l < i + 1; ++l)
937 {
938 csum = csum + crk[l - 1][k] * A[i][l];
939 }
940 crk[i][k] = B[i][k] + csum;
941 }
942 }
943
944 c.resize(stages);
945 c[0] = 0.;
946 for(unsigned int k = 1; k < stages; ++k)
947 {
948 c[k] = 0.;
949 for(unsigned int l = 0; l < k + 1; ++l)
950 {
951 c[k] = c[k] + crk[k - 1][l];
952 }
953 }
954}
955
956} // namespace ExaDG
957
958#endif /* INCLUDE_EXADG_TIME_INTEGRATION_SSP_RUNGE_KUTTA_H_ */
ExaDG
Definition driver.cpp:33
Generated by doxygen 1.13.2