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, 249 insertions, 0 deletions

diff --git a/src/RK65.c b/src/RK65.c
new file mode 100644
index 0000000..22bd26b
--- /dev/null
+++ b/src/RK65.c
@@ -0,0 +1,249 @@
+ /*@@
+   @file      RK65.c
+   @date      Sun May 26 03:47:15 2002
+   @author    Peter Diener (based on RK45.c by Ian Hawke)
+   @desc 
+   RK65 following P. J. Prince and J. R. Dormand
+   Journal of Computational and Applied Mathematics, volume 7, no 1, 1981
+   @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_RK65_c);
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_RK65Add(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    MoL_RK65Add
+   @date       Sun May 26 03:50:44 2002
+   @author     Peter Diener (based on MoL_RK45Add by Ian Hawke)
+   @desc 
+   Performs a single step of a Runge-Kutta 65 type time
+   integration, storing the error estimate.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_RK65Add(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       * restrict UpdateVar;
+  CCTK_REAL const * restrict RHSVar;
+  CCTK_REAL       * restrict ScratchVar;
+  CCTK_REAL       * restrict ErrorVar;
+  CCTK_REAL const * restrict OldVar;
+
+  CCTK_INT arrayscratchlocation;
+
+  CCTK_REAL beta, gamma, gamma_error;
+
+  static const CCTK_REAL beta_array[7][7] = {
+    { 1.0/10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+    { -2.0/81.0, 20.0/81.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+    { 615.0/1372.0, -270.0/343.0, 1053.0/1372.0, 0.0, 0.0, 0.0, 0.0 },
+    { 3243.0/5500.0, -54.0/55.0, 50949.0/71500.0, 4998.0/17875.0, 0.0, 0.0, 0.0 },
+    { -26492.0/37125.0, 72.0/55.0, 2808.0/23375.0, -24206.0/37125.0, 338.0/459.0, 0.0, 0.0 },
+    { 5561.0/2376.0, -35.0/11.0, -24117.0/31603.0, 899983.0/200772.0, -5225.0/1836.0, 3925.0/4056.0, 0.0 },
+    { 465467.0/266112.0, -2945.0/1232.0, -5610201.0/14158144.0, 10513573.0/3212352.0, -424325.0/205632.0, 376225.0/454272.0, 0.0 }
+  };
+
+  static const CCTK_REAL gamma_array[8] = 
+    { 61.0/864.0, 
+      0.0, 
+      98415.0/321776.0, 
+      16807.0/146016.0, 
+      1375.0/7344.0, 
+      1375.0/5408.0,
+      -37.0/1120.0,
+      1.0/10.0
+    };
+
+  static const CCTK_REAL gammastar_array[8] = 
+    { 821.0/10800.0,
+      0.0, 
+      19683.0/71825.0, 
+      175273.0/912600.0, 
+      395.0/3672.0, 
+      785.0/2704.0,
+      3.0/50.0,
+      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++)
+  {
+    const CCTK_REAL tmp = (*Original_Delta_Time) / cctkGH->cctk_timefac;
+
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL const *)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] = tmp * RHSVar[index];
+    }
+  }
+
+
+  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  {
+    OldVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 1, 
+						 EvolvedVariableIndex[var]);
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL const *)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);
+
+    if (*MoL_Intermediate_Step - 1)
+    { 
+
+      for (index = 0; index < totalsize; index++)
+      {
+        UpdateVar[index] = OldVar[index];
+      }
+
+      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 (index = 0; index < totalsize; index++)
+      {
+        UpdateVar[index] = OldVar[index];
+        ErrorVar[index] = 0;
+      }
+    
+      for (scratchstep = 0; scratchstep < 8; 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, "Peter has been too lazy to write the RK65 routine "
+              "out for array variables. Better send him an email...");
+
+  }
+        
+  return;
+}
+
+/********************************************************************
+ *********************     Local Routines   *************************
+ ********************************************************************/
+