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

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++)
   {