Add the Classic RK3 method (as a generic method, so use Generic_Type).

Agrees with other RK3's to floating point round off (except at
boundaries) for linear case. Uses more storage and is slower than
standard RK3 so I don't recommend it.

This showed up (so I fixed) a bug with the generic methods when used
with Carpet.


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

6 files changed, 108 insertions, 33 deletions

diff --git a/param.ccl b/param.ccl
index 74752ad..e982acb 100644
--- a/param.ccl
+++ b/param.ccl
@@ -94,6 +94,7 @@ KEYWORD Generic_Type "If using the generic method, which sort"
   "RK"		:: "One of the standard TVD Runga-Kutta methods"
   "ICN"		:: "Iterative Crank Nicholson as a generic method"
   "Table"       :: "Given from the generic method descriptor parameter"
+  "Classic RK3"	:: "Efficient RK3 - classical version"
 } "RK"
 
 CCTK_INT MoL_Intermediate_Steps "Number of intermediate steps taken by the ODE method"
diff --git a/src/GenericRK.c b/src/GenericRK.c
index 6c052af..942866b 100644
--- a/src/GenericRK.c
+++ b/src/GenericRK.c
@@ -108,7 +108,7 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
 
 #ifdef MOLDOESCOMPLEX
   
-  Complex_Delta_Time = CCTK_Cmplx((*Original_Delta_Time), 0);
+  Complex_Delta_Time = CCTK_Cmplx((*Original_Delta_Time) / cctkGH->cctk_timefac, 0);
   Complex_beta = CCTK_Cmplx(beta, 0);
 
 #endif
@@ -122,22 +122,26 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
                                               EvolvedVariableIndex[var]);
     RHSVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
                                            RHSVariableIndex[var]);
-/* #define MOLDEBUG */
+/* #define MOLDEBUG 1 */
 #ifdef MOLDEBUG
-    printf("In generic RK. Variable %d (%s). RHS %d (%s).\n",
+    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]));
+           CCTK_VarName(RHSVariableIndex[var]),
+           beta);
 #endif
 
     for (index = 0; index < totalsize; index++)
     {
-      UpdateVar[index] = (*Original_Delta_Time) * beta * RHSVar[index];
+      UpdateVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * beta * RHSVar[index];
 #ifdef MOLDEBUG
-      printf("Variable: %d. Index: %d. dt: %f. beta %f. RHS: %f. q: %f.\n",
-             var, index, (*Original_Delta_Time), beta, RHSVar[index], 
-             UpdateVar[index]);
+      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
     }
     
@@ -150,6 +154,15 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
       scratchindex = scratchstep - 1;
 
       alpha = RKAlphaCoefficients[alphaindex];
+#ifdef MOLDEBUG
+      printf("In generic RK. Variable %d (%s). RHS %d (%s). step %d. alpha %g.\n",
+             EvolvedVariableIndex[var],
+             CCTK_VarName(EvolvedVariableIndex[var]),
+             RHSVariableIndex[var],
+             CCTK_VarName(RHSVariableIndex[var]),
+             scratchstep,
+             alpha);
+#endif
 
       if (scratchstep) 
       {
@@ -166,9 +179,12 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
                                       + var 
                                       + MoL_Num_Evolved_Vars * scratchindex);
 #ifdef MOLDEBUG
-        printf("Reading from scratch space, initial address %ld index %d\n", 
-               ScratchVar, (var * MoL_Num_Scratch_Levels + 
-                            scratchindex) * totalsize);
+        if (CCTK_EQUALS(verbose,"extreme"))
+        {
+          printf("Reading from scratch space, initial address %ld index %d\n", 
+                 ScratchVar, (var * MoL_Num_Scratch_Levels + 
+                              scratchindex) * totalsize);
+        }
 #endif
       }
       else
@@ -183,10 +199,13 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
         {
           UpdateVar[index] += alpha * ScratchVar[index];
 #ifdef MOLDEBUG
-          printf("Variable: %d. Index: %d. step: %d. "
-                 "alpha: %f. Scratch: %f. q: %f.\n",
-                 var, index, (*MoL_Intermediate_Step), alpha, 
-                 ScratchVar[index], UpdateVar[index]);
+          if (CCTK_EQUALS(verbose,"extreme"))
+          {
+            printf("Variable: %d. Index: %d. step: %d. "
+                   "alpha: %f. Scratch: %f. q: %f.\n",
+                   var, index, (*MoL_Intermediate_Step), alpha, 
+                   ScratchVar[index], UpdateVar[index]);
+          }
 #endif
         }
       }
@@ -225,9 +244,12 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
       {
         ScratchVar[index] = UpdateVar[index];
 #ifdef MOLDEBUG
-        printf("Variable: %d. Index: %d. step: %d. Scratch: %f.\n",
-               var, index, (*MoL_Intermediate_Step), ScratchVar[index]);
-        fflush(stdout);
+        if (CCTK_EQUALS(verbose,"extreme"))
+        {
+          printf("Variable: %d. Index: %d. step: %d. Scratch: %f.\n",
+                 var, index, (*MoL_Intermediate_Step), ScratchVar[index]);
+          fflush(stdout);
+        }
 #endif        
       }
     }
@@ -344,7 +366,7 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
 
     for (index = 0; index < arraytotalsize; index++)
     {
-      UpdateVar[index] = (*Original_Delta_Time) * beta * RHSVar[index];
+      UpdateVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * beta * RHSVar[index];
     }
     
     for (scratchstep = 0; 
@@ -363,9 +385,12 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
                                        (MoL_Max_Evolved_Array_Size+1) + 
                                         arrayscratchlocation];
 #ifdef MOLDEBUG
-        printf("Reading from scratch space, initial address %ld index %d\n", 
-               ScratchVar, scratchindex*(MoL_Max_Evolved_Array_Size+1) + 
-               arrayscratchlocation);
+        if (CCTK_EQUALS(verbose,"extreme"))
+        {
+          printf("Reading from scratch space, initial address %ld index %d\n", 
+                 ScratchVar, scratchindex*(MoL_Max_Evolved_Array_Size+1) + 
+                 arrayscratchlocation);
+        }
 #endif
       }
       else
@@ -380,10 +405,13 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
         {
           UpdateVar[index] += alpha * ScratchVar[index];
 #ifdef MOLDEBUG
-          printf("Variable: %d. Index: %d. step: %d. "
-                 "alpha: %f. Scratch: %f. q: %f.\n",
-                 var, index, (*MoL_Intermediate_Step), 
-                 alpha, ScratchVar[index], UpdateVar[index]);
+          if (CCTK_EQUALS(verbose,"extreme"))
+          {
+            printf("Variable: %d. Index: %d. step: %d. "
+                   "alpha: %f. Scratch: %f. q: %f.\n",
+                   var, index, (*MoL_Intermediate_Step), 
+                   alpha, ScratchVar[index], UpdateVar[index]);
+          }
 #endif
         }
       }
@@ -433,8 +461,11 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
       {
         ScratchVar[index] = UpdateVar[index];
 #ifdef MOLDEBUG
-        printf("Variable: %d. Index: %d. step: %d. Scratch: %f.\n",
-               var, index, (*MoL_Intermediate_Step), ScratchVar[index]);
+        if (CCTK_EQUALS(verbose,"extreme"))
+        {
+          printf("Variable: %d. Index: %d. step: %d. Scratch: %f.\n",
+                 var, index, (*MoL_Intermediate_Step), ScratchVar[index]);
+        }
 #endif        
       }
       arrayscratchlocation += arraytotalsize;
diff --git a/src/IndexArrays.c b/src/IndexArrays.c
index e0df297..096a745 100644
--- a/src/IndexArrays.c
+++ b/src/IndexArrays.c
@@ -282,6 +282,10 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
       sprintf(infoline, "Generic method, options:\n %s\n",
               Generic_Method_Descriptor);
     }
+    else if (CCTK_EQUALS(Generic_Type,"Classic RK3")) 
+    {
+      sprintf(infoline, "Classic Runge-Kutta 3");
+    }
     else
     {
       CCTK_WARN(0, "Generic_Type not recognized!");
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index 47c573a..5316c9f 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -75,6 +75,13 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
                      "of scratch levels must be at least the "
                      "number of intermediate steps - 1");
     }
+    if ( (CCTK_Equals(Generic_Type, "Classic RK3"))&&
+         ((!(MoL_Intermediate_Steps == 3))||(!(MoL_Num_Scratch_Levels > 1))) )
+    {
+      CCTK_PARAMWARN("When using the classic RK3 evolver the "
+                     "number of intermediate steps must be 3 "
+                     "and the number of scratch levels at least 2");
+    }
     if (CCTK_Equals(Generic_Type, "Table"))
     {
       options_table =
@@ -121,6 +128,7 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
                    "number of intermediate steps must be 3");
   }
 
+
   return 0;
   
 }
diff --git a/src/RKCoefficients.c b/src/RKCoefficients.c
index 905fb60..57ac38a 100644
--- a/src/RKCoefficients.c
+++ b/src/RKCoefficients.c
@@ -73,7 +73,36 @@ int MoL_SetupRKCoefficients(CCTK_ARGUMENTS)
 
   CCTK_INT ierr, options_table;
 
-  if (CCTK_Equals(Generic_Type,"ICN")) 
+  if (CCTK_EQUALS(Generic_Type,"Classic RK3"))
+  {
+    if (MoL_Num_Scratch_Levels != 2)
+    {
+      CCTK_WARN(0, "For Classic RK3, MoL_Num_Scratch_Levels "
+                "should be at least 2");
+    }
+    if (MoL_Intermediate_Steps != 3)
+    {
+      CCTK_WARN(0, "For Classic RK3, MoL_Intermediate_Steps "
+                "should be at least 3");
+    }
+    for (i = 0; i < MoL_Intermediate_Steps; i++) 
+    {
+      for (j = 0; j < MoL_Num_Scratch_Levels + 1; j++) 
+      {  
+        RKAlphaCoefficients[i * MoL_Intermediate_Steps + j] = 0.0;
+      }
+      RKBetaCoefficients[i] = 0.0;
+    }
+    RKAlphaCoefficients[0] = 1.0;
+    RKAlphaCoefficients[3] = 1.0;
+    RKAlphaCoefficients[6] = 1.0 / 9.0;
+    RKAlphaCoefficients[7] = 4.0 / 9.0;
+    RKAlphaCoefficients[8] = 4.0 / 9.0;
+    RKBetaCoefficients[0] = 0.5;
+    RKBetaCoefficients[1] = 0.75;
+    RKBetaCoefficients[2] = 4.0 / 9.0;    
+  }
+  else if (CCTK_EQUALS(Generic_Type,"ICN")) 
   {
     for (i = 0; i < MoL_Intermediate_Steps; i++) 
     {
@@ -92,7 +121,7 @@ int MoL_SetupRKCoefficients(CCTK_ARGUMENTS)
       } 
     }
   }
-  else if (CCTK_Equals(Generic_Type,"RK")) 
+  else if (CCTK_EQUALS(Generic_Type,"RK")) 
   {
     if (MoL_Num_Scratch_Levels < MoL_Intermediate_Steps - 1)
     {
@@ -152,7 +181,7 @@ int MoL_SetupRKCoefficients(CCTK_ARGUMENTS)
               "with MoL_Intermediate_Steps greater than 4");
     }
   }
-  else if (CCTK_Equals(Generic_Type,"Table")) 
+  else if (CCTK_EQUALS(Generic_Type,"Table")) 
   {
     if (MoL_Num_Scratch_Levels < MoL_Intermediate_Steps - 1)
     {
diff --git a/src/SetTime.c b/src/SetTime.c
index 76a06a4..94380b1 100644
--- a/src/SetTime.c
+++ b/src/SetTime.c
@@ -19,6 +19,8 @@ static const char *rcsid = "$Header$";
 
 CCTK_FILEVERSION(CactusBase_MoL_SetTime_c);
 
+/* #define MOLDEBUG 1 */
+
 /********************************************************************
  *********************     Local Data Types   ***********************
  ********************************************************************/
@@ -192,7 +194,7 @@ int MoL_ResetTime(CCTK_ARGUMENTS)
     }
     else if (*MoL_Intermediate_Step == 1)
     {
-      cctkGH->cctk_time = (*Original_Time) + 
+      cctkGH->cctk_time = (*Original_Time) - 
         0.5*(*Original_Delta_Time)/cctkGH->cctk_timefac;
     }
   }
@@ -254,7 +256,7 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS)
   {
     cctkGH->cctk_delta_time = RKBetaCoefficients[MoL_Intermediate_Steps - 
                                                 (*MoL_Intermediate_Step)] *
-      (*Original_Delta_Time)/cctkGH->cctk_timefac;
+      (*Original_Delta_Time);
   }
   else if (CCTK_EQUALS(ODE_Method,"RK2"))
   {