Runge-Kutta (Fehlberg) 45 with error estimation.

Fourth order accurate evolution with an additional fifth order step for error estimation. How much sense the error makes is unclear to me, but hey.

For the moment the error is stored in an internal MoL array ErrorEstimate; there is one per evolved variable. At a later point this may be moved out to user thorns who can register their own etc.

As the implementation uses 6 evaluations of the RHS (necessary) and 6 levels of scratch space (one more than necessary - laziness kicked in) then this is very expensive.

This is a partial fix for PR/1840.


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

4 files changed, 249 insertions, 0 deletions

diff --git a/src/IndexArrays.c b/src/IndexArrays.c
index 096a745..60cd140 100644
--- a/src/IndexArrays.c
+++ b/src/IndexArrays.c
@@ -299,6 +299,10 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
   {
     sprintf(infoline, "Runge-Kutta 3");
   }
+  else if (CCTK_EQUALS(ODE_Method,"RK45")) 
+  {
+    sprintf(infoline, "Runge-Kutta 45");
+  }
   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 5316c9f..bd91e5e 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -129,6 +129,14 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
   }
 
 
+  if ( (CCTK_Equals(ODE_Method, "RK45")) &&
+       ( !((MoL_Intermediate_Steps == 6)||(MoL_Num_Scratch_Levels > 5)) ) )
+  {
+    CCTK_PARAMWARN("When using the RK45 evolver the "
+                   "number of intermediate steps must be 6 " 
+                   "and the number of scratch levels at least 6.");
+  }
+
   return 0;
   
 }
diff --git a/src/RK45.c b/src/RK45.c
new file mode 100644
index 0000000..41a0462
--- /dev/null
+++ b/src/RK45.c
@@ -0,0 +1,236 @@
+ /*@@
+   @file      RK45.c
+   @date      Sun May 26 03:47:15 2002
+   @author    Ian Hawke
+   @desc 
+   RK45 following Forsythe, Malcolm and Moler
+   (Computer Methods for Mathematical Computations).
+   @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_RK45_c);
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_RK45Add(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    MoL_RK45Add
+   @date       Sun May 26 03:50:44 2002
+   @author     Ian Hawke
+   @desc 
+   Performs a single step of a Runge-Kutta 45 type time
+   integration, storing the error estimate.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_RK45Add(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 *UpdateVar;
+  CCTK_REAL *RHSVar;
+  CCTK_REAL *ScratchVar;
+  CCTK_REAL *ErrorVar;
+  CCTK_REAL *OldVar;
+
+  CCTK_INT arrayscratchlocation;
+
+  CCTK_REAL beta, gamma, gamma_error;
+
+  static const CCTK_REAL beta_array[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 },
+    { 2.03240740740740740740740740741, -8.0, 7.17348927875243664717348927875, -0.205896686159844054580896686160, 0.0 },
+    { -0.296296296296296296296296296296, 2.0, -1.38167641325536062378167641326, 0.452972709551656920077972709552, -0.275 }
+  };
+
+  static const CCTK_REAL gamma_array[6] = 
+    { 0.118518518518518518518518518519, 
+      0.0, 
+      0.518986354775828460038986354776, 
+      0.506131490342016657806131490342, 
+      -0.18, 
+      0.0363636363636363636363636363636
+    };
+
+  static const CCTK_REAL gammastar_array[6] = 
+    { 0.115740740740740740740740740741, 
+      0.0, 
+      0.548927875243664717348927875244, 
+      0.535331384015594541910331384016, 
+      -0.2, 
+      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++)
+  {
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL *)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] = (*Original_Delta_Time) / cctkGH->cctk_timefac *
+        RHSVar[index];
+    }
+  }
+
+
+  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  {
+    OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, 
+                                           EvolvedVariableIndex[var]);
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL *)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);
+    
+    for (index = 0; index < totalsize; index++)
+    {
+      UpdateVar[index] = OldVar[index];
+      ErrorVar[index] = 0;
+    }
+    
+    if (*MoL_Intermediate_Step - 1)
+    { 
+
+      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 (scratchstep = 0; scratchstep < 6; 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, "Ian has been too lazy to write the RK45 routine "
+              "out for array variables. Better send him an email...");
+
+  }
+        
+  return;
+}
+
+/********************************************************************
+ *********************     Local Routines   *************************
+ ********************************************************************/
+
diff --git a/src/make.code.defn b/src/make.code.defn
index adfbe27..544c7bd 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -11,6 +11,7 @@ SRCS = 	ChangeType.c \
 	ParamCheck.c \
 	RK2.c \
 	RK3.c \
+	RK45.c \
 	Registration.c \
 	RKCoefficients.c \
 	RHSNaNCheck.c \