Peter Diener's RK65 and RK87 adaptive timestep integrators.

As yet not altered to do grid arrays.

As with RK45, adaptive timestepping does not work with mesh refinement.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@106 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b

7 files changed, 629 insertions, 5 deletions

diff --git a/src/IndexArrays.c b/src/IndexArrays.c
index 59f3950..505b1ab 100644
--- a/src/IndexArrays.c
+++ b/src/IndexArrays.c
@@ -308,6 +308,14 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
   {
     sprintf(infoline, "Runge-Kutta 45");
   }
+  else if (CCTK_EQUALS(ODE_Method,"RK65")) 
+  {
+    sprintf(infoline, "Runge-Kutta 65");
+  }
+  else if (CCTK_EQUALS(ODE_Method,"RK87")) 
+  {
+    sprintf(infoline, "Runge-Kutta 87");
+  }
   else if (CCTK_EQUALS(ODE_Method,"ICN"))
   {
     sprintf(infoline, "Iterative Crank Nicholson with %i iterations", 
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index 3bd1c21..272117b 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -136,15 +136,31 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
                    "and the number of scratch levels at least 6.");
   }
 
+  if ( (CCTK_Equals(ODE_Method, "RK65")) &&
+       ( !((MoL_Intermediate_Steps == 8)||(MoL_Num_Scratch_Levels > 7)) ) )
+  {
+    CCTK_PARAMWARN("When using the RK65 evolver the "
+                   "number of intermediate steps must be 8 " 
+                   "and the number of scratch levels at least 8.");
+  }
+
+  if ( (CCTK_Equals(ODE_Method, "RK87")) &&
+       ( !((MoL_Intermediate_Steps == 13)||(MoL_Num_Scratch_Levels > 12)) ) )
+  {
+    CCTK_PARAMWARN("When using the RK87 evolver the "
+                   "number of intermediate steps must be 13 " 
+                   "and the number of scratch levels at least 13.");
+  }
+
   if (adaptive_stepsize)
   {
-    if (CCTK_Equals(ODE_Method, "RK45"))
+    if (CCTK_Equals(ODE_Method, "RK45")||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 solver");
+      CCTK_PARAMWARN("Adaptive time step sizes are only possible with the RK45, RK65 and RK87 solvers");
     }
   }
 
diff --git a/src/RK65.c b/src/RK65.c
new file mode 100644
index 0000000..22bd26b
--- /dev/null
+++ b/src/RK65.c
@@ -0,0 +1,249 @@
+ /*@@
+   @file      RK65.c
+   @date      Sun May 26 03:47:15 2002
+   @author    Peter Diener (based on RK45.c by Ian Hawke)
+   @desc 
+   RK65 following P. J. Prince and J. R. Dormand
+   Journal of Computational and Applied Mathematics, volume 7, no 1, 1981
+   @enddesc 
+   @version   $Header$
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "ExternalVariables.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusBase_MoL_RK65_c);
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_RK65Add(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    MoL_RK65Add
+   @date       Sun May 26 03:50:44 2002
+   @author     Peter Diener (based on MoL_RK45Add by Ian Hawke)
+   @desc 
+   Performs a single step of a Runge-Kutta 65 type time
+   integration, storing the error estimate.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_RK65Add(CCTK_ARGUMENTS)
+{
+
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+    
+  cGroupDynamicData arraydata;
+  CCTK_INT groupindex, ierr;
+  CCTK_INT arraytotalsize, arraydim;
+
+  CCTK_INT index, var, scratchstep, alphaindex, scratchindex;
+  CCTK_INT totalsize;
+
+  CCTK_REAL       * restrict UpdateVar;
+  CCTK_REAL const * restrict RHSVar;
+  CCTK_REAL       * restrict ScratchVar;
+  CCTK_REAL       * restrict ErrorVar;
+  CCTK_REAL const * restrict OldVar;
+
+  CCTK_INT arrayscratchlocation;
+
+  CCTK_REAL beta, gamma, gamma_error;
+
+  static const CCTK_REAL beta_array[7][7] = {
+    { 1.0/10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+    { -2.0/81.0, 20.0/81.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+    { 615.0/1372.0, -270.0/343.0, 1053.0/1372.0, 0.0, 0.0, 0.0, 0.0 },
+    { 3243.0/5500.0, -54.0/55.0, 50949.0/71500.0, 4998.0/17875.0, 0.0, 0.0, 0.0 },
+    { -26492.0/37125.0, 72.0/55.0, 2808.0/23375.0, -24206.0/37125.0, 338.0/459.0, 0.0, 0.0 },
+    { 5561.0/2376.0, -35.0/11.0, -24117.0/31603.0, 899983.0/200772.0, -5225.0/1836.0, 3925.0/4056.0, 0.0 },
+    { 465467.0/266112.0, -2945.0/1232.0, -5610201.0/14158144.0, 10513573.0/3212352.0, -424325.0/205632.0, 376225.0/454272.0, 0.0 }
+  };
+
+  static const CCTK_REAL gamma_array[8] = 
+    { 61.0/864.0, 
+      0.0, 
+      98415.0/321776.0, 
+      16807.0/146016.0, 
+      1375.0/7344.0, 
+      1375.0/5408.0,
+      -37.0/1120.0,
+      1.0/10.0
+    };
+
+  static const CCTK_REAL gammastar_array[8] = 
+    { 821.0/10800.0,
+      0.0, 
+      19683.0/71825.0, 
+      175273.0/912600.0, 
+      395.0/3672.0, 
+      785.0/2704.0,
+      3.0/50.0,
+      0.0 
+    };
+
+  totalsize = 1;
+  for (arraydim = 0; arraydim < cctk_dim; arraydim++)
+  {
+    totalsize *= cctk_lsh[arraydim];
+  }  
+
+  /* Real GFs */
+
+  /* First store (dt times) the rhs in the scratch array. */
+
+  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  {
+    const CCTK_REAL tmp = (*Original_Delta_Time) / cctkGH->cctk_timefac;
+
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+						 RHSVariableIndex[var]);
+    ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+                                              + var 
+                                              + MoL_Num_Evolved_Vars * 
+                                              (MoL_Intermediate_Steps - 
+                                               (*MoL_Intermediate_Step)));
+    for (index = 0; index < totalsize; index++)
+    {
+      ScratchVar[index] = tmp * RHSVar[index];
+    }
+  }
+
+
+  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  {
+    OldVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 1, 
+						 EvolvedVariableIndex[var]);
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+						 CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+						 + var 
+						 + MoL_Num_Evolved_Vars * 
+						 (MoL_Intermediate_Steps - 
+						  (*MoL_Intermediate_Step)));
+    ErrorVar =  (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
+                                             CCTK_FirstVarIndex("MOL::ERRORESTIMATE")
+                                             + var);
+
+    if (*MoL_Intermediate_Step - 1)
+    { 
+
+      for (index = 0; index < totalsize; index++)
+      {
+        UpdateVar[index] = OldVar[index];
+      }
+
+      for (scratchstep = 0; 
+           scratchstep < MoL_Intermediate_Steps - (*MoL_Intermediate_Step) + 1;
+           scratchstep++)
+      {
+        
+        ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                  CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+                                                  + var 
+                                                  + MoL_Num_Evolved_Vars * scratchstep);
+        
+        beta = beta_array[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)][scratchstep];
+        
+        if ( (beta > MoL_Tiny)||(beta < -MoL_Tiny) )
+        {
+          for (index = 0; index < totalsize; index++)
+          {
+            UpdateVar[index] += beta * ScratchVar[index];
+          }
+        }
+        
+      }
+      
+    }
+    else
+    {
+
+      for (index = 0; index < totalsize; index++)
+      {
+        UpdateVar[index] = OldVar[index];
+        ErrorVar[index] = 0;
+      }
+    
+      for (scratchstep = 0; scratchstep < 8; scratchstep++)
+      {
+        
+        ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                  CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+                                                  + var 
+                                                  + MoL_Num_Evolved_Vars * scratchstep);
+        
+        gamma = gamma_array[scratchstep];
+        gamma_error = gamma - gammastar_array[scratchstep];
+        
+        if ( (gamma > MoL_Tiny)||(gamma < -MoL_Tiny) )
+        {
+          for (index = 0; index < totalsize; index++)
+          {
+            UpdateVar[index] += gamma * ScratchVar[index];
+            ErrorVar[index] += gamma_error * ScratchVar[index];
+          }
+        }
+      }
+      
+    }
+    
+  }
+
+  /* Real arrays */
+
+  arrayscratchlocation = 0;
+
+  for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
+  {
+
+    CCTK_WARN(0, "Peter has been too lazy to write the RK65 routine "
+              "out for array variables. Better send him an email...");
+
+  }
+        
+  return;
+}
+
+/********************************************************************
+ *********************     Local Routines   *************************
+ ********************************************************************/
+
diff --git a/src/RK87.c b/src/RK87.c
new file mode 100644
index 0000000..3ce19b6
--- /dev/null
+++ b/src/RK87.c
@@ -0,0 +1,272 @@
+ /*@@
+   @file      RK87.c
+   @date      Sun May 26 03:47:15 2002
+   @author    Peter Diener (based on RK45.c by Ian Hawke)
+   @desc 
+   RK87 following P. J. Prince and J. R. Dormand
+   Journal of Computational and Applied Mathematics, volume 7, no 1, 1981
+   @enddesc 
+   @version   $Header$
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "ExternalVariables.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusBase_MoL_RK87_c);
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_RK87Add(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    MoL_RK87Add
+   @date       Sun May 26 03:50:44 2002
+   @author     Peter Diener (based on MoL_RK45Add by Ian Hawke)
+   @desc 
+   Performs a single step of a Runge-Kutta 87 type time
+   integration, storing the error estimate.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_RK87Add(CCTK_ARGUMENTS)
+{
+
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+    
+  cGroupDynamicData arraydata;
+  CCTK_INT groupindex, ierr, i, j;
+  CCTK_INT arraytotalsize, arraydim;
+
+  static CCTK_INT scratchspace_firstindex = -99;
+  CCTK_INT index, var, scratchstep, alphaindex, scratchindex;
+  CCTK_INT totalsize;
+
+  CCTK_REAL       * restrict UpdateVar;
+  CCTK_REAL const * restrict RHSVar;
+  CCTK_REAL       * restrict ScratchVar;
+  CCTK_REAL       * restrict ErrorVar;
+  CCTK_REAL const * restrict OldVar;
+
+  CCTK_INT arrayscratchlocation;
+
+  CCTK_REAL beta, gamma, gamma_error;
+
+  static const CCTK_REAL beta_array[12][12] = {
+    { 1.0/18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+    { 1.0/48.0, 1.0/16.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+    { 1.0/32.0, 0.0, 3.0/32.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+    { 5.0/16.0, 0.0, -75.0/64.0, 75.0/64.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+    { 3.0/80.0, 0.0, 0.0, 3.0/16.0, 3.0/20.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+    { 29443841.0/614563906.0, 0.0, 0.0, 77736538.0/692538347.0, -28693883.0/1125000000.0, 23124283.0/1800000000.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+    { 16016141.0/946692911.0, 0.0, 0.0, 61564180.0/158732637.0, 22789713.0/633445777.0, 545815736.0/2771057229.0, -180193667.0/1043307555.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+    { 39632708.0/573591083.0, 0.0, 0.0, -433636366.0/683701615.0, -421739975.0/2616292301.0, 100302831.0/723423059.0, 790204164.0/839813087.0, 800635310.0/3783071287.0, 0.0, 0.0, 0.0, 0.0 },
+    { 246121993.0/1340847787.0, 0.0, 0.0, -37695042795.0/15268766246.0, -309121744.0/1061227803.0, -12992083.0/490766935.0, 6005943493.0/2108947869.0, 393006217.0/1396673457.0, 123872331.0/1001029789.0, 0.0, 0.0, 0.0 },
+    { -1028468189.0/846180014.0, 0.0, 0.0, 8478235783.0/508512852.0, 1311729495.0/1432422823.0, -10304129995.0/1701304382.0, -48777925059.0/3047939560.0, 15336726248.0/1032824649.0, -45442868181.0/3398467696.0, 3065993473.0/597172653.0, 0.0, 0.0 },
+    { 185892177.0/718116043.0, 0.0, 0.0, -3185094517.0/667107341.0, -477755414.0/1098053517.0, -703635378.0/230739211.0, 5731566787.0/1027545527.0, 5232866602.0/850066563.0, -4093664535.0/808688257.0, 3962137247.0/1805957418.0, 65686358.0/487910083.0, 0.0 },
+    { 403863854.0/491063109.0, 0.0, 0.0, -5068492393.0/434740067.0, -411421997.0/543043805.0, 652783627.0/914296604.0, 11173962825.0/925320556.0, -13158990841.0/6184727034.0, 3936647629.0/1978049680.0, -160528059.0/685178525.0, 248638103.0/1413531060.0, 0.0 }
+  };
+
+  static const CCTK_REAL gamma_array[13] = 
+    { 14005451.0/335480064.0,
+      0.0,
+      0.0,
+      0.0,
+      0.0,
+      -59238493.0/1068277825.0,
+      181606767.0/758867731.0,
+      561292985.0/797845732.0,
+      -1041891430.0/1371343529.0,
+      760417239.0/1151165299.0,
+      118820643.0/751138087.0,
+      -528747749.0/2220607170.0,
+      1.0/4.0
+    };
+
+  static const CCTK_REAL gammastar_array[13] = 
+    { 13451932.0/455176623.0,
+      0.0,
+      0.0,
+      0.0,
+      0.0,
+      -808719846.0/976000145.0,
+      1757004468.0/5645159321.0,
+      656045339.0/265891186.0,
+      -3867574721.0/1518517206.0,
+      465885868.0/322736535.0,
+      53011238.0/667516719.0,
+      2.0/45.0,
+      0.0
+    };
+
+  totalsize = 1;
+  for (arraydim = 0; arraydim < cctk_dim; arraydim++)
+  {
+    totalsize *= cctk_lsh[arraydim];
+  }  
+
+  if (scratchspace_firstindex == -99)
+  {
+    scratchspace_firstindex = CCTK_FirstVarIndex("MOL::SCRATCHSPACE");
+  }
+
+  /* Real GFs */
+
+  /* First store (dt times) the rhs in the scratch array. */
+
+  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  {
+    const CCTK_REAL tmp = (*Original_Delta_Time) / cctkGH->cctk_timefac;
+
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+						 RHSVariableIndex[var]);
+    ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+                                              + var 
+                                              + MoL_Num_Evolved_Vars * 
+                                              (MoL_Intermediate_Steps - 
+                                               (*MoL_Intermediate_Step)));
+    for (index = 0; index < totalsize; index++)
+    {
+      ScratchVar[index] = tmp * RHSVar[index];
+    }
+
+  }
+
+
+  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  {
+    OldVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 1, 
+						 EvolvedVariableIndex[var]);
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+						 CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+						 + var 
+						 + MoL_Num_Evolved_Vars * 
+						 (MoL_Intermediate_Steps - 
+						  (*MoL_Intermediate_Step)));
+    ErrorVar =  (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
+                                             CCTK_FirstVarIndex("MOL::ERRORESTIMATE")
+                                             + var);
+
+    if (*MoL_Intermediate_Step - 1)
+    { 
+
+      for (index = 0; index < totalsize; index++)
+      {
+        UpdateVar[index] = OldVar[index];
+      }
+
+      for (scratchstep = 0; 
+           scratchstep < MoL_Intermediate_Steps - (*MoL_Intermediate_Step) + 1;
+           scratchstep++)
+      {
+        
+        ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                  CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+                                                  + var 
+                                                  + MoL_Num_Evolved_Vars * scratchstep);
+        
+        beta = beta_array[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)][scratchstep];
+        
+        if ( (beta > MoL_Tiny)||(beta < -MoL_Tiny) )
+        {
+          for (index = 0; index < totalsize; index++)
+          {
+            UpdateVar[index] += beta * ScratchVar[index];
+          }
+        }
+        
+      }
+      
+    }
+    else
+    {
+
+      for (index = 0; index < totalsize; index++)
+      {
+        UpdateVar[index] = OldVar[index];
+        ErrorVar[index] = 0;
+      }
+    
+      for (scratchstep = 0; scratchstep < 13; scratchstep++)
+      {
+        
+        ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                  CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+                                                  + var 
+                                                  + MoL_Num_Evolved_Vars * scratchstep);
+        
+        gamma = gamma_array[scratchstep];
+        gamma_error = gamma - gammastar_array[scratchstep];
+        
+        if ( (gamma > MoL_Tiny)||(gamma < -MoL_Tiny) )
+        {
+          for (index = 0; index < totalsize; index++)
+          {
+            UpdateVar[index] += gamma * ScratchVar[index];
+            ErrorVar[index] += gamma_error * ScratchVar[index];
+          }
+        }
+
+      }
+      
+    }
+    
+  }
+
+  /* Real arrays */
+
+  arrayscratchlocation = 0;
+
+  for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
+  {
+
+    CCTK_WARN(0, "Peter has been too lazy to write the RK87 routine "
+              "out for array variables. Better send him an email...");
+
+  }
+        
+  return;
+}
+
+/********************************************************************
+ *********************     Local Routines   *************************
+ ********************************************************************/
+
diff --git a/src/SetTime.c b/src/SetTime.c
index af044f4..1d1596d 100644
--- a/src/SetTime.c
+++ b/src/SetTime.c
@@ -57,6 +57,35 @@ static const CCTK_REAL alpha_array[6] = {
   1.0/2.0,
 };
 
+/* RK65 coefficients */
+static const CCTK_REAL alpha_array65[8] = {
+  0.0,
+  1.0/10.0,
+  2.0/9.0,
+  3.0/7.0,
+  3.0/5.0,
+  4.0/5.0,
+  1.0,
+  1.0
+};
+
+/* RK87 coefficients */
+ static const CCTK_REAL alpha_array87[13] = {
+  0.0,
+  1.0/18.0,
+  1.0/12.0,
+  1.0/8.0,
+  5.0/16.0,
+  3.0/8.0,
+  59.0/400.0,
+  93.0/200.0,
+  5490023248.0/9719169821.0,
+  13.0/20.0,
+  1201146811.0/1299019798.0,
+  1.0,
+  1.0
+};
+
 /********************************************************************
  *********************     External Routines   **********************
  ********************************************************************/
@@ -216,6 +245,22 @@ int MoL_ResetTime(CCTK_ARGUMENTS)
          + ((alpha_array[substep] - 1)
             * (* Original_Delta_Time) / cctkGH->cctk_timefac));
   }
+  else if (CCTK_EQUALS(ODE_Method,"RK65"))
+  {
+    const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+    cctkGH->cctk_time
+      = ((* Original_Time)
+         + ((alpha_array65[substep] - 1)
+            * (* Original_Delta_Time) / cctkGH->cctk_timefac));
+  }
+  else if (CCTK_EQUALS(ODE_Method,"RK87"))
+  {
+    const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+    cctkGH->cctk_time
+      = ((* Original_Time)
+         + ((alpha_array87[substep] - 1)
+            * (* Original_Delta_Time) / cctkGH->cctk_timefac));
+  }
 #ifdef MOLDEBUG
   printf("MoL has once more reset t (%d): %f.\n", 
          *MoL_Intermediate_Step, cctkGH->cctk_time);
@@ -301,6 +346,20 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS)
       = ((alpha_array[substep + 1] - alpha_array[substep])
          * (* Original_Delta_Time));
   }
+  else if (CCTK_EQUALS(ODE_Method,"RK65"))
+  {
+    const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+    cctkGH->cctk_delta_time
+      = ((alpha_array65[substep + 1] - alpha_array65[substep])
+         * (* Original_Delta_Time));
+  }
+  else if (CCTK_EQUALS(ODE_Method,"RK87"))
+  {
+    const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+    cctkGH->cctk_delta_time
+      = ((alpha_array87[substep + 1] - alpha_array87[substep])
+         * (* Original_Delta_Time));
+  }
 #ifdef MOLDEBUG
   printf("MoL has once more reset dt (%d): %f.\n", 
          *MoL_Intermediate_Step, 
diff --git a/src/StepSize.c b/src/StepSize.c
index 03d6fb9..cb5499e 100644
--- a/src/StepSize.c
+++ b/src/StepSize.c
@@ -213,7 +213,7 @@ void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS)
   int redop;
 
   CCTK_REAL red_local[2], red_global[2];
-
+  CCTK_REAL p1, p2;
   int ierr;
 
   /* Get global result over all processors */
@@ -235,6 +235,24 @@ void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS)
     CCTK_VInfo (CCTK_THORNSTRING, "Integration accuracy quotient is %g", (double)*Error);
   }
 
+  if ( CCTK_EQUALS(ODE_Method,"RK45") )
+  {
+     p1 = 0.2;
+     p2 = 0.25;
+  }
+
+  if ( CCTK_EQUALS(ODE_Method,"RK65") )
+  {
+     p1 = 1.0/6.0;
+     p2 = 0.2;
+  }
+
+  if ( CCTK_EQUALS(ODE_Method,"RK87") )
+  {
+     p1 = 0.125;
+     p2 = 1.0/7.0;
+  }
+
   /* Decide whether to accept this step */
   *MoL_Stepsize_Bad = *Error > 1;
 
@@ -249,7 +267,7 @@ void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS)
     }
 
     cctkGH->cctk_delta_time
-      = safety_factor * (*Original_Delta_Time) / pow(*Error, 0.2);
+      = safety_factor * (*Original_Delta_Time) / pow(*Error, p1);
     if (! CCTK_EQUALS(verbose, "none"))
     {
       CCTK_VInfo (CCTK_THORNSTRING, "Setting time step to %g", (double)cctkGH->cctk_delta_time);
@@ -271,7 +289,7 @@ void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS)
   {
     /* The error is acceptable; estimate the next step size */
     *EstimatedDt = 
-      safety_factor * (*Original_Delta_Time) / pow(*Error, 0.25); 
+      safety_factor * (*Original_Delta_Time) / pow(*Error, p2);
 
     if (*EstimatedDt > (*Original_Delta_Time) * maximum_increase)
     {
diff --git a/src/make.code.defn b/src/make.code.defn
index 7683c09..f62902e 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -12,6 +12,8 @@ SRCS = 	ChangeType.c \
 	RK2.c \
 	RK3.c \
 	RK45.c \
+	RK65.c \
+	RK87.c \
 	Registration.c \
 	RKCoefficients.c \
 	RHSNaNCheck.c \