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

3 files changed, 361 insertions, 2 deletions

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 =