Efficient RK4, as provided by Yosef Zlochower.

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

7 files changed, 356 insertions, 0 deletions

diff --git a/param.ccl b/param.ccl
index ea83504..bd8e86f 100644
--- a/param.ccl
+++ b/param.ccl
@@ -87,6 +87,7 @@ KEYWORD ODE_Method "The ODE method use by MoL to do time integration"
   "ICN-avg"	:: "Iterative Crank Nicholson with averaging"
   "RK2"		:: "Efficient RK2"
   "RK3"		:: "Efficient RK3"
+  "RK4"		:: "Efficient RK4"
   "RK45"        :: "RK45 (Fehlberg) with error estimation"
   "RK45CK"      :: "RK45CK (Cash-Karp) with error estimation"
   "RK65"        :: "RK65 with error estimation"
diff --git a/schedule.ccl b/schedule.ccl
index aff477d..b54905c 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -558,6 +558,13 @@ else if (CCTK_Equals(ODE_Method,"RK3"))
     LANG: C
   } "Updates calculated with the efficient Runge-Kutta 3 method"
 }
+else if (CCTK_Equals(ODE_Method,"RK4"))
+{
+  schedule MoL_RK4Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
+  {
+    LANG: C
+  } "Updates calculated with the efficient Runge-Kutta 4 method"
+}
 else if (CCTK_Equals(ODE_Method,"RK45") || CCTK_Equals(ODE_Method,"RK45CK"))
 {
   STORAGE: ErrorEstimate ErrorScalars
@@ -625,6 +632,29 @@ schedule GROUP MoL_PostStep AT PostRestrict
 {
 } "Ensure that everything is correct after restriction"
 
+############################################################
+###  Additional boundary condition bins as requested by  ###
+###  Yosef Zlochower.                                    ###
+############################################################
+
+# schedule GROUP MoL_OldBdry_Wrap IN MoL_PostStep
+# {
+# } "Wrapper group, do not schedule directly into this group"
+# 
+# schedule GROUP MoL_OldStyleBoundaries in MoL_OldBdry_Wrap
+# {
+# } "Place old style boundary routines here"
+# 
+# schedule MoL_OldBdry_SetDt IN MoL_OldBdry_Wrap BEFORE MoL_OldStyleBoundaries
+# {
+#   LANGUAGE: C
+# } "Store and change dt"
+# 
+# schedule MoL_OldBdry_ResetDt IN MoL_OldBdry_Wrap AFTER MoL_OldStyleBoundaries
+# {
+#   LANGUAGE: C
+# } "Reset dt"
+
 #################################################
 ###  Final internal MoL stuff; decrement the  ###
 ###  counter, change time and timestep        ###
diff --git a/src/IndexArrays.c b/src/IndexArrays.c
index 5e59734..e30d169 100644
--- a/src/IndexArrays.c
+++ b/src/IndexArrays.c
@@ -304,6 +304,10 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
   {
     sprintf(infoline, "Runge-Kutta 3");
   }
+  else if (CCTK_EQUALS(ODE_Method,"RK4")) 
+  {
+    sprintf(infoline, "Runge-Kutta 4");
+  }
   else if (CCTK_EQUALS(ODE_Method,"RK45")) 
   {
     sprintf(infoline, "Runge-Kutta 45 (Fehlberg)");
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index c2346b0..112c530 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -128,6 +128,14 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
                    "number of intermediate steps must be 3");
   }
 
+  if ( (CCTK_Equals(ODE_Method, "RK4"))&&(!(MoL_Intermediate_Steps == 4))
+       && (!(MoL_Num_Scratch_Levels > 1)))
+  {
+    CCTK_PARAMWARN("When using the efficient RK4 evolver the "
+                   "number of intermediate steps must be 4, and"
+                   " the number of scratch levels at least 1");
+  }
+
   if ( (CCTK_Equals(ODE_Method, "RK45") || CCTK_Equals(ODE_Method, "RK45CK")) &&
        ( !((MoL_Intermediate_Steps == 6)||(MoL_Num_Scratch_Levels > 5)) ) )
   {
diff --git a/src/RK4.c b/src/RK4.c
new file mode 100644
index 0000000..940402b
--- /dev/null
+++ b/src/RK4.c
@@ -0,0 +1,296 @@
+ /*@@
+   @file      RK4.c
+   @date      Fri July 14, 2006
+   @author    Yosef Zlochower
+   @desc 
+   A routine to perform RK4 evolution. Mostly copied from
+   genericRK.c
+   @enddesc 
+   @version   $Header$
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include <stdio.h>
+#include "ExternalVariables.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusBase_MoL_RK4_c);
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_RK4Add(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    MoL_RK4Add
+   @date       
+   @author     
+   @desc 
+   Performs a single step of a RK4 type time
+   integration.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_RK4Add(CCTK_ARGUMENTS)
+{
+
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+    
+  cGroupDynamicData arraydata;
+  CCTK_INT groupindex, ierr;
+  CCTK_INT arraytotalsize, arraydim;
+
+  /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+
+CCTK_WARN(0, "not implemented");
+  CCTK_COMPLEX Complex_alpha, Complex_beta, Complex_Delta_Time;
+  CCTK_COMPLEX * restrict OldComplexVar;
+  CCTK_COMPLEX * restrict UpdateComplexVar;
+  CCTK_COMPLEX const * restrict RHSComplexVar;
+  CCTK_COMPLEX * restrict ScratchComplexVar;
+
+#endif
+
+  static CCTK_INT scratchspace_firstindex = -99;
+  CCTK_INT index, var, scratchstep, alphaindex, scratchindex;
+  CCTK_INT totalsize, singlearraysize;
+  CCTK_REAL alpha, beta;
+  CCTK_REAL * restrict UpdateVar;
+  CCTK_REAL * restrict OldVar;
+  CCTK_REAL const * restrict RHSVar;
+  CCTK_REAL * restrict ScratchVar;
+
+  CCTK_INT arrayscratchlocation;
+
+  totalsize = 1;
+  for (arraydim = 0; arraydim < cctk_dim; arraydim++)
+  {
+    totalsize *= cctk_lsh[arraydim];
+  }
+
+  if (scratchspace_firstindex == -99)
+  {
+    scratchspace_firstindex = CCTK_FirstVarIndex("MOL::SCRATCHSPACE");
+  }
+
+  switch (MoL_Intermediate_Steps - (*MoL_Intermediate_Step))
+  {
+    case 0:
+      alpha = 1.0 / 3.0;
+      beta  = 0.5;
+      break;
+    case 1:
+      alpha = 2.0 / 3.0;
+      beta  = 0.5;
+      break;
+    case 2:
+      alpha = 1.0 / 3.0;
+      beta  = 1.0;
+      break;
+    case 3:
+      alpha = 1.0;
+      beta  = 1.0 / 6.0;
+  }
+
+  /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+  
+  Complex_Delta_Time = CCTK_Cmplx((*Original_Delta_Time) / cctkGH->cctk_timefac, 0);
+  Complex_beta = CCTK_Cmplx(beta, 0);
+
+#endif
+
+  /* Real GFs */
+
+  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  {
+    
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                 RHSVariableIndex[var]);
+/* #define MOLDEBUG 1 */
+#ifdef MOLDEBUG
+    printf("In generic RK. Variable %d (%s). RHS %d (%s). beta %g.\n",
+           EvolvedVariableIndex[var],
+           CCTK_VarName(EvolvedVariableIndex[var]),
+           RHSVariableIndex[var],
+           CCTK_VarName(RHSVariableIndex[var]),
+           beta);
+#endif
+
+    for (index = 0; index < totalsize; index++)
+    {
+      UpdateVar[index] = OldVar[index] + 
+        (*Original_Delta_Time) / cctkGH->cctk_timefac * beta * RHSVar[index];
+#ifdef MOLDEBUG
+      if (CCTK_EQUALS(verbose,"extreme"))
+      {
+        printf("Variable: %d. Index: %d. dt: %f. beta %f. RHS: %f. q: %f.\n",
+               var, index, (*Original_Delta_Time) / cctkGH->cctk_timefac, beta, RHSVar[index], 
+               UpdateVar[index]);
+      }
+#endif
+    }
+    ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                      scratchspace_firstindex
+                                      + var );
+    
+    /* scratch storage */
+    if ((*MoL_Intermediate_Step) == MoL_Intermediate_Steps)
+    {
+      for (index = 0; index < totalsize; index++)
+      {
+        ScratchVar[index] = 0;
+      }
+    }
+
+    if ((*MoL_Intermediate_Step)>1)
+    {
+      for (index = 0; index < totalsize; index++)
+      {
+        ScratchVar[index] += alpha * UpdateVar[index];
+      }
+    }
+    else
+    {
+      for (index = 0; index < totalsize; index++)
+      {
+        UpdateVar[index] += ScratchVar[index] - 4.0 / 3.0 * OldVar[index];
+      }
+    }
+
+  }
+
+
+  /* Complex GFs */
+
+  /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+CCTK_WARN(0, "not done");
+#endif
+
+  /* Real arrays */
+
+  arrayscratchlocation = 0;
+
+  arraytotalsize = MoL_Max_Evolved_Array_Size;
+  if (MoL_Num_Scratch_Levels)
+  {
+    groupindex = CCTK_GroupIndex("MOL::ARRAYSCRATCHSPACE");
+    ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
+                                 &arraydata);
+    if (ierr)
+    {
+      CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, 
+                 "The driver does not return group information "
+                 "for group '%s' (%d).", 
+                 CCTK_GroupName(groupindex), ierr);
+    }
+    arraytotalsize = 1;
+    for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
+    {
+      arraytotalsize *= arraydata.lsh[arraydim];
+    }
+  }
+  singlearraysize = arraytotalsize / MoLNumEvolvedArrayVariables;
+
+  for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
+  {
+    
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedArrayVariableIndex[var]);
+    OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, 
+                                              EvolvedArrayVariableIndex[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                 RHSArrayVariableIndex[var]);
+    
+    groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
+    ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
+                                 &arraydata);
+    if (ierr)
+    {
+      CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, 
+                 "The driver does not return group information "
+                 "for group '%s'.", 
+                 CCTK_GroupName(groupindex));
+    }
+    arraytotalsize = 1;
+    for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
+    {
+      arraytotalsize *= arraydata.lsh[arraydim];
+    }
+
+    ScratchVar = &ArrayScratchSpace[arrayscratchlocation];
+
+    for (index = 0; index < arraytotalsize; index++)
+    {
+      UpdateVar[index] = OldVar[index] +
+        (*Original_Delta_Time) / cctkGH->cctk_timefac * beta * RHSVar[index];
+    }
+
+    if ((*MoL_Intermediate_Step) == MoL_Intermediate_Steps)
+    {
+      for (index = 0; index < arraytotalsize; index++)
+      {
+        ScratchVar[index] = 0;
+      }
+    }
+    
+    if ((*MoL_Intermediate_Step)>1)
+    {
+      for (index = 0; index < arraytotalsize; index++)
+      {
+        ScratchVar[index] += alpha * UpdateVar[index];
+      }
+    }
+    else
+    {
+      for (index = 0; index < arraytotalsize; index++)
+      {
+        UpdateVar[index] += ScratchVar[index] - 4.0 / 3.0 * OldVar[index];
+      }
+    }
+    arrayscratchlocation +=  arraytotalsize;
+  }
+
+  return;
+}
diff --git a/src/SetTime.c b/src/SetTime.c
index 3e9ae3f..5149eb2 100644
--- a/src/SetTime.c
+++ b/src/SetTime.c
@@ -247,6 +247,22 @@ int MoL_ResetTime(CCTK_ARGUMENTS)
         0.5*(*Original_Delta_Time)/cctkGH->cctk_timefac;
     }
   }
+  else if (CCTK_EQUALS(ODE_Method,"RK4"))
+  {
+    const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+
+    CCTK_REAL dt = (*Original_Delta_Time)/cctkGH->cctk_timefac;
+    switch (substep)
+    {
+      case 1:
+      case 2:
+        dt *= 0.5;
+        break;
+      default:
+        dt = 0;
+    }
+    cctkGH->cctk_time = (*Original_Time) - dt;
+  }
   else if (CCTK_EQUALS(ODE_Method,"RK45") || CCTK_EQUALS(ODE_Method,"RK45CK"))
   {
     const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
diff --git a/src/make.code.defn b/src/make.code.defn
index f62902e..926dbbb 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -11,6 +11,7 @@ SRCS = 	ChangeType.c \
 	ParamCheck.c \
 	RK2.c \
 	RK3.c \
+	RK4.c \
 	RK45.c \
 	RK65.c \
 	RK87.c \