Skip to content

Improve precision of some RK coefficients #675

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Mar 21, 2025
Merged

Improve precision of some RK coefficients #675

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
93eb45d
Improve precision of some RK coefficients
Steven-Roberts Mar 20, 2025
0b4b3de
update outputs
gardner48 Mar 21, 2025
684a74d
Update single precision tests
Steven-Roberts Mar 21, 2025
59aec10
update answers submodule
gardner48 Mar 21, 2025
dd9e0f7
remove trailing whitespace
gardner48 Mar 21, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,10 @@

### New Features and Enhancements

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`.

The Soderlind time step adaptivity controller now starts with an I controller
until there is sufficient history of past time steps and errors.

Expand Down
4 changes: 2 additions & 2 deletions doc/arkode/guide/source/Butcher.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -692,8 +692,8 @@ Accessible via the string ``"ARKODE_SAYFY_ABURUB_6_3_4"`` to
\frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 \\
1 & -1 & 2 & 0 & 0 & 0 & 0 \\
1 & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} & 0 & 0 & 0 \\
\frac{1}{2} & 0.137 & 0.226 & 0.137 & 0 & 0 & 0 \\
1 & 0.452 & -0.904 & -0.548 & 0 & 2 & 0 \\
\frac{1}{2} & \frac{137}{1000} & \frac{113}{500} & \frac{137}{1000} & 0 & 0 & 0 \\
1 & \frac{113}{250} & -\frac{113}{125} & -\frac{137}{250} & 0 & 2 & 0 \\
\hline
4 & \frac{1}{6} & \frac{1}{3} & \frac{1}{12} & 0 & \frac{1}{3} & \frac{1}{12} \\
3 & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} & 0 & 0 & 0
Expand Down
4 changes: 4 additions & 0 deletions doc/shared/RecentChanges.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,10 @@

**New Features and Enhancements**

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``.

The Soderlind time step adaptivity controller now starts with an I controller
until there is sufficient history of past time steps and errors.

Expand Down
20 changes: 10 additions & 10 deletions examples/arkode/CXX_serial/ark_kpr_Mt_2_8_0_-10.out
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ Nonlinear Kvaerno-Prothero-Robinson test problem with mass matrix:
-2.600000 0.756013 1.218400 8.59e-10 7.76e-11
-2.500000 0.774227 1.183861 4.35e-10 5.23e-11
-2.400000 0.794546 1.150885 7.76e-11 2.72e-11
-2.300000 0.816616 1.119953 2.03e-10 1.67e-12
-2.300000 0.816616 1.119953 2.03e-10 1.66e-12
-2.200000 0.840089 1.091560 4.11e-10 2.58e-11
-2.100000 0.864625 1.066204 5.56e-10 5.68e-11
-2.000000 0.889903 1.044367 6.53e-10 9.23e-11
Expand Down Expand Up @@ -55,20 +55,20 @@ Nonlinear Kvaerno-Prothero-Robinson test problem with mass matrix:
1.200000 1.086821 1.712320 3.71e-10 7.19e-12
1.300000 1.064777 1.721499 3.96e-10 4.90e-12
1.400000 1.041625 1.727845 4.23e-10 2.64e-12
1.500000 1.017531 1.731328 4.53e-10 4.22e-13
1.500000 1.017531 1.731328 4.53e-10 4.20e-13
1.600000 0.992673 1.731928 4.82e-10 1.70e-12
1.700000 0.967253 1.729643 5.11e-10 3.62e-12
1.800000 0.941488 1.724485 5.35e-10 5.17e-12
1.700000 0.967253 1.729643 5.11e-10 3.63e-12
1.800000 0.941488 1.724485 5.35e-10 5.18e-12
1.900000 0.915617 1.716479 5.49e-10 6.09e-12
2.000000 0.889903 1.705666 5.46e-10 5.99e-12
2.100000 0.864625 1.692102 5.18e-10 4.38e-12
2.200000 0.840089 1.675857 4.51e-10 6.23e-13
2.300000 0.816616 1.657017 3.34e-10 6.00e-12
2.200000 0.840089 1.675857 4.51e-10 6.27e-13
2.300000 0.816616 1.657017 3.34e-10 5.99e-12
2.400000 0.794546 1.635684 1.52e-10 1.62e-11
2.500000 0.774227 1.611978 1.06e-10 3.05e-11
2.600000 0.756013 1.586033 4.44e-10 4.92e-11
2.700000 0.740246 1.558005 8.51e-10 7.15e-11
2.800000 0.727247 1.528067 1.30e-09 9.60e-11
2.800000 0.727247 1.528067 1.30e-09 9.59e-11
2.900000 0.717301 1.496412 1.75e-09 1.20e-10
3.000000 0.710636 1.463257 2.13e-09 1.41e-10
3.100000 0.707412 1.428839 2.39e-09 1.55e-10
Expand All @@ -95,7 +95,7 @@ Nonlinear Kvaerno-Prothero-Robinson test problem with mass matrix:
5.200000 1.110972 1.056667 4.85e-10 1.32e-10
5.300000 1.130127 1.080617 4.29e-10 8.39e-11
5.400000 1.147757 1.107807 3.81e-10 4.08e-11
5.500000 1.163759 1.137743 3.41e-10 5.62e-12
5.500000 1.163759 1.137743 3.41e-10 5.63e-12
5.600000 1.178042 1.169929 3.10e-10 2.06e-11
5.700000 1.190528 1.203875 2.87e-10 3.82e-11
5.800000 1.201149 1.239112 2.70e-10 4.87e-11
Expand All @@ -104,7 +104,7 @@ Nonlinear Kvaerno-Prothero-Robinson test problem with mass matrix:
6.100000 1.221325 1.348272 2.48e-10 5.30e-11
6.200000 1.224039 1.384525 2.46e-10 4.98e-11
6.300000 1.224716 1.420146 2.47e-10 4.58e-11
6.400000 1.223353 1.454836 2.49e-10 4.15e-11
6.400000 1.223353 1.454836 2.49e-10 4.14e-11
6.500000 1.219956 1.488328 2.52e-10 3.71e-11
6.600000 1.214544 1.520375 2.56e-10 3.30e-11
6.700000 1.207142 1.550758 2.62e-10 2.91e-11
Expand All @@ -120,4 +120,4 @@ Final Solver Statistics:
Total mass matrix setups = 1314
Total mass matrix solves = 1716
Total mass times evals = 101
Errors: u = 7.987e-10, v = 1.04573e-10, total = 5.69586e-10
Errors: u = 7.98701e-10, v = 1.04573e-10, total = 5.69587e-10
98 changes: 45 additions & 53 deletions src/arkode/arkode_butcher_dirk.def
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@
ARKODE_BILLINGTON_3_3_2 SDIRK N N Y
ARKODE_TRBDF2_3_3_2 ESDIRK N N Y
ARKODE_KVAERNO_4_2_3 ESDIRK Y Y Y
ARKODE_ARK324L2SA_DIRK_4_2_3* ESDIRK Y Y N
ARKODE_ARK324L2SA_DIRK_4_2_3* ESDIRK Y Y Y
ARKODE_ESDIRK324L2SA_4_2_3 ESDIRK Y Y Y
ARKODE_ESDIRK325L2SA_5_2_3 ESDIRK Y Y Y
ARKODE_ESDIRK32I5L2SA_5_2_3 ESDIRK Y Y Y
Expand Down Expand Up @@ -246,28 +246,28 @@ ARK_BUTCHER_TABLE(ARKODE_ARK324L2SA_DIRK_4_2_3, { /* ARK3(2)4L[2]SA-ESDIRK */
ARKodeButcherTable B = ARKodeButcherTable_Alloc(4, SUNTRUE);
B->q = 3;
B->p = 2;
B->A[1][0] = SUN_RCONST(1767732205903.0)/SUN_RCONST(4055673282236.0);
B->A[1][1] = SUN_RCONST(1767732205903.0)/SUN_RCONST(4055673282236.0);
B->A[2][0] = SUN_RCONST(2746238789719.0)/SUN_RCONST(10658868560708.0);
B->A[2][1] = SUN_RCONST(-640167445237.0)/SUN_RCONST(6845629431997.0);
B->A[2][2] = SUN_RCONST(1767732205903.0)/SUN_RCONST(4055673282236.0);
B->A[3][0] = SUN_RCONST(1471266399579.0)/SUN_RCONST(7840856788654.0);
B->A[3][1] = SUN_RCONST(-4482444167858.0)/SUN_RCONST(7529755066697.0);
B->A[3][2] = SUN_RCONST(11266239266428.0)/SUN_RCONST(11593286722821.0);
B->A[3][3] = SUN_RCONST(1767732205903.0)/SUN_RCONST(4055673282236.0);

B->b[0] = SUN_RCONST(1471266399579.0)/SUN_RCONST(7840856788654.0);
B->b[1] = SUN_RCONST(-4482444167858.0)/SUN_RCONST(7529755066697.0);
B->b[2] = SUN_RCONST(11266239266428.0)/SUN_RCONST(11593286722821.0);
B->b[3] = SUN_RCONST(1767732205903.0)/SUN_RCONST(4055673282236.0);

B->d[0] = SUN_RCONST(2756255671327.0)/SUN_RCONST(12835298489170.0);
B->d[1] = SUN_RCONST(-10771552573575.0)/SUN_RCONST(22201958757719.0);
B->d[2] = SUN_RCONST(9247589265047.0)/SUN_RCONST(10645013368117.0);
B->d[3] = SUN_RCONST(2193209047091.0)/SUN_RCONST(5459859503100.0);

B->c[1] = SUN_RCONST(1767732205903.0)/SUN_RCONST(2027836641118.0);
B->c[2] = SUN_RCONST(3.0)/SUN_RCONST(5.0);
B->A[1][0] = SUN_RCONST(0.4358665215084589994160194511935568425293);
B->A[1][1] = SUN_RCONST(0.4358665215084589994160194511935568425293);
B->A[2][0] = SUN_RCONST(0.2576482460664272457999960162840797092643);
B->A[2][1] = SUN_RCONST(-0.09351476757488624521601546747763655179361);
B->A[2][2] = SUN_RCONST(0.4358665215084589994160194511935568425293);
B->A[3][0] = SUN_RCONST(0.1876410243467238251612921441668043913795);
B->A[3][1] = SUN_RCONST(-0.5952974735769549480478230275858851737782);
B->A[3][2] = SUN_RCONST(0.9717899277217721234705114322255239398694);
B->A[3][3] = SUN_RCONST(0.4358665215084589994160194511935568425293);

B->b[0] = SUN_RCONST(0.1876410243467238251612921441668043913795);
B->b[1] = SUN_RCONST(-0.5952974735769549480478230275858851737782);
B->b[2] = SUN_RCONST(0.9717899277217721234705114322255239398694);
B->b[3] = SUN_RCONST(0.4358665215084589994160194511935568425293);

B->d[0] = SUN_RCONST(0.2147402862233891404862383406484193714659);
B->d[1] = SUN_RCONST(-0.4851622638849390928209050808398155895845);
B->d[2] = SUN_RCONST(0.86872500252038755116621237682951240796);
B->d[3] = SUN_RCONST(0.4016969751411624011684543633618838101586);

B->c[1] = SUN_RCONST(0.8717330430169179988320389023871136850586);
B->c[2] = SUN_RCONST(0.6);
B->c[3] = SUN_RCONST(1.0);
return B;
})
Expand Down Expand Up @@ -693,40 +693,32 @@ ARK_BUTCHER_TABLE(ARKODE_ARK548L2SAb_DIRK_8_4_5, { /* ARK5(4)8L[2]SAb-ESDIRK */
})

ARK_BUTCHER_TABLE(ARKODE_ESDIRK324L2SA_4_2_3, { /* ESDIRK3(2)4L[2]SA (A,L stable) */
const sunrealtype g = SUN_RCONST(0.435866521508458999416019451193556843);
const sunrealtype g2 = g * g;
const sunrealtype g3 = g2 * g;
const sunrealtype g4 = g3 * g;
const sunrealtype g5 = g4 * g;
const sunrealtype c3 = SUN_RCONST(0.6);

ARKodeButcherTable B = ARKodeButcherTable_Alloc(4, SUNTRUE);
B->q = 3;
B->p = 2;

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));
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));
B->b[3] = g;
B->b[0] = SUN_RCONST(1.0) - g - B->b[1] - B->b[2];

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))
+ (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));
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));
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);
B->d[0] = SUN_RCONST(1.0) - B->d[1] - B->d[2] - B->d[3];

B->A[1][0] = g;
B->A[1][1] = g;
B->A[2][1] = c3*(c3-SUN_RCONST(2.0)*g)/(SUN_RCONST(4.0)*g);
B->A[2][0] = c3 - g - B->A[2][1];
B->A[2][2] = g;
B->A[3][0] = B->b[0];
B->A[3][1] = B->b[1];
B->A[3][2] = B->b[2];
B->A[3][3] = B->b[3];

B->c[1] = SUN_RCONST(2.0)*g;
B->c[2] = SUN_RCONST(3.0)/SUN_RCONST(5.0);
B->A[1][0] = SUN_RCONST(0.4358665215084589994160194511935568425293);
B->A[1][1] = SUN_RCONST(0.4358665215084589994160194511935568425293);
B->A[2][0] = SUN_RCONST(0.2576482460664272457999960162840797092643);
B->A[2][1] = SUN_RCONST(-0.09351476757488624521601546747763655179361);
B->A[2][2] = SUN_RCONST(0.4358665215084589994160194511935568425293);
B->A[3][0] = SUN_RCONST(0.1876410243467238251612921441668043913795);
B->A[3][1] = SUN_RCONST(-0.5952974735769549480478230275858851737782);
B->A[3][2] = SUN_RCONST(0.9717899277217721234705114322255239398694);
B->A[3][3] = SUN_RCONST(0.4358665215084589994160194511935568425293);

B->b[0] = SUN_RCONST(0.1876410243467238251612921441668043913795);
B->b[1] = SUN_RCONST(-0.5952974735769549480478230275858851737782);
B->b[2] = SUN_RCONST(0.9717899277217721234705114322255239398694);
B->b[3] = SUN_RCONST(0.4358665215084589994160194511935568425293);

B->d[0] = SUN_RCONST(0.1088966176158644541561307380704960821824);
B->d[1] = SUN_RCONST(-0.9153258118707127534816380978168183454991);
B->d[2] = SUN_RCONST(1.271273597302152167844715894135642876535);
B->d[3] = SUN_RCONST(0.5351555969526961314807914656106793867813);

B->c[1] = SUN_RCONST(0.8717330430169179988320389023871136850586);
B->c[2] = SUN_RCONST(0.6);
B->c[3] = SUN_RCONST(1.0);
return B;
})
Expand Down
Loading