1 files changed, 43 insertions, 68 deletions

diff --git a/src/RK4.c b/src/RK4.c
index 1d790ab..4e72389 100644
--- a/src/RK4.c
+++ b/src/RK4.c
@@ -13,8 +13,8 @@
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
-#include <stdio.h>
 #include "ExternalVariables.h"
+#include "Operators.h"
 
 static const char *rcsid = "$Header$";
 
@@ -99,14 +99,14 @@ CCTK_WARN(0, "not implemented");
   /* Keep a running total of alpha as we perform the substeps, so that
      we know the "real" alpha (including round-off errors) when we
      calculate the final result. */
-  static CCTK_REAL sum_alpha;
-  /* the last MoL intermediate step that was summed */
-  static CCTK_INT last_sum_step;
-
+  CCTK_REAL const time_rhs = 1.0;
+  CCTK_REAL const old_time = 0.0;
+  static CCTK_REAL time, scratch_time;
+  
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }
 
   if (scratchspace_firstindex == -99)
@@ -131,21 +131,14 @@ CCTK_WARN(0, "not implemented");
     case 3:
       alpha = 1.0;
       beta  = 1.0 / 6.0;
+      break;
   }
 
-  /* Initialise alpha before the first intermediate step */
   if (MoL_Intermediate_Steps == (*MoL_Intermediate_Step))
   {
-    sum_alpha = 0.0;
-    last_sum_step = -1;
-  }
-  /* Add alpha once per intermediate step */
-  if (last_sum_step != (*MoL_Intermediate_Step))
-  {
-    sum_alpha += alpha;
-    last_sum_step = (*MoL_Intermediate_Step);
+    time = 0.0;
   }
-
+  
   /* FIXME */
 
 #ifdef MOLDOESCOMPLEX
@@ -159,71 +152,53 @@ CCTK_WARN(0, "not implemented");
 
   for (var = 0; var < MoLNumEvolvedVariables; var++)
   {
-    
-    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
-                                              EvolvedVariableIndex[var]);
-    OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, 
-                                              EvolvedVariableIndex[var]);
-    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                 RHSVariableIndex[var]);
-/* #define MOLDEBUG 1 */
-#ifdef MOLDEBUG
-    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]),
-           beta);
-#endif
-
-#pragma omp parallel for
-    for (index = 0; index < totalsize; index++)
     {
-      UpdateVar[index] = OldVar[index] + 
-        (*Original_Delta_Time) / cctkGH->cctk_timefac * beta * RHSVar[index];
-#ifdef MOLDEBUG
-      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
+      int const nsrcs = 2;
+      CCTK_INT const srcs[] =
+        {EvolvedVariableIndex[var], RHSVariableIndex[var]};
+      CCTK_INT const tls[] = {1, 0};
+      CCTK_REAL const facts[] =
+        {1.0, (*Original_Delta_Time) / cctkGH->cctk_timefac * beta};
+      MoL_LinearCombination(cctkGH,
+                            EvolvedVariableIndex[var], 0.0,
+                            srcs, tls, facts, nsrcs);
+      time = facts[0] * old_time + facts[1] * time_rhs;
     }
-    ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
-                                      scratchspace_firstindex
-                                      + var );
-    
+
     /* scratch storage */
     if ((*MoL_Intermediate_Step) == MoL_Intermediate_Steps)
     {
-#pragma omp parallel for
-      for (index = 0; index < totalsize; index++)
-      {
-        ScratchVar[index] = 0;
-      }
+      MoL_LinearCombination(cctkGH,
+                            scratchspace_firstindex + var, 0.0,
+                            NULL, NULL, NULL, 0);
+      scratch_time = 0.0;
     }
 
     if ((*MoL_Intermediate_Step)>1)
     {
-#pragma omp parallel for
-      for (index = 0; index < totalsize; index++)
-      {
-        ScratchVar[index] += alpha * UpdateVar[index];
-      }
+      int const nsrcs = 1;
+      CCTK_INT const srcs[] = {EvolvedVariableIndex[var]};
+      CCTK_INT const tls[] = {0};
+      CCTK_REAL const facts[] = {alpha};
+      MoL_LinearCombination(cctkGH,
+                            scratchspace_firstindex + var, 1.0,
+                            srcs, tls, facts, nsrcs);
+      scratch_time += facts[0] * time;
     }
     else
     {
-#pragma omp parallel for
-      for (index = 0; index < totalsize; index++)
-      {
-        UpdateVar[index] += ScratchVar[index] - (sum_alpha - 1.0) * OldVar[index];
-      }
+      int const nsrcs = 2;
+      CCTK_INT const srcs[] =
+        {scratchspace_firstindex + var, EvolvedVariableIndex[var]};
+      CCTK_INT const tls[] = {0, 1};
+      CCTK_REAL const facts[] = {1.0, -4.0/3.0};
+      MoL_LinearCombination(cctkGH,
+                            EvolvedVariableIndex[var], 1.0,
+                            srcs, tls, facts, nsrcs);
+      time += facts[0] * scratch_time + facts[1] * old_time;
     }
-
   }
 
-
   /* Complex GFs */
 
   /* FIXME */
@@ -259,7 +234,7 @@ CCTK_WARN(0, "not done");
     arraytotalsize = 1;
     for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
     {
-      arraytotalsize *= arraydata.lsh[arraydim];
+      arraytotalsize *= arraydata.ash[arraydim];
     }
 
     ScratchVar = &ArrayScratchSpace[arrayscratchlocation];
@@ -293,7 +268,7 @@ CCTK_WARN(0, "not done");
 #pragma omp parallel for
       for (index = 0; index < arraytotalsize; index++)
       {
-        UpdateVar[index] += ScratchVar[index] - (sum_alpha - 1.0) * OldVar[index];
+        UpdateVar[index] += ScratchVar[index] - 4.0 / 3.0 * OldVar[index];
       }
     }
     arrayscratchlocation += arraytotalsize;