aboutsummaryrefslogtreecommitdiff
path: root/src/RK45.c
diff options
context:
space:
mode:
authorschnetter <schnetter@578cdeb0-5ea1-4b81-8215-5a3b8777ee0b>2006-02-19 00:31:54 +0000
committerschnetter <schnetter@578cdeb0-5ea1-4b81-8215-5a3b8777ee0b>2006-02-19 00:31:54 +0000
commitc4896403254a6de126501a42f9cb046538a44727 (patch)
tree5044935d014dab7bca78235a52a54f64dad45402 /src/RK45.c
parenta7134900411d8adb1cb1f4d6962dfcf892983326 (diff)
Implement RK45 Cash-Karp integrator in MoL.
This integrator is similar to the existing RK45 integrator. It also supports adaptive time stepping, but uses slightly different coefficients. The Numerical Recipes say that it has "slightly better error properties". git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@110 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b
Diffstat (limited to 'src/RK45.c')
-rw-r--r--src/RK45.c61
1 files changed, 54 insertions, 7 deletions
diff --git a/src/RK45.c b/src/RK45.c
index f92159a..3ca2456 100644
--- a/src/RK45.c
+++ b/src/RK45.c
@@ -67,11 +67,9 @@ void MoL_RK45Add(CCTK_ARGUMENTS)
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
- cGroupDynamicData arraydata;
- CCTK_INT groupindex, ierr;
- CCTK_INT arraytotalsize, arraydim;
+ CCTK_INT arraydim;
- CCTK_INT index, var, scratchstep, alphaindex, scratchindex;
+ CCTK_INT index, var, scratchstep;
CCTK_INT totalsize;
CCTK_REAL * restrict UpdateVar;
@@ -84,7 +82,8 @@ void MoL_RK45Add(CCTK_ARGUMENTS)
CCTK_REAL beta, gamma, gamma_error;
- static const CCTK_REAL beta_array[5][5] = {
+ /* Fehlberg coefficients */
+ static const CCTK_REAL beta_array_F[5][5] = {
{ 0.25, 0.0, 0.0, 0.0, 0.0 },
{ 0.09375, 0.28125, 0.0, 0.0, 0.0 },
{ 0.879380974055530268548020027310, -3.27719617660446062812926718252, 3.32089212562585343650432407829, 0.0, 0.0 },
@@ -92,7 +91,7 @@ void MoL_RK45Add(CCTK_ARGUMENTS)
{ -0.296296296296296296296296296296, 2.0, -1.38167641325536062378167641326, 0.452972709551656920077972709552, -0.275 }
};
- static const CCTK_REAL gamma_array[6] =
+ static const CCTK_REAL gamma_array_F[6] =
{ 0.118518518518518518518518518519,
0.0,
0.518986354775828460038986354776,
@@ -101,7 +100,7 @@ void MoL_RK45Add(CCTK_ARGUMENTS)
0.0363636363636363636363636363636
};
- static const CCTK_REAL gammastar_array[6] =
+ static const CCTK_REAL gammastar_array_F[6] =
{ 0.115740740740740740740740740741,
0.0,
0.548927875243664717348927875244,
@@ -110,6 +109,54 @@ void MoL_RK45Add(CCTK_ARGUMENTS)
0.0
};
+ /* Cash-Karp coefficients */
+ static const CCTK_REAL beta_array_CK[5][5] = {
+ { 1.0/5.0, 0.0, 0.0, 0.0, 0.0, },
+ { 3.0/40.0, 9.0/40.0, 0.0, 0.0, 0.0, },
+ { 3.0/10.0, -9.0/10.0, 6.0/5.0, 0.0, 0.0, },
+ { -11.0/54.0, 5.0/2.0, -70.0/27.0, 35.0/27.0, 0.0, },
+ { 1631.0/55296.0, 175.0/512.0, 575.0/13824.0, 44275.0/110592.0, 253.0/4096.0, },
+ };
+
+ static const CCTK_REAL gamma_array_CK[6] = {
+ 37.0/378.0,
+ 0.0,
+ 250.0/621.0,
+ 125.0/594.0,
+ 0.0,
+ 512.0/1771.0,
+ };
+
+ static const CCTK_REAL gammastar_array_CK[6] = {
+ 2825.0/27648.0,
+ 0.0,
+ 18575.0/48384.0,
+ 13525.0/55296.0,
+ 277.0/14336.0,
+ 1.0/4.0,
+ };
+
+ const CCTK_REAL (* restrict beta_array)[5];
+ const CCTK_REAL * restrict gamma_array;
+ const CCTK_REAL * restrict gammastar_array;
+
+ if (CCTK_EQUALS(ODE_Method, "RK45"))
+ {
+ beta_array = beta_array_F;
+ gamma_array = gamma_array_F;
+ gammastar_array = gammastar_array_F;
+ }
+ else if (CCTK_EQUALS(ODE_Method, "RK45CK"))
+ {
+ beta_array = beta_array_CK;
+ gamma_array = gamma_array_CK;
+ gammastar_array = gammastar_array_CK;
+ }
+ else
+ {
+ CCTK_WARN (0, "internal error");
+ }
+
totalsize = 1;
for (arraydim = 0; arraydim < cctk_dim; arraydim++)
{