Adaptive step size control using RK45. Due to Erik Schnetter.

Note that if you want to use this with Carpet you currently have to
use the development (darcs) version of Carpet together with the
parameter

carpet::adaptive_stepsize = "yes".


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

6 files changed, 480 insertions, 20 deletions

diff --git a/interface.ccl b/interface.ccl
index f9f45f5..40e1c79 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -194,6 +194,7 @@ CCTK_REAL RKBetaCoefficients TYPE=ARRAY DIM=1 SIZE=MoL_Intermediate_Steps DISTRI
 CCTK_INT MoL_Counters TYPE = SCALAR
 {
   MoL_Intermediate_Step
+  MoL_Stepsize_Bad
 } "The counter for the time integration method"
 
 CCTK_REAL MoL_Original_Time TYPE = SCALAR
@@ -221,6 +222,13 @@ CCTK_REAL ArraySandRScratchSpace TYPE = ARRAY DIM = 1 SIZE = MoL_Max_Evolved_Arr
 
 #CCTK_COMPLEX ComplexArraySandRScratchSpace TYPE = ARRAY DIM = 1 SIZE = MoL_Max_Evolved_ComplexArray_Size Timelevels = 1
 
-#Error vector for RK45
+#Error vector and scalars for RK45
 
 CCTK_REAL ErrorEstimate[MoL_Num_Evolved_Vars] TYPE=GF Timelevels = 1 tags='Prolongation="None"'
+
+CCTK_REAL ErrorScalars TYPE=SCALAR
+{
+  Error
+  Count
+  EstimatedDt
+} "Global error estimate"
diff --git a/param.ccl b/param.ccl
index a65216a..d9b6265 100644
--- a/param.ccl
+++ b/param.ccl
@@ -132,6 +132,44 @@ BOOLEAN disable_prolongation "If Mesh refinement is enabled should we use buffer
 {
 } "yes"
 
+
+
+BOOLEAN adaptive_stepsize "Choose the time step size adaptively"
+{
+} "no"
+
+REAL maximum_absolute_error "Maximum allowed absolute error for adaptive stepsize control"
+{
+  (0.0:*) :: ""
+} 1.0e-6
+
+REAL maximum_relative_error "Maximum allowed relative error for adaptive stepsize control"
+{
+  (0.0:*) :: ""
+} 1.0e-6
+
+REAL RHS_error_weight "Weight of the RHS in the relative error calculation"
+{
+  0.0:* :: "should be between zero and one"
+} 1.0
+
+REAL safety_factor "Safety factor for stepsize control"
+{
+  (0.0:*) :: "should be less than one"
+} 0.9
+
+REAL maximum_decrease "Maximum stepsize decrease factor"
+{
+  (0.0:*) :: "should be larger than one"
+} 10.0
+
+REAL maximum_increase "Maximum stepsize increase factor"
+{
+  (0.0:*) :: "should be larger than one"
+} 5.0
+
+
+
 KEYWORD verbose "How verbose should MoL be?"
 {
   "none"     :: "No output at all (not implemented)"
diff --git a/schedule.ccl b/schedule.ccl
index 43dc1f9..4e0ae19 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -130,6 +130,22 @@ if (initial_data_is_crap)
   } "A bad routine. Fills all previous timelevels with data copied from the current."
 }
 
+##########################################
+###  Initialise the step size control  ###
+##########################################
+
+schedule MoL_StartLoop AT Evol BEFORE MoL_Evolution
+{
+  LANG: C
+  OPTIONS: LEVEL
+} "Initialise the step size control"
+
+schedule MoL_StartLoop AT Initial
+{
+  LANG: C
+  OPTIONS: LEVEL
+} "Initialise the step size control"
+
 ######################################################
 ###  The evolution step. This is almost a self     ###
 ###  contained EVOL step with PRE and POST steps   ###
@@ -149,7 +165,7 @@ if (MoL_Num_Scratch_Levels)
     {  
       if (MoL_Max_Evolved_Array_Size)
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: ScratchSpace
 	  STORAGE: SandRScratchSpace
@@ -163,7 +179,7 @@ if (MoL_Num_Scratch_Levels)
       }
       else
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: ScratchSpace
 	  STORAGE: SandRScratchSpace
@@ -178,7 +194,7 @@ if (MoL_Num_Scratch_Levels)
     {
       if (MoL_Max_Evolved_Array_Size)
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: ScratchSpace
 	  STORAGE: ArrayScratchSpace
@@ -191,7 +207,7 @@ if (MoL_Num_Scratch_Levels)
       }
       else
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: ScratchSpace
 	  #   STORAGE: ComplexScratchSpace
@@ -208,7 +224,7 @@ if (MoL_Num_Scratch_Levels)
     {
       if (MoL_Max_Evolved_Array_Size)
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: SandRScratchSpace
 	  STORAGE: ArrayScratchSpace
@@ -221,7 +237,7 @@ if (MoL_Num_Scratch_Levels)
       }
       else
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: SandRScratchSpace
 	  #   STORAGE: ComplexScratchSpace
@@ -235,7 +251,7 @@ if (MoL_Num_Scratch_Levels)
     {
       if (MoL_Max_Evolved_Array_Size)
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: ArrayScratchSpace
 	  STORAGE: ArraySandRScratchSpace
@@ -247,7 +263,7 @@ if (MoL_Num_Scratch_Levels)
       }
       else
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  #   STORAGE: ComplexScratchSpace
 	  #   STORAGE: ComplexSandRScratchSpace
@@ -266,7 +282,7 @@ else
     {
       if (MoL_Max_Evolved_Array_Size)
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: SandRScratchSpace
 	  STORAGE: ArraySandRScratchSpace
@@ -278,7 +294,7 @@ else
       }
       else
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: SandRScratchSpace
 	  #   STORAGE: ComplexScratchSpace
@@ -292,7 +308,7 @@ else
     {
       if (MoL_Max_Evolved_Array_Size)
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: ArraySandRScratchSpace
 	  #   STORAGE: ComplexScratchSpace
@@ -303,7 +319,7 @@ else
       }
       else
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  #   STORAGE: ComplexScratchSpace
 	  #   STORAGE: ComplexSandRScratchSpace
@@ -319,7 +335,7 @@ else
     {
       if (MoL_Max_Evolved_Array_Size)
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: SandRScratchSpace
 	  STORAGE: ArraySandRScratchSpace
@@ -331,7 +347,7 @@ else
       }
       else
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: SandRScratchSpace
 	  #   STORAGE: ComplexScratchSpace
@@ -345,7 +361,7 @@ else
     {
       if (MoL_Max_Evolved_Array_Size)
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  STORAGE: ArraySandRScratchSpace
 	  #   STORAGE: ComplexScratchSpace
@@ -356,7 +372,7 @@ else
       }
       else
       {
-	schedule GROUP MoL_Evolution AT Evol
+	schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad
 	{
 	  #   STORAGE: ComplexScratchSpace
 	  #   STORAGE: ComplexSandRScratchSpace
@@ -515,7 +531,7 @@ else if (CCTK_Equals(ODE_Method,"RK3"))
 }
 else if (CCTK_Equals(ODE_Method,"RK45")) 
 {
-  STORAGE: ErrorEstimate
+  STORAGE: ErrorEstimate ErrorScalars
 
   schedule MoL_RK45Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
   {
@@ -595,6 +611,45 @@ schedule MoL_RestoreSandR IN MoL_Evolution AFTER MoL_PostStep
   LANG: C
 } "Restoring the Save and Restore variables to the original state"
 
+###################################################
+###  Loop until the step size was small enough  ###
+###################################################
+
+if (adaptive_stepsize)
+{
+  # Adaptive step size control
+  schedule MoL_InitAdaptiveError IN MoL_Evolution AFTER MoL_PostStep
+  {
+    LANG: C
+    OPTIONS: LEVEL
+  } "Control the step size: initialize error check variables"
+  
+  schedule MoL_FindAdaptiveError IN MoL_Evolution AFTER MoL_InitAdaptiveError
+  {
+    LANG: C
+  } "Control the step size: compute error check variables"
+  
+  schedule MoL_ReduceAdaptiveError IN MoL_Evolution AFTER MoL_FindAdaptiveError
+  {
+    LANG: C
+    OPTIONS: LEVEL
+  } "Control the step size: reduce error check variables"
+
+  schedule MoL_SetEstimatedDt AT POSTSTEP 
+  {
+    LANG: C
+    OPTIONS: LEVEL
+  } "Control the step size: set the new timestep"
+}
+else
+{
+  schedule MoL_FinishLoop IN MoL_Evolution AFTER MoL_PostStep
+  {
+    LANG: C
+    OPTIONS: LEVEL
+  } "Control the step size"
+}
+
 ################################################################
 ###  At the end (but before driver terminate to avoid those  ###
 ###  irritating segfaults) free the index arrays.            ###
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index bd91e5e..3bd1c21 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -128,7 +128,6 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
                    "number of intermediate steps must be 3");
   }
 
-
   if ( (CCTK_Equals(ODE_Method, "RK45")) &&
        ( !((MoL_Intermediate_Steps == 6)||(MoL_Num_Scratch_Levels > 5)) ) )
   {
@@ -137,6 +136,18 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
                    "and the number of scratch levels at least 6.");
   }
 
+  if (adaptive_stepsize)
+  {
+    if (CCTK_Equals(ODE_Method, "RK45"))
+    {
+      /* everything is fine, do nothing */
+    }
+    else
+    {
+      CCTK_PARAMWARN("Adaptive time step sizes are only possible with the RK45 solver");
+    }
+  }
+
   return 0;
   
 }
diff --git a/src/StepSize.c b/src/StepSize.c
new file mode 100644
index 0000000..61266fc
--- /dev/null
+++ b/src/StepSize.c
@@ -0,0 +1,347 @@
+ /*@@
+   @file      StepSize.c
+   @date      Tue Sep 07 2004
+   @author    Erik Schnetter
+   @desc 
+   Control the time step size.
+   @enddesc 
+   @version   $Header$
+ @@*/
+
+#include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "ExternalVariables.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusBase_MoL_StepSize_c);
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_StartLoop(CCTK_ARGUMENTS);
+
+void MoL_InitAdaptiveError(CCTK_ARGUMENTS);
+void MoL_FindAdaptiveError(CCTK_ARGUMENTS);
+void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS);
+
+void MoL_FinishLoop(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    MoL_StartLoop
+   @date       Tue Sep 07 2004
+   @author     Erik Schnetter
+   @desc 
+   Start the step size control loop, so that at least one iteration is done.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void
+MoL_StartLoop(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  *MoL_Stepsize_Bad = 1;
+
+  if (adaptive_stepsize)
+  {
+    *EstimatedDt = cctkGH->cctk_delta_time;
+  }
+  
+}
+
+ /*@@
+   @routine    MoL_InitAdaptiveError
+   @date       Tue Sep 07 2004
+   @author     Erik Schnetter
+   @desc 
+   Initialize error counters for adaptive stepsize control
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+static inline CCTK_REAL
+square (CCTK_REAL const x)
+{
+  return x * x;
+}
+
+void MoL_InitAdaptiveError(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  /* Initialise global error */
+  *Error = 0;
+  *Count = 0;
+}
+
+ /*@@
+   @routine    MoL_FindAdaptiveError
+   @date       Thu Jan 27 10:22:26 2005
+   @author     Erik Schnetter
+   @desc 
+   Compute the error local to this component/patch/...
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_FindAdaptiveError(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_REAL const * restrict UpdateVar;
+  CCTK_REAL const * restrict RHSVar;
+  CCTK_REAL const * restrict ErrorVar;
+
+  int imin[3], imax[3];
+
+  int var;
+  int i, j, k;
+  int d;
+
+  int ierr;
+
+  assert (cctk_dim <= 3);
+  for (d = 0; d < cctk_dim; d++)
+  {
+    imin[d] = cctk_bbox[2*d] ? 0 : cctk_nghostzones[d];
+    imax[d] = cctk_lsh[d] - (cctk_bbox[2*d+1] ? 0 : cctk_nghostzones[d]);
+  }
+
+  /* Calculate absolute error */
+  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  {
+
+    UpdateVar = CCTK_VarDataPtrI(cctkGH, 0, EvolvedVariableIndex[var]);
+    RHSVar = CCTK_VarDataPtrI(cctkGH, 0, RHSVariableIndex[var]);
+    ErrorVar
+      = CCTK_VarDataPtrI(cctkGH, 0, 
+                         CCTK_FirstVarIndex("MOL::ERRORESTIMATE") + var);
+
+    assert (cctk_dim == 3);
+    for (k = imin[2]; k < imax[2]; k++)
+    {
+      for (j = imin[1]; j < imax[1]; j++)
+      {
+        for (i = imin[0]; i < imax[0]; i++)
+        {
+          int const index = CCTK_GFINDEX3D(cctkGH, i, j, k);
+          CCTK_REAL const scale
+            = (square(maximum_absolute_error)
+               + square(maximum_relative_error * UpdateVar[index])
+               + square(maximum_relative_error * RHS_error_weight
+                        * (*Original_Delta_Time) * RHSVar[index]));
+          *Error += square(ErrorVar[index]) / scale;
+        }
+      }
+    }
+
+  } /* for var */
+
+  *Count
+    += (MoLNumEvolvedVariables
+        * (imax[0] - imin[0]) * (imax[1] - imin[1]) * (imax[0] - imin[0]));
+}
+
+ /*@@
+   @routine    MoL_ReduceAdaptiveError
+   @date       Thu Jan 27 10:23:14 2005
+   @author     Erik Schnetter
+   @desc 
+   Find the global error estimate. 
+   Change the timestep based on the error estimate.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  int redop;
+
+  CCTK_REAL red_local[2], red_global[2];
+
+  int ierr;
+
+  /* Get global result over all processors */
+  redop = CCTK_ReductionHandle ("sum");
+  assert (redop >= 0);
+
+  red_local[0] = *Error;
+  red_local[1] = *Count;
+  ierr = CCTK_ReduceLocArrayToArray1D
+    (cctkGH, -1, redop, red_local, red_global, 2, CCTK_VARIABLE_REAL);
+  assert (ierr == 0);
+  *Error = red_global[0];
+  *Count = red_global[1];
+
+  /* Calculate L2-norm */
+  *Error = sqrt(*Error / *Count);
+  if (! CCTK_EQUALS(verbose, "none"))
+  {
+    CCTK_VInfo (CCTK_THORNSTRING, "Integration accuracy quotient is %g", (double)*Error);
+  }
+
+  /* Decide whether to accept this step */
+  *MoL_Stepsize_Bad = *Error > 1;
+
+  if (*MoL_Stepsize_Bad)
+  {
+    /* The error is too large; reject the time step and reduce the
+       step size */
+    cctkGH->cctk_time -= cctkGH->cctk_delta_time;
+    if (! CCTK_EQUALS(verbose, "none"))
+    {
+      CCTK_VInfo (CCTK_THORNSTRING, "*** REJECTING TIME STEP ***");
+    }
+
+    cctkGH->cctk_delta_time
+      = safety_factor * (*Original_Delta_Time) / pow(*Error, 0.2);
+    if (! CCTK_EQUALS(verbose, "none"))
+    {
+      CCTK_VInfo (CCTK_THORNSTRING, "Setting time step to %g", (double)cctkGH->cctk_delta_time);
+    }
+
+    if (cctkGH->cctk_delta_time < (*Original_Delta_Time) / maximum_decrease)
+    {
+      /* No more than a factor of 10 decrease */
+      cctkGH->cctk_delta_time = (*Original_Delta_Time) / maximum_decrease;
+      if (! CCTK_EQUALS(verbose, "none"))
+      {
+        CCTK_VInfo (CCTK_THORNSTRING, "   Time step reduction too large; clamping time step to %g", (double)cctkGH->cctk_delta_time);
+      }
+    }
+
+    cctkGH->cctk_time += cctkGH->cctk_delta_time;
+  }
+  else
+  {
+    /* The error is acceptable; estimate the next step size */
+    *EstimatedDt = 
+      safety_factor * (*Original_Delta_Time) / pow(*Error, 0.25); 
+
+    if (*EstimatedDt > (*Original_Delta_Time) * maximum_increase)
+    {
+      /* No more than a factor of 5 increase */
+      *EstimatedDt = (*Original_Delta_Time) * maximum_increase;
+      if (! CCTK_EQUALS(verbose, "none"))
+      {
+        CCTK_VInfo (CCTK_THORNSTRING, "   Time step increase too large; clamping time step to %g", (double)(*EstimatedDt));
+      }
+    }
+  }
+}
+
+ /*@@
+   @routine    MoL_SetEstimatedDt
+   @date       Thu Jan 27 14:05:08 2005
+   @author     Ian Hawke
+   @desc 
+   Actually set the timestep in PostStep. 
+   This avoids problems when the timestep is changed in the middle of the
+   evolution loop.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_SetEstimatedDt(CCTK_ARGUMENTS)
+{
+
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  cctkGH->cctk_delta_time = *EstimatedDt;
+
+  if (! CCTK_EQUALS(verbose, "none"))
+  {
+    CCTK_VInfo (CCTK_THORNSTRING, "Setting time step to %g", 
+                (double)cctkGH->cctk_delta_time);
+  }
+
+}
+
+ /*@@
+   @routine    MoL_FinishLoop
+   @date       Thu Jan 27 10:24:19 2005
+   @author     Erik Schnetter
+   @desc 
+   Loop control if adaptive timestepping is not used.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_FinishLoop(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  /* Keep time step size unchanged */
+  *MoL_Stepsize_Bad = 0;
+}
+
+/********************************************************************
+ *********************     Local Routines   *************************
+ ********************************************************************/
diff --git a/src/make.code.defn b/src/make.code.defn
index 544c7bd..7683c09 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -17,7 +17,8 @@ SRCS = 	ChangeType.c \
 	RHSNaNCheck.c \
 	SandR.c \
 	SetTime.c \
-	Startup.c
+	Startup.c \
+	StepSize.c
 
 # Subdirectories containing source files
 SUBDIRS =