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

1 files changed, 21 insertions, 3 deletions

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)
     {