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

4 files changed, 103 insertions, 16 deletions

diff --git a/src/IndexArrays.c b/src/IndexArrays.c
index 505b1ab..5e59734 100644
--- a/src/IndexArrays.c
+++ b/src/IndexArrays.c
@@ -306,7 +306,11 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
   }
   else if (CCTK_EQUALS(ODE_Method,"RK45")) 
   {
-    sprintf(infoline, "Runge-Kutta 45");
+    sprintf(infoline, "Runge-Kutta 45 (Fehlberg)");
+  }
+  else if (CCTK_EQUALS(ODE_Method,"RK45CK")) 
+  {
+    sprintf(infoline, "Runge-Kutta 45 (Cash-Karp)");
   }
   else if (CCTK_EQUALS(ODE_Method,"RK65")) 
   {
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index 272117b..c2346b0 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -128,10 +128,10 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
                    "number of intermediate steps must be 3");
   }
 
-  if ( (CCTK_Equals(ODE_Method, "RK45")) &&
+  if ( (CCTK_Equals(ODE_Method, "RK45") || CCTK_Equals(ODE_Method, "RK45CK")) &&
        ( !((MoL_Intermediate_Steps == 6)||(MoL_Num_Scratch_Levels > 5)) ) )
   {
-    CCTK_PARAMWARN("When using the RK45 evolver the "
+    CCTK_PARAMWARN("When using the RK45 or RK45CK evolver, the "
                    "number of intermediate steps must be 6 " 
                    "and the number of scratch levels at least 6.");
   }
@@ -154,13 +154,13 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
 
   if (adaptive_stepsize)
   {
-    if (CCTK_Equals(ODE_Method, "RK45")||CCTK_Equals(ODE_Method, "RK65")||CCTK_Equals(ODE_Method, "RK87"))
+    if (CCTK_Equals(ODE_Method, "RK45")||CCTK_Equals(ODE_Method, "RK45CK")||CCTK_Equals(ODE_Method, "RK65")||CCTK_Equals(ODE_Method, "RK87"))
     {
       /* everything is fine, do nothing */
     }
     else
     {
-      CCTK_PARAMWARN("Adaptive time step sizes are only possible with the RK45, RK65 and RK87 solvers");
+      CCTK_PARAMWARN("Adaptive time step sizes are only possible with the RK45, RK45CK, RK65, and RK87 solvers");
     }
   }
 
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++)
   {
diff --git a/src/SetTime.c b/src/SetTime.c
index 1d1596d..3e9ae3f 100644
--- a/src/SetTime.c
+++ b/src/SetTime.c
@@ -47,8 +47,8 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS);
  *********************     Local Data   *****************************
  ********************************************************************/
 
-/* RK45 coefficients */
-static const CCTK_REAL alpha_array[6] = {
+/* RK45 Fehlberg coefficients */
+static const CCTK_REAL alpha_array_F[6] = {
   0.0,
   1.0/4.0,
   3.0/8.0,
@@ -57,6 +57,16 @@ static const CCTK_REAL alpha_array[6] = {
   1.0/2.0,
 };
 
+/* RK45 Cash-Karp coefficients */
+static const CCTK_REAL alpha_array_CK[6] = {
+  0.0,
+  1.0/5.0,
+  3.0/10.0,
+  3.0/5.0,
+  1.0,
+  7.0/8.0,
+};
+
 /* RK65 coefficients */
 static const CCTK_REAL alpha_array65[8] = {
   0.0,
@@ -237,9 +247,22 @@ int MoL_ResetTime(CCTK_ARGUMENTS)
         0.5*(*Original_Delta_Time)/cctkGH->cctk_timefac;
     }
   }
-  else if (CCTK_EQUALS(ODE_Method,"RK45"))
+  else if (CCTK_EQUALS(ODE_Method,"RK45") || CCTK_EQUALS(ODE_Method,"RK45CK"))
   {
     const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+    const CCTK_REAL * alpha_array;
+    if (CCTK_EQUALS(ODE_Method, "RK45"))
+    {
+      alpha_array = alpha_array_F;
+    }
+    else if (CCTK_EQUALS(ODE_Method, "RK45CK"))
+    {
+      alpha_array = alpha_array_CK;
+    }
+    else
+    {
+      CCTK_WARN (0, "internal error");
+    }
     cctkGH->cctk_time
       = ((* Original_Time)
          + ((alpha_array[substep] - 1)
@@ -339,9 +362,22 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS)
       cctkGH->cctk_delta_time = 2.0/3.0*(*Original_Delta_Time);
     }
   }
-  else if (CCTK_EQUALS(ODE_Method,"RK45"))
+  else if (CCTK_EQUALS(ODE_Method,"RK45") || CCTK_EQUALS(ODE_Method,"RK45CK"))
   {
     const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+    const CCTK_REAL * alpha_array;
+    if (CCTK_EQUALS(ODE_Method, "RK45"))
+    {
+      alpha_array = alpha_array_F;
+    }
+    else if (CCTK_EQUALS(ODE_Method, "RK45CK"))
+    {
+      alpha_array = alpha_array_CK;
+    }
+    else
+    {
+      CCTK_WARN (0, "internal error");
+    }
     cctkGH->cctk_delta_time
       = ((alpha_array[substep + 1] - alpha_array[substep])
          * (* Original_Delta_Time));