Skip to content

Commit 1ef089c

Browse files
Steven-Robertsgardner48
Steven-Roberts
and
gardner48
committed
Feature: Improve precision of some RK coefficients (#675)
Improved the precision of the coefficients for `ARKODE_ARK324L2SA_ERK_4_2_3`, `ARKODE_VERNER_9_5_6`, `ARKODE_VERNER_10_6_7`, `ARKODE_VERNER_13_7_8`, `ARKODE_ARK324L2SA_DIRK_4_2_3`, and `ARKODE_ESDIRK324L2SA_4_2_3`. --------- Co-authored-by: David J. Gardner <gardner48@llnl.gov>
1 parent c20033e commit 1ef089c

File tree

7 files changed

+291
-290
lines changed

7 files changed

+291
-290
lines changed

‎CHANGELOG.md‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,10 @@
66

77
### New Features and Enhancements
88

9+
Improved the precision of the coefficients for `ARKODE_ARK324L2SA_ERK_4_2_3`,
10+
`ARKODE_VERNER_9_5_6`, `ARKODE_VERNER_10_6_7`, `ARKODE_VERNER_13_7_8`,
11+
`ARKODE_ARK324L2SA_DIRK_4_2_3`, and `ARKODE_ESDIRK324L2SA_4_2_3`.
12+
913
The Soderlind time step adaptivity controller now starts with an I controller
1014
until there is sufficient history of past time steps and errors.
1115

‎doc/arkode/guide/source/Butcher.rst‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -692,8 +692,8 @@ Accessible via the string ``"ARKODE_SAYFY_ABURUB_6_3_4"`` to
692692
\frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 \\
693693
1 & -1 & 2 & 0 & 0 & 0 & 0 \\
694694
1 & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} & 0 & 0 & 0 \\
695-
\frac{1}{2} & 0.137 & 0.226 & 0.137 & 0 & 0 & 0 \\
696-
1 & 0.452 & -0.904 & -0.548 & 0 & 2 & 0 \\
695+
\frac{1}{2} & \frac{137}{1000} & \frac{113}{500} & \frac{137}{1000} & 0 & 0 & 0 \\
696+
1 & \frac{113}{250} & -\frac{113}{125} & -\frac{137}{250} & 0 & 2 & 0 \\
697697
\hline
698698
4 & \frac{1}{6} & \frac{1}{3} & \frac{1}{12} & 0 & \frac{1}{3} & \frac{1}{12} \\
699699
3 & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} & 0 & 0 & 0

‎doc/shared/RecentChanges.rst‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,10 @@
22

33
**New Features and Enhancements**
44

5+
Improved the precision of the coefficients for ``ARKODE_ARK324L2SA_ERK_4_2_3``,
6+
``ARKODE_VERNER_9_5_6``, ``ARKODE_VERNER_10_6_7``, ``ARKODE_VERNER_13_7_8``,
7+
``ARKODE_ARK324L2SA_DIRK_4_2_3``, and ``ARKODE_ESDIRK324L2SA_4_2_3``.
8+
59
The Soderlind time step adaptivity controller now starts with an I controller
610
until there is sufficient history of past time steps and errors.
711

‎examples/arkode/CXX_serial/ark_kpr_Mt_2_8_0_-10.out‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ Nonlinear Kvaerno-Prothero-Robinson test problem with mass matrix:
1717
-2.600000 0.756013 1.218400 8.59e-10 7.76e-11
1818
-2.500000 0.774227 1.183861 4.35e-10 5.23e-11
1919
-2.400000 0.794546 1.150885 7.76e-11 2.72e-11
20-
-2.300000 0.816616 1.119953 2.03e-10 1.67e-12
20+
-2.300000 0.816616 1.119953 2.03e-10 1.66e-12
2121
-2.200000 0.840089 1.091560 4.11e-10 2.58e-11
2222
-2.100000 0.864625 1.066204 5.56e-10 5.68e-11
2323
-2.000000 0.889903 1.044367 6.53e-10 9.23e-11
@@ -55,20 +55,20 @@ Nonlinear Kvaerno-Prothero-Robinson test problem with mass matrix:
5555
1.200000 1.086821 1.712320 3.71e-10 7.19e-12
5656
1.300000 1.064777 1.721499 3.96e-10 4.90e-12
5757
1.400000 1.041625 1.727845 4.23e-10 2.64e-12
58-
1.500000 1.017531 1.731328 4.53e-10 4.22e-13
58+
1.500000 1.017531 1.731328 4.53e-10 4.20e-13
5959
1.600000 0.992673 1.731928 4.82e-10 1.70e-12
60-
1.700000 0.967253 1.729643 5.11e-10 3.62e-12
61-
1.800000 0.941488 1.724485 5.35e-10 5.17e-12
60+
1.700000 0.967253 1.729643 5.11e-10 3.63e-12
61+
1.800000 0.941488 1.724485 5.35e-10 5.18e-12
6262
1.900000 0.915617 1.716479 5.49e-10 6.09e-12
6363
2.000000 0.889903 1.705666 5.46e-10 5.99e-12
6464
2.100000 0.864625 1.692102 5.18e-10 4.38e-12
65-
2.200000 0.840089 1.675857 4.51e-10 6.23e-13
66-
2.300000 0.816616 1.657017 3.34e-10 6.00e-12
65+
2.200000 0.840089 1.675857 4.51e-10 6.27e-13
66+
2.300000 0.816616 1.657017 3.34e-10 5.99e-12
6767
2.400000 0.794546 1.635684 1.52e-10 1.62e-11
6868
2.500000 0.774227 1.611978 1.06e-10 3.05e-11
6969
2.600000 0.756013 1.586033 4.44e-10 4.92e-11
7070
2.700000 0.740246 1.558005 8.51e-10 7.15e-11
71-
2.800000 0.727247 1.528067 1.30e-09 9.60e-11
71+
2.800000 0.727247 1.528067 1.30e-09 9.59e-11
7272
2.900000 0.717301 1.496412 1.75e-09 1.20e-10
7373
3.000000 0.710636 1.463257 2.13e-09 1.41e-10
7474
3.100000 0.707412 1.428839 2.39e-09 1.55e-10
@@ -95,7 +95,7 @@ Nonlinear Kvaerno-Prothero-Robinson test problem with mass matrix:
9595
5.200000 1.110972 1.056667 4.85e-10 1.32e-10
9696
5.300000 1.130127 1.080617 4.29e-10 8.39e-11
9797
5.400000 1.147757 1.107807 3.81e-10 4.08e-11
98-
5.500000 1.163759 1.137743 3.41e-10 5.62e-12
98+
5.500000 1.163759 1.137743 3.41e-10 5.63e-12
9999
5.600000 1.178042 1.169929 3.10e-10 2.06e-11
100100
5.700000 1.190528 1.203875 2.87e-10 3.82e-11
101101
5.800000 1.201149 1.239112 2.70e-10 4.87e-11
@@ -104,7 +104,7 @@ Nonlinear Kvaerno-Prothero-Robinson test problem with mass matrix:
104104
6.100000 1.221325 1.348272 2.48e-10 5.30e-11
105105
6.200000 1.224039 1.384525 2.46e-10 4.98e-11
106106
6.300000 1.224716 1.420146 2.47e-10 4.58e-11
107-
6.400000 1.223353 1.454836 2.49e-10 4.15e-11
107+
6.400000 1.223353 1.454836 2.49e-10 4.14e-11
108108
6.500000 1.219956 1.488328 2.52e-10 3.71e-11
109109
6.600000 1.214544 1.520375 2.56e-10 3.30e-11
110110
6.700000 1.207142 1.550758 2.62e-10 2.91e-11
@@ -120,4 +120,4 @@ Final Solver Statistics:
120120
Total mass matrix setups = 1314
121121
Total mass matrix solves = 1716
122122
Total mass times evals = 101
123-
Errors: u = 7.987e-10, v = 1.04573e-10, total = 5.69586e-10
123+
Errors: u = 7.98701e-10, v = 1.04573e-10, total = 5.69587e-10

‎src/arkode/arkode_butcher_dirk.def‎

Lines changed: 45 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,7 @@
4242
ARKODE_BILLINGTON_3_3_2 SDIRK N N Y
4343
ARKODE_TRBDF2_3_3_2 ESDIRK N N Y
4444
ARKODE_KVAERNO_4_2_3 ESDIRK Y Y Y
45-
ARKODE_ARK324L2SA_DIRK_4_2_3* ESDIRK Y Y N
45+
ARKODE_ARK324L2SA_DIRK_4_2_3* ESDIRK Y Y Y
4646
ARKODE_ESDIRK324L2SA_4_2_3 ESDIRK Y Y Y
4747
ARKODE_ESDIRK325L2SA_5_2_3 ESDIRK Y Y Y
4848
ARKODE_ESDIRK32I5L2SA_5_2_3 ESDIRK Y Y Y
@@ -246,28 +246,28 @@ ARK_BUTCHER_TABLE(ARKODE_ARK324L2SA_DIRK_4_2_3, { /* ARK3(2)4L[2]SA-ESDIRK */
246246
ARKodeButcherTable B = ARKodeButcherTable_Alloc(4, SUNTRUE);
247247
B->q = 3;
248248
B->p = 2;
249-
B->A[1][0] = SUN_RCONST(1767732205903.0)/SUN_RCONST(4055673282236.0);
250-
B->A[1][1] = SUN_RCONST(1767732205903.0)/SUN_RCONST(4055673282236.0);
251-
B->A[2][0] = SUN_RCONST(2746238789719.0)/SUN_RCONST(10658868560708.0);
252-
B->A[2][1] = SUN_RCONST(-640167445237.0)/SUN_RCONST(6845629431997.0);
253-
B->A[2][2] = SUN_RCONST(1767732205903.0)/SUN_RCONST(4055673282236.0);
254-
B->A[3][0] = SUN_RCONST(1471266399579.0)/SUN_RCONST(7840856788654.0);
255-
B->A[3][1] = SUN_RCONST(-4482444167858.0)/SUN_RCONST(7529755066697.0);
256-
B->A[3][2] = SUN_RCONST(11266239266428.0)/SUN_RCONST(11593286722821.0);
257-
B->A[3][3] = SUN_RCONST(1767732205903.0)/SUN_RCONST(4055673282236.0);
258-
259-
B->b[0] = SUN_RCONST(1471266399579.0)/SUN_RCONST(7840856788654.0);
260-
B->b[1] = SUN_RCONST(-4482444167858.0)/SUN_RCONST(7529755066697.0);
261-
B->b[2] = SUN_RCONST(11266239266428.0)/SUN_RCONST(11593286722821.0);
262-
B->b[3] = SUN_RCONST(1767732205903.0)/SUN_RCONST(4055673282236.0);
263-
264-
B->d[0] = SUN_RCONST(2756255671327.0)/SUN_RCONST(12835298489170.0);
265-
B->d[1] = SUN_RCONST(-10771552573575.0)/SUN_RCONST(22201958757719.0);
266-
B->d[2] = SUN_RCONST(9247589265047.0)/SUN_RCONST(10645013368117.0);
267-
B->d[3] = SUN_RCONST(2193209047091.0)/SUN_RCONST(5459859503100.0);
268-
269-
B->c[1] = SUN_RCONST(1767732205903.0)/SUN_RCONST(2027836641118.0);
270-
B->c[2] = SUN_RCONST(3.0)/SUN_RCONST(5.0);
249+
B->A[1][0] = SUN_RCONST(0.4358665215084589994160194511935568425293);
250+
B->A[1][1] = SUN_RCONST(0.4358665215084589994160194511935568425293);
251+
B->A[2][0] = SUN_RCONST(0.2576482460664272457999960162840797092643);
252+
B->A[2][1] = SUN_RCONST(-0.09351476757488624521601546747763655179361);
253+
B->A[2][2] = SUN_RCONST(0.4358665215084589994160194511935568425293);
254+
B->A[3][0] = SUN_RCONST(0.1876410243467238251612921441668043913795);
255+
B->A[3][1] = SUN_RCONST(-0.5952974735769549480478230275858851737782);
256+
B->A[3][2] = SUN_RCONST(0.9717899277217721234705114322255239398694);
257+
B->A[3][3] = SUN_RCONST(0.4358665215084589994160194511935568425293);
258+
259+
B->b[0] = SUN_RCONST(0.1876410243467238251612921441668043913795);
260+
B->b[1] = SUN_RCONST(-0.5952974735769549480478230275858851737782);
261+
B->b[2] = SUN_RCONST(0.9717899277217721234705114322255239398694);
262+
B->b[3] = SUN_RCONST(0.4358665215084589994160194511935568425293);
263+
264+
B->d[0] = SUN_RCONST(0.2147402862233891404862383406484193714659);
265+
B->d[1] = SUN_RCONST(-0.4851622638849390928209050808398155895845);
266+
B->d[2] = SUN_RCONST(0.86872500252038755116621237682951240796);
267+
B->d[3] = SUN_RCONST(0.4016969751411624011684543633618838101586);
268+
269+
B->c[1] = SUN_RCONST(0.8717330430169179988320389023871136850586);
270+
B->c[2] = SUN_RCONST(0.6);
271271
B->c[3] = SUN_RCONST(1.0);
272272
return B;
273273
})
@@ -693,40 +693,32 @@ ARK_BUTCHER_TABLE(ARKODE_ARK548L2SAb_DIRK_8_4_5, { /* ARK5(4)8L[2]SAb-ESDIRK */
693693
})
694694

695695
ARK_BUTCHER_TABLE(ARKODE_ESDIRK324L2SA_4_2_3, { /* ESDIRK3(2)4L[2]SA (A,L stable) */
696-
const sunrealtype g = SUN_RCONST(0.435866521508458999416019451193556843);
697-
const sunrealtype g2 = g * g;
698-
const sunrealtype g3 = g2 * g;
699-
const sunrealtype g4 = g3 * g;
700-
const sunrealtype g5 = g4 * g;
701-
const sunrealtype c3 = SUN_RCONST(0.6);
702-
703696
ARKodeButcherTable B = ARKodeButcherTable_Alloc(4, SUNTRUE);
704697
B->q = 3;
705698
B->p = 2;
706699

707-
B->b[1] = (-SUN_RCONST(2.0)+SUN_RCONST(3.0)*c3+SUN_RCONST(6.0)*g*(SUN_RCONST(1.0)-c3))/(SUN_RCONST(12.0)*g*(c3-SUN_RCONST(2.0)*g));
708-
B->b[2] = (SUN_RCONST(1.0)-SUN_RCONST(6.0)*g+SUN_RCONST(6.0)*g2)/(SUN_RCONST(3.0)*c3*(c3-SUN_RCONST(2.0)*g));
709-
B->b[3] = g;
710-
B->b[0] = SUN_RCONST(1.0) - g - B->b[1] - B->b[2];
711-
712-
B->d[1] = c3*(-SUN_RCONST(1.0)+SUN_RCONST(6.0)*g-SUN_RCONST(24.0)*g3+SUN_RCONST(12.0)*g4-SUN_RCONST(6.0)*g5)/(SUN_RCONST(4.0)*g*(SUN_RCONST(2.0)*g-c3)*(SUN_RCONST(1.0)-SUN_RCONST(6.0)*g+SUN_RCONST(6.0)*g2))
713-
+ (SUN_RCONST(3.0)-SUN_RCONST(27.0)*g+SUN_RCONST(68.0)*g2-SUN_RCONST(55.0)*g3+SUN_RCONST(21.0)*g4-SUN_RCONST(6.0)*g5)/(SUN_RCONST(2.0)*(SUN_RCONST(2.0)*g-c3)*(SUN_RCONST(1.0)-SUN_RCONST(6.0)*g+SUN_RCONST(6.0)*g2));
714-
B->d[2] = -g*(-SUN_RCONST(2.0)+SUN_RCONST(21.0)*g-SUN_RCONST(68.0)*g2+SUN_RCONST(79.0)*g3-SUN_RCONST(33.0)*g4+SUN_RCONST(12.0)*g5)/(c3*(c3-SUN_RCONST(2.0)*g)*(SUN_RCONST(1.0)-SUN_RCONST(6.0)*g+SUN_RCONST(6.0)*g2));
715-
B->d[3] = -SUN_RCONST(3.0)*g2*(-SUN_RCONST(1.0)+SUN_RCONST(4.0)*g-SUN_RCONST(2.0)*g2+g3)/(SUN_RCONST(1.0)-SUN_RCONST(6.0)*g+SUN_RCONST(6.0)*g2);
716-
B->d[0] = SUN_RCONST(1.0) - B->d[1] - B->d[2] - B->d[3];
717-
718-
B->A[1][0] = g;
719-
B->A[1][1] = g;
720-
B->A[2][1] = c3*(c3-SUN_RCONST(2.0)*g)/(SUN_RCONST(4.0)*g);
721-
B->A[2][0] = c3 - g - B->A[2][1];
722-
B->A[2][2] = g;
723-
B->A[3][0] = B->b[0];
724-
B->A[3][1] = B->b[1];
725-
B->A[3][2] = B->b[2];
726-
B->A[3][3] = B->b[3];
727-
728-
B->c[1] = SUN_RCONST(2.0)*g;
729-
B->c[2] = SUN_RCONST(3.0)/SUN_RCONST(5.0);
700+
B->A[1][0] = SUN_RCONST(0.4358665215084589994160194511935568425293);
701+
B->A[1][1] = SUN_RCONST(0.4358665215084589994160194511935568425293);
702+
B->A[2][0] = SUN_RCONST(0.2576482460664272457999960162840797092643);
703+
B->A[2][1] = SUN_RCONST(-0.09351476757488624521601546747763655179361);
704+
B->A[2][2] = SUN_RCONST(0.4358665215084589994160194511935568425293);
705+
B->A[3][0] = SUN_RCONST(0.1876410243467238251612921441668043913795);
706+
B->A[3][1] = SUN_RCONST(-0.5952974735769549480478230275858851737782);
707+
B->A[3][2] = SUN_RCONST(0.9717899277217721234705114322255239398694);
708+
B->A[3][3] = SUN_RCONST(0.4358665215084589994160194511935568425293);
709+
710+
B->b[0] = SUN_RCONST(0.1876410243467238251612921441668043913795);
711+
B->b[1] = SUN_RCONST(-0.5952974735769549480478230275858851737782);
712+
B->b[2] = SUN_RCONST(0.9717899277217721234705114322255239398694);
713+
B->b[3] = SUN_RCONST(0.4358665215084589994160194511935568425293);
714+
715+
B->d[0] = SUN_RCONST(0.1088966176158644541561307380704960821824);
716+
B->d[1] = SUN_RCONST(-0.9153258118707127534816380978168183454991);
717+
B->d[2] = SUN_RCONST(1.271273597302152167844715894135642876535);
718+
B->d[3] = SUN_RCONST(0.5351555969526961314807914656106793867813);
719+
720+
B->c[1] = SUN_RCONST(0.8717330430169179988320389023871136850586);
721+
B->c[2] = SUN_RCONST(0.6);
730722
B->c[3] = SUN_RCONST(1.0);
731723
return B;
732724
})

0 commit comments

Comments
 (0)