MoL Update

New integrator Euler. This is an explicit, first-order method, mostly
useful for debugging.

New integrators AB (Adams-Bashforth) with various orders. These are
explicit integrators using several past timelevels to provide
higher-order integration with a single RHS evaluation each.

Introduce LinearCombination, a generic routine to calculate linear
combinations of grid functions. This simplifies existing code, and can
be overloaded if MoL should run on a device (e.g. with OpenCL).

Replace cctk_lsh with cctk_ash where necessary.


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

24 files changed, 1073 insertions, 520 deletions

diff --git a/interface.ccl b/interface.ccl
index d2d14f9..94a4bb7 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -206,10 +206,28 @@ PROVIDES FUNCTION MoLNumIntegratorSubsteps WITH MoL_NumIntegratorSubsteps \
 # PROVIDES FUNCTION MoLChangeToNoneComplexArray WITH \
 #                   MoL_ChangeToNoneComplexArray LANGUAGE C
 
-##################################################
-### Functions provided by devices              ###
-##################################################
-
+##############################################################
+### Aliased functions that can be provided by other thorns ###
+### to override low-level grid variable operations.        ###
+##############################################################
+
+CCTK_INT FUNCTION LinearCombination             \
+    (CCTK_POINTER_TO_CONST IN cctkGH,           \
+     CCTK_INT              IN var,              \
+     CCTK_REAL             IN scale,            \
+     CCTK_INT  ARRAY       IN srcs,             \
+     CCTK_INT  ARRAY       IN tls,              \
+     CCTK_REAL ARRAY       IN facts,            \
+     CCTK_INT              IN nsrcs)
+USES FUNCTION LinearCombination
+
+CCTK_INT FUNCTION Accelerator_NotifyDataModified \
+  (CCTK_POINTER_TO_CONST IN cctkGH,              \
+   CCTK_INT ARRAY        IN variables,           \
+   CCTK_INT ARRAY        IN tls,                 \
+   CCTK_INT              IN nvariables,          \
+   CCTK_INT              IN on_device)
+USES FUNCTION Accelerator_NotifyDataModified
 
 private:
 
diff --git a/param.ccl b/param.ccl
index 43a1fd4..cce525f 100644
--- a/param.ccl
+++ b/param.ccl
@@ -91,6 +91,7 @@ KEYWORD ODE_Method "The ODE method use by MoL to do time integration"
   "Generic"	:: "Generic Shu-Osher Runga-Kutta type"
   "ICN"		:: "Iterative Crank Nicholson"
   "ICN-avg"	:: "Iterative Crank Nicholson with averaging"
+  "Euler"	:: "Euler"
   "RK2"		:: "Efficient RK2"
   "RK3"		:: "Efficient RK3"
   "RK4"		:: "Efficient RK4"
@@ -98,6 +99,7 @@ KEYWORD ODE_Method "The ODE method use by MoL to do time integration"
   "RK45CK"      :: "RK45CK (Cash-Karp) with error estimation"
   "RK65"        :: "RK65 with error estimation"
   "RK87"        :: "RK87 with error estimation"
+  "AB"          :: "Adams-Bashforth"
   "RK2-MR-2:1"  :: "2nd order 2:1 multirate RK scheme based on RK2 due to Schlegel et al 2009"
   "RK4-MR-2:1"  :: "3rd order 2:1 multirate RK scheme based on RK43 due to Schlegel et al 2009"
   "RK4-RK2"     :: "RK4 as fast method and RK2 as slow method"
@@ -120,6 +122,19 @@ BOOLEAN ICN_avg_swapped "Use swapped averages in ICN method?"
 {
 } "no"
 
+KEYWORD AB_Type "If using the the AB method, which sort"
+{
+  "1" :: "same as forward Euler"
+  "2" :: "second order"
+  "3" :: "third order"
+  "4" :: "fourth order"
+  "5" :: "fifth order"
+} "1"
+
+BOOLEAN AB_initially_reduce_order "Reduce order of accuracy initially so that no past timelevels of initial data are required"
+{
+} "yes"
+
 CCTK_INT MoL_Intermediate_Steps "Number of intermediate steps taken by the ODE method"
 {
   1:*		:: "Anything greater than 1"
@@ -152,6 +167,11 @@ BOOLEAN copy_ID_after_MoL_PostStep \
 {
 } "no"
 
+BOOLEAN set_ID_boundaries "Should boundaries be overwritten (via synchronization, prolongation, boundary conditions) by MoL?"
+{
+} "yes"
+
+
 
 # The default for this parameter corresponds to generic RK2
 STRING Generic_Method_Descriptor "A string used to create a table containing the description of the generic method"
@@ -221,6 +241,8 @@ KEYWORD verbose "How verbose should MoL be?"
   "extreme"  :: "Everything you never wanted to know"
 } "normal"
 
-BOOLEAN set_ID_boundaries "Should boundaries be overwritten (via synchronization, prolongation, boundary conditions) by MoL?"
-{
-} "yes"
+
+
+shares: Cactus
+
+USES CCTK_REAL cctk_initial_time
diff --git a/schedule.ccl b/schedule.ccl
index 48283fd..21896e3 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -551,6 +551,13 @@ if (CCTK_Equals(ODE_Method,"Generic"))
     LANG: C
   } "Updates calculated with a generic method"
 }
+else if (CCTK_Equals(ODE_Method,"Euler")) 
+{
+  schedule MoL_EulerAdd AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
+  {
+    LANG: C
+  } "Updates calculated with the Euler method"
+}
 else if (CCTK_Equals(ODE_Method,"RK2")) 
 {
   schedule MoL_RK2Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
@@ -613,6 +620,13 @@ else if (CCTK_Equals(ODE_Method,"ICN-avg"))
     LANG: C
   } "Updates calculated with the averaging ICN method"
 }
+else if (CCTK_Equals(ODE_Method,"AB"))
+{
+  schedule MoL_ABAdd AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
+  {
+    LANG: C
+  } "Updates calculated with the Adams-Bashforth"
+}
 else if (CCTK_Equals(ODE_Method,"RK2-MR-2:1")) 
 {
   schedule MoL_RK2_MR_2_1_Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
@@ -765,7 +779,7 @@ schedule MoL_ResetDeltaTime IN MoL_Step AFTER (MoL_PostStep MoL_PostStepModify)
 ###  variables to their original state.        ###
 ##################################################
 
-schedule MoL_RestoreSandR IN MoL_Evolution AFTER MoL_FinishLoop
+schedule MoL_RestoreSandR IN MoL_Evolution AFTER (MoL_ReduceAdaptiveError MoL_FinishLoop)
 {
   LANG: C
 } "Restoring the Save and Restore variables to the original state"
diff --git a/src/AB.c b/src/AB.c
new file mode 100644
index 0000000..ee7fdad
--- /dev/null
+++ b/src/AB.c
@@ -0,0 +1,234 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include <assert.h>
+#include <stdlib.h>
+
+#include "ExternalVariables.h"
+
+
+
+/* Coefficients taken from
+ * http://en.wikipedia.org/wiki/Linear_multistep_method, which cites
+ * (Hairer, Nørsett & Wanner 1993, §III.1; Butcher 2003, p. 103).
+ * 
+ * Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (1993),
+ * Solving ordinary differential equations I: Nonstiff problems (2nd
+ * ed.), Berlin: Springer Verlag, ISBN 978-3-540-56670-0.
+ *
+ * Butcher, John C. (2003), Numerical Methods for Ordinary
+ * Differential Equations, John Wiley, ISBN 978-0-471-96758-3.
+ */
+
+/* The following Mathematic expression (see
+ * http://en.wikipedia.org/wiki/Linear_multistep_method) calculates
+ * Adams-Bashforth coefficients for arbitrary orders:
+ *
+ * b = Table[
+ *   Table[(-1)^j/(j! (s - j - 1)!) Integrate[
+ *      Product[If[j == i, 1, u + i], {i, 0, s - 1}], {u, 0, 1}], {j, 0, 
+ *     s - 1}], {s, 1, 8}]
+ *
+ * where s is the desired order.  The coefficients up to order 8 are:
+ * 
+ * {1}
+ *
+ * {3/2, -(1/2)}
+ *
+ * {23/12, -(4/3), 5/12}
+ *
+ * {55/24, -(59/24), 37/24, -(3/8)}
+ *
+ * {1901/720, -(1387/360), 109/30, -(637/360), 251/720}
+ *
+ * {4277/1440, -(2641/480), 4991/720, -(3649/720), 959/480, -(95/288)}
+ *
+ * {198721/60480, -(18637/2520), 235183/20160, -(10754/945),
+ * 135713/20160, -(5603/2520), 19087/60480}
+ *
+ * {16083/4480, -(1152169/120960), 242653/13440, -(296053/13440),
+ * 2102243/120960, -(115747/13440), 32863/13440, -(5257/17280)}
+ */
+
+
+
+static void order1 (CCTK_REAL* restrict const UpdateVar,
+                    CCTK_REAL const* restrict *restrict const RHSVars,
+                    CCTK_REAL const dt,
+                    int const totalsize)
+{
+  CCTK_REAL const* restrict const RHSVar0 = RHSVars[0];
+#pragma omp parallel for
+  for (int index = 0; index < totalsize; index++) {
+    UpdateVar[index] += dt * RHSVar0[index];
+  }
+}
+
+static void order2 (CCTK_REAL* restrict const UpdateVar,
+                    CCTK_REAL const* restrict *restrict const RHSVars,
+                    CCTK_REAL const dt,
+                    int const totalsize)
+{
+  CCTK_REAL const* restrict const RHSVar0 = RHSVars[0];
+  CCTK_REAL const* restrict const RHSVar1 = RHSVars[1];
+  CCTK_REAL const f0 = + (3.0/2.0) * dt;
+  CCTK_REAL const f1 = - (1.0/2.0) * dt;
+#pragma omp parallel for
+  for (int index = 0; index < totalsize; index++) {
+    UpdateVar[index] += f0 * RHSVar0[index] + f1 * RHSVar1[index];
+  }
+}
+
+static void order3 (CCTK_REAL* restrict const UpdateVar,
+                    CCTK_REAL const* restrict *restrict const RHSVars,
+                    CCTK_REAL const dt,
+                    int const totalsize)
+{
+  CCTK_REAL const* restrict const RHSVar0 = RHSVars[0];
+  CCTK_REAL const* restrict const RHSVar1 = RHSVars[1];
+  CCTK_REAL const* restrict const RHSVar2 = RHSVars[2];
+  CCTK_REAL const f0 = + (23.0/12.0) * dt;
+  CCTK_REAL const f1 = - ( 4.0/ 3.0) * dt;
+  CCTK_REAL const f2 = + ( 5.0/12.0) * dt;
+#pragma omp parallel for
+  for (int index = 0; index < totalsize; index++) {
+    UpdateVar[index] +=
+      f0 * RHSVar0[index] + f1 * RHSVar1[index] + f2 * RHSVar2[index];
+  }
+}
+
+static void order4 (CCTK_REAL* restrict const UpdateVar,
+                    CCTK_REAL const* restrict *restrict const RHSVars,
+                    CCTK_REAL const dt,
+                    int const totalsize)
+{
+  CCTK_REAL const* restrict const RHSVar0 = RHSVars[0];
+  CCTK_REAL const* restrict const RHSVar1 = RHSVars[1];
+  CCTK_REAL const* restrict const RHSVar2 = RHSVars[2];
+  CCTK_REAL const* restrict const RHSVar3 = RHSVars[3];
+  CCTK_REAL const f0 = + (55.0/24.0) * dt;
+  CCTK_REAL const f1 = - (59.0/24.0) * dt;
+  CCTK_REAL const f2 = + (37.0/24.0) * dt;
+  CCTK_REAL const f3 = - ( 3.0/ 8.0) * dt;
+#pragma omp parallel for
+  for (int index = 0; index < totalsize; index++) {
+    UpdateVar[index] +=
+      f0 * RHSVar0[index] + f1 * RHSVar1[index] + f2 * RHSVar2[index] +
+      f3 * RHSVar3[index];
+  }
+}
+
+static void order5 (CCTK_REAL* restrict const UpdateVar,
+                    CCTK_REAL const* restrict *restrict const RHSVars,
+                    CCTK_REAL const dt,
+                    int const totalsize)
+{
+  CCTK_REAL const* restrict const RHSVar0 = RHSVars[0];
+  CCTK_REAL const* restrict const RHSVar1 = RHSVars[1];
+  CCTK_REAL const* restrict const RHSVar2 = RHSVars[2];
+  CCTK_REAL const* restrict const RHSVar3 = RHSVars[3];
+  CCTK_REAL const* restrict const RHSVar4 = RHSVars[4];
+  CCTK_REAL const f0 = + (1901.0/720.0) * dt;
+  CCTK_REAL const f1 = - (1387.0/360.0) * dt;
+  CCTK_REAL const f2 = + ( 109.0/ 30.0) * dt;
+  CCTK_REAL const f3 = - ( 637.0/360.0) * dt;
+  CCTK_REAL const f4 = + ( 251.0/720.0) * dt;
+#pragma omp parallel for
+  for (int index = 0; index < totalsize; index++) {
+    UpdateVar[index] +=
+      f0 * RHSVar0[index] + f1 * RHSVar1[index] + f2 * RHSVar2[index] +
+      f3 * RHSVar3[index] + f4 * RHSVar4[index];
+  }
+}
+
+/* Array of function pointers */
+static
+void (* const orders[]) (CCTK_REAL* restrict const UpdateVar,
+                         CCTK_REAL const* restrict *restrict const RHSVars,
+                         CCTK_REAL const dt,
+                         int const totalsize)
+= { order1, order2, order3, order4, order5 };
+static int const max_order = sizeof orders / sizeof *orders;
+
+
+
+void MoL_ABAdd(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_REAL const dt = CCTK_DELTA_TIME;
+  
+  /* Determine the order of accuracy */
+  int order;
+  if (CCTK_EQUALS(AB_Type,"1")) {
+    order = 1 ;
+  } else if (CCTK_EQUALS(AB_Type,"2")) {
+    order = 2;
+  } else if (CCTK_EQUALS(AB_Type,"3")) {
+    order = 3;
+  } else if (CCTK_EQUALS(AB_Type,"4")) {
+    order = 4;
+  } else if (CCTK_EQUALS(AB_Type,"5")) {
+    order = 5;
+  } else {
+    abort();
+  }
+  if (AB_initially_reduce_order) {
+    /* Reduce the order for the first time steps */
+    /* int const iteration = cctk_iteration; */
+    int const iteration = 1 + floor((cctk_time - cctk_initial_time) / dt + 0.5);
+    if (order > iteration) {
+      order = iteration;
+      CCTK_VInfo (CCTK_THORNSTRING,
+                  "Reducing Adams-Bashforth order to %d", order);
+    }
+  }
+/* printf ("MoL AB: iter=%d, order=%d\n", cctk_iteration, order); */
+  assert (order >= 1 && order <= max_order);
+  
+  int totalsize = 1;
+  for (int arraydim = 0; arraydim < cctk_dim; arraydim++) {
+    totalsize *= cctk_ash[arraydim];
+  }
+  
+  for (int var = 0; var < MoLNumEvolvedVariables; var++) {
+    CCTK_REAL* restrict const UpdateVar =
+      CCTK_VarDataPtrI(cctkGH, 0, EvolvedVariableIndex[var]);
+    CCTK_REAL const* restrict RHSVars[order];
+    for (int tl = 0; tl < order; tl++) {
+      RHSVars[tl] = CCTK_VarDataPtrI(cctkGH, tl, RHSVariableIndex[var]);
+    }
+    
+    /* Add RHS */
+    (orders[order-1]) (UpdateVar, RHSVars, dt, totalsize);
+  } /* var */
+  
+  for (int var = 0; var < MoLNumEvolvedArrayVariables; var++) {
+    CCTK_REAL* restrict const UpdateVar =
+      CCTK_VarDataPtrI(cctkGH, 0, EvolvedArrayVariableIndex[var]);
+    CCTK_REAL const* restrict RHSVars[order];
+    for (int tl = 0; tl < order; tl++) {
+      RHSVars[tl] = CCTK_VarDataPtrI(cctkGH, tl, RHSArrayVariableIndex[var]);
+    }
+    
+    int const groupindex =
+      CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
+    cGroupDynamicData arraydata;
+    int const ierr = CCTK_GroupDynamicData(cctkGH, groupindex, &arraydata);
+    if (ierr) {
+      CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, 
+                 "The driver does not return group information for group '%s'.",
+                 CCTK_GroupName(groupindex));
+    }
+    int arraytotalsize = 1;
+    for (int arraydim = 0; arraydim < arraydata.dim; arraydim++) {
+      arraytotalsize *= arraydata.ash[arraydim];
+    }
+    
+    /* Add RHS */
+    (orders[order-1]) (UpdateVar, RHSVars, dt, arraytotalsize);
+  } /* var */
+
+}
diff --git a/src/Euler.c b/src/Euler.c
new file mode 100644
index 0000000..0033cd1
--- /dev/null
+++ b/src/Euler.c
@@ -0,0 +1,116 @@
+ /*@@
+   @file      Euler.c
+   @date      2012-12-05
+   @author    Erik Schnetter
+   @desc 
+   An explicit Euler time integrator. This method is unstable in most
+   cases, but is sometimes useful for debugging.
+   @enddesc 
+   @version   $Header$
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "ExternalVariables.h"
+#include "Operators.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusBase_MoL_Euler_c);
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_EulerAdd(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    MoL_EulerAdd
+   @date       2012-12-05
+   @author     Erik Schnetter
+   @desc 
+   Performs first order Euler time integration.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_EulerAdd(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  int totalsize = 1;
+  for (int arraydim = 0; arraydim < cctk_dim; arraydim++)
+  {
+    totalsize *= cctk_ash[arraydim];
+  }
+  
+  switch (*MoL_Intermediate_Step)
+  {
+  case 1:
+    {
+      for (int var = 0; var < MoLNumEvolvedVariables; var++)
+      {
+        int const nsrcs = 2;
+        CCTK_INT const srcs[] =
+          {EvolvedVariableIndex[var], RHSVariableIndex[var]};
+        CCTK_INT const tls[] = {1, 0};
+        CCTK_REAL const facts[] = {1.0, CCTK_DELTA_TIME};
+        MoL_LinearCombination(cctkGH,
+                              EvolvedVariableIndex[var], 0.0,
+                              srcs, tls, facts, nsrcs);
+      }
+      
+      for (int var = 0; var < MoLNumEvolvedArrayVariables; var++)
+      {
+        int const nsrcs = 2;
+        CCTK_INT const srcs[] =
+          {EvolvedArrayVariableIndex[var], RHSArrayVariableIndex[var]};
+        CCTK_INT const tls[] = {1, 0};
+        CCTK_REAL const facts[] = {1.0, CCTK_DELTA_TIME};
+        MoL_LinearCombination(cctkGH,
+                              EvolvedArrayVariableIndex[var], 0.0,
+                              srcs, tls, facts, nsrcs);
+      }
+      
+      break;
+    }
+  default:
+    {
+      CCTK_WARN(CCTK_WARN_ABORT,
+                "Euler expects MoL_Intermediate_Step to be in [1,1]. "
+                "This should be caught at ParamCheck - bug Ian!");
+    }
+  }
+}
+
+/********************************************************************
+ *********************     Local Routines   *************************
+ ********************************************************************/
diff --git a/src/GenericRK.c b/src/GenericRK.c
index 4a7b854..6c88ca6 100644
--- a/src/GenericRK.c
+++ b/src/GenericRK.c
@@ -99,7 +99,7 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }
 
   if (scratchspace_firstindex == -99)
diff --git a/src/ICN.c b/src/ICN.c
index fb71783..af074e3 100644
--- a/src/ICN.c
+++ b/src/ICN.c
@@ -105,7 +105,7 @@ void MoL_ICNAdd(CCTK_ARGUMENTS)
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }
 
 #ifdef MOLDEBUG
@@ -149,7 +149,7 @@ void MoL_ICNAdd(CCTK_ARGUMENTS)
     arraytotalsize = 1;
     for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
     {
-      arraytotalsize *= arraydata.lsh[arraydim];
+      arraytotalsize *= arraydata.ash[arraydim];
     }
 
 /*     CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,  */
@@ -256,7 +256,7 @@ void MoL_ICNAverage(CCTK_ARGUMENTS)
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }
 
 #ifdef MOLDEBUG
@@ -298,7 +298,7 @@ void MoL_ICNAverage(CCTK_ARGUMENTS)
     arraytotalsize = 1;
     for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
     {
-      arraytotalsize *= arraydata.lsh[arraydim];
+      arraytotalsize *= arraydata.ash[arraydim];
     }
 
 #pragma omp parallel for
diff --git a/src/IndexArrays.c b/src/IndexArrays.c
index 485a4bd..65a4dac 100644
--- a/src/IndexArrays.c
+++ b/src/IndexArrays.c
@@ -12,6 +12,7 @@
 #include "cctk.h"
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
+#include "util_String.h"
 
 #include <stdlib.h>
 #include <stdio.h>
@@ -313,6 +314,10 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
       CCTK_WARN(0, "Generic_Type not recognized!");
     }
   }
+  else if (CCTK_EQUALS(ODE_Method,"Euler")) 
+  {
+    sprintf(infoline, "Euler");
+  }
   else if (CCTK_EQUALS(ODE_Method,"RK2")) 
   {
     sprintf(infoline, "Runge-Kutta 2");
@@ -352,6 +357,38 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
             "Averaging iterative Crank Nicholson with %i iterations", 
             MoL_Intermediate_Steps);
   }
+  else if (CCTK_EQUALS(ODE_Method,"AB"))
+  {
+    if (CCTK_EQUALS(AB_Type,"1"))
+    {
+      Util_snprintf(infoline, sizeof infoline,
+                    "Adams-Bashforth of order 1");
+    }
+    else if (CCTK_EQUALS(AB_Type,"2"))
+    {
+      Util_snprintf(infoline, sizeof infoline,
+                    "Adams-Bashforth of order 2");
+    }
+    else if (CCTK_EQUALS(AB_Type,"3"))
+    {
+      Util_snprintf(infoline, sizeof infoline,
+                    "Adams-Bashforth of order 3");
+    }
+    else if (CCTK_EQUALS(AB_Type,"4"))
+    {
+      Util_snprintf(infoline, sizeof infoline,
+                    "Adams-Bashforth of order 4");
+    }
+    else if (CCTK_EQUALS(AB_Type,"5"))
+    {
+      Util_snprintf(infoline, sizeof infoline,
+                    "Adams-Bashforth of order 5");
+    }
+    else
+    {
+      CCTK_WARN(0, "AB_Type not recognized!");
+    }
+  }
   else if (CCTK_EQUALS(ODE_Method,"RK2-MR-2:1")) 
   {
     sprintf(infoline, "Multi-rate 2:1 Runge-Kutta 2");
diff --git a/src/InitialCopy.c b/src/InitialCopy.c
index b421908..5cef13c 100644
--- a/src/InitialCopy.c
+++ b/src/InitialCopy.c
@@ -20,6 +20,7 @@
 #include "cctk_Parameters.h"
 
 #include "ExternalVariables.h"
+#include "Operators.h"
 
 static const char *rcsid = "$Header$";
 
@@ -104,11 +105,15 @@ void MoL_InitialCopy(CCTK_ARGUMENTS)
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }
 
   for (var = 0; var < MoLNumEvolvedVariables; var++)
   {
+    const int       nsrc = 1;
+    const int       srcs[1] = {EvolvedVariableIndex[var]};
+    const int       tls[1] = {1};
+    const CCTK_REAL facts[1] = {1.0};
     
     StorageOn = CCTK_QueryGroupStorageI(cctkGH,
                                        CCTK_GroupIndexFromVarI(EvolvedVariableIndex[var]));
@@ -130,20 +135,8 @@ void MoL_InitialCopy(CCTK_ARGUMENTS)
       CCTK_WARN(0, "The grid function does not have storage assigned.");
     }
     
-    PreviousVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 1, 
-                                               EvolvedVariableIndex[var]);
-    CurrentVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                              EvolvedVariableIndex[var]);
-    if (PreviousVar && CurrentVar)
-    {      
-      memcpy(CurrentVar, PreviousVar, totalsize * sizeof(CCTK_REAL));
-    }
-    else
-    {
-      CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,"Null pointer for variable %s", 
-                 CCTK_VarName(EvolvedVariableIndex[var]));
-    }
-
+    MoL_LinearCombination(cctkGH, EvolvedVariableIndex[var], 0.0,
+                          srcs, tls, facts, nsrc);
   }
     
   /* Set up the array sizes */
@@ -186,7 +179,7 @@ void MoL_InitialCopy(CCTK_ARGUMENTS)
     arraytotalsize = 1;
     for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
     {
-      arraytotalsize *= arraydata.lsh[arraydim];
+      arraytotalsize *= arraydata.ash[arraydim];
     }
 
     ArrayScratchSizes[var] = arraytotalsize;
@@ -526,12 +519,11 @@ void MoL_InitRHS(CCTK_ARGUMENTS)
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }
 
   for (var = 0; var < MoLNumEvolvedVariables; var++)
   {
-    
     StorageOn = CCTK_QueryGroupStorageI(cctkGH,
                                        CCTK_GroupIndexFromVarI(RHSVariableIndex[var]));
     
@@ -552,22 +544,8 @@ void MoL_InitRHS(CCTK_ARGUMENTS)
       CCTK_WARN(0, "The grid function does not have storage assigned.");
     }
     
-    RHSVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                          RHSVariableIndex[var]);
-    if (RHSVar)
-    {
-#pragma omp parallel for
-      for (index = 0; index < totalsize; index++)
-      {
-        RHSVar[index] = 0;
-      }
-    }
-    else
-    {
-      CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,"Null pointer for variable %s", 
-                 CCTK_VarName(RHSVariableIndex[var]));
-    }
-
+    MoL_LinearCombination(cctkGH, RHSVariableIndex[var], 0.0,
+                          NULL, NULL, NULL, 0);
   }
     
   for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
@@ -588,7 +566,7 @@ void MoL_InitRHS(CCTK_ARGUMENTS)
     arraytotalsize = 1;
     for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
     {
-      arraytotalsize *= arraydata.lsh[arraydim];
+      arraytotalsize *= arraydata.ash[arraydim];
     }
 
     if (arraytotalsize)
@@ -690,7 +668,7 @@ void MoL_FillAllLevels(CCTK_ARGUMENTS)
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }
 
   for (var = 0; var < MoLNumEvolvedVariables; var++)
diff --git a/src/Operators.c b/src/Operators.c
new file mode 100644
index 0000000..e66312e
--- /dev/null
+++ b/src/Operators.c
@@ -0,0 +1,312 @@
+#include "Operators.h"
+#include <assert.h>
+#include <stddef.h>
+#include <stdlib.h>
+
+/* These are MoL's low-level operators. If they are overloaded as
+   aliased functions, these aliased functions are called; otherwise, a
+   default implementation is used. */
+
+/* The aliased functions should never be called from MoL's time
+   integrators, because they may not exist. Instead, the operators
+   defined here (with the MoL_ prefix) should be used. */
+
+static
+void
+error_no_storage(int const var, int const tl)
+{
+  char *const fullname = CCTK_FullName(var);
+  CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+             "Variable %s, timelevel %d has no storage", fullname, tl);
+  free(fullname);
+  return;
+}
+
+
+
+/* Some common special cases to improve performance */
+static
+void
+op_real_set_0(CCTK_REAL *restrict const varptr,
+              ptrdiff_t const npoints)
+{
+#pragma omp parallel for
+  for (ptrdiff_t i=0; i<npoints; ++i) {
+    varptr[i] = 0.0;
+  }
+}
+
+static
+void
+op_real_set_1(CCTK_REAL *restrict const varptr,
+              CCTK_REAL const *restrict const srcptr0,
+              CCTK_REAL const fact0,
+              ptrdiff_t const npoints)
+{
+#pragma omp parallel for
+  for (ptrdiff_t i=0; i<npoints; ++i) {
+    varptr[i] = fact0 * srcptr0[i];
+  }
+}
+
+static
+void
+op_real_set_2(CCTK_REAL *restrict const varptr,
+              CCTK_REAL const *restrict const srcptr0,
+              CCTK_REAL const fact0,
+              CCTK_REAL const *restrict const srcptr1,
+              CCTK_REAL const fact1,
+              ptrdiff_t const npoints)
+{
+#pragma omp parallel for
+  for (ptrdiff_t i=0; i<npoints; ++i) {
+    varptr[i] = fact0 * srcptr0[i] + fact1 * srcptr1[i];
+  }
+}
+
+static
+void
+op_real_set_3(CCTK_REAL *restrict const varptr,
+              CCTK_REAL const *restrict const srcptr0,
+              CCTK_REAL const fact0,
+              CCTK_REAL const *restrict const srcptr1,
+              CCTK_REAL const fact1,
+              CCTK_REAL const *restrict const srcptr2,
+              CCTK_REAL const fact2,
+              ptrdiff_t const npoints)
+{
+#pragma omp parallel for
+  for (ptrdiff_t i=0; i<npoints; ++i) {
+    varptr[i] = fact0 * srcptr0[i] + fact1 * srcptr1[i] + fact2 * srcptr2[i];
+  }
+}
+
+static
+void
+op_real_update_0(CCTK_REAL *restrict const varptr,
+                 CCTK_REAL const scale,
+                 ptrdiff_t const npoints)
+{
+#pragma omp parallel for
+  for (ptrdiff_t i=0; i<npoints; ++i) {
+    varptr[i] = scale * varptr[i];
+  }
+}
+
+static
+void
+op_real_update_1(CCTK_REAL *restrict const varptr,
+                 CCTK_REAL const scale,
+                 CCTK_REAL const *restrict const srcptr0,
+                 CCTK_REAL const fact0,
+                 ptrdiff_t const npoints)
+{
+#pragma omp parallel for
+  for (ptrdiff_t i=0; i<npoints; ++i) {
+    varptr[i] = scale * varptr[i] + fact0 * srcptr0[i];
+  }
+}
+
+static
+void
+op_real_update_2(CCTK_REAL *restrict const varptr,
+                 CCTK_REAL const scale,
+                 CCTK_REAL const *restrict const srcptr0,
+                 CCTK_REAL const fact0,
+                 CCTK_REAL const *restrict const srcptr1,
+                 CCTK_REAL const fact1,
+                 ptrdiff_t const npoints)
+{
+#pragma omp parallel for
+  for (ptrdiff_t i=0; i<npoints; ++i) {
+    varptr[i] = scale * varptr[i] + fact0 * srcptr0[i] + fact1 * srcptr1[i];
+  }
+}
+
+static
+void
+op_real_update_3(CCTK_REAL *restrict const varptr,
+                 CCTK_REAL const scale,
+                 CCTK_REAL const *restrict const srcptr0,
+                 CCTK_REAL const fact0,
+                 CCTK_REAL const *restrict const srcptr1,
+                 CCTK_REAL const fact1,
+                 CCTK_REAL const *restrict const srcptr2,
+                 CCTK_REAL const fact2,
+                 ptrdiff_t const npoints)
+{
+#pragma omp parallel for
+  for (ptrdiff_t i=0; i<npoints; ++i) {
+    varptr[i] =
+      scale * varptr[i] +
+      fact0 * srcptr0[i] + fact1 * srcptr1[i] + fact2 * srcptr2[i];
+  }
+}
+
+
+
+CCTK_INT
+MoL_LinearCombination(cGH const *const cctkGH,
+                      CCTK_INT   const var,
+                      CCTK_REAL  const scale,
+                      CCTK_INT   const srcs[],
+                      CCTK_INT   const tls[],
+                      CCTK_REAL  const facts[],
+                      CCTK_INT   const nsrcs)
+{
+  // Forward call to aliased function, if it is defined
+  static int is_aliased = -1;
+  if (is_aliased < 0) {
+    is_aliased = CCTK_IsFunctionAliased("LinearCombination");
+  }
+  if (is_aliased) {
+    return LinearCombination(cctkGH, var, scale, srcs, tls, facts, nsrcs);
+  }
+  
+  // Determine grid variable size
+  int const dim = CCTK_GroupDimFromVarI(var);
+  int ash[dim];
+  int const ierr = CCTK_GroupashVI(cctkGH, dim, ash, var);
+  assert(!ierr);
+  // TODO: check that all src variables have the same ash
+  ptrdiff_t npoints = 1;
+  for (int d=0; d<dim; ++d) {
+    npoints *= ash[d];
+  }
+  
+  switch (CCTK_VarTypeI(var)) {
+    
+  case CCTK_VARIABLE_REAL: {
+    // Obtain pointer to variable data
+    // TODO: check that all variable types are CCTK_REAL
+    CCTK_REAL *restrict const varptr = CCTK_VarDataPtrI(cctkGH, 0, var);
+    if (!varptr) error_no_storage(var, 0);
+    CCTK_REAL const *restrict srcptrs[nsrcs];
+    for (int n=0; n<nsrcs; ++n) {
+      srcptrs[n] = CCTK_VarDataPtrI(cctkGH, tls[n], srcs[n]);
+      if (!srcptrs[n]) error_no_storage(srcs[n], tls[n]);
+    }
+    
+    if (scale == 0.0) {
+      // Set (overwrite) target variable
+      
+      // Introduce special cases for some common cases to improve
+      // performance
+      switch (nsrcs) {
+      case 0:
+        op_real_set_0(varptr, npoints);
+        break;
+      case 1:
+        op_real_set_1(varptr, srcptrs[0], facts[0], npoints);
+        break;
+      case 2:
+        op_real_set_2(varptr,
+                      srcptrs[0], facts[0], srcptrs[1], facts[1], npoints);
+        break;
+      case 3:
+        op_real_set_3(varptr,
+                      srcptrs[0], facts[0], srcptrs[1], facts[1],
+                      srcptrs[2], facts[2], npoints);
+        break;
+      default:
+        // Loop over all grid points
+#pragma omp parallel for
+        for (ptrdiff_t i=0; i<npoints; ++i) {
+          CCTK_REAL tmp = 0.0;
+          for (int n=0; n<nsrcs; ++n) {
+            tmp += facts[n] * srcptrs[n][i];
+          }
+          varptr[i] = tmp;
+        }
+        break;
+      }
+      
+    } else {
+      // Update (add to) target variable
+      
+      // Introduce special cases for some common cases to improve
+      // performance
+      switch (nsrcs) {
+      case 0:
+        op_real_update_0(varptr, scale, npoints);
+        break;
+      case 1:
+        op_real_update_1(varptr, scale, srcptrs[0], facts[0], npoints);
+        break;
+      case 2:
+        op_real_update_2(varptr, scale,
+                         srcptrs[0], facts[0], srcptrs[1], facts[1], npoints);
+        break;
+      case 3:
+        op_real_update_3(varptr, scale,
+                         srcptrs[0], facts[0], srcptrs[1], facts[1],
+                         srcptrs[2], facts[2], npoints);
+        break;
+      default:
+        // Loop over all grid points
+#pragma omp parallel for
+        for (ptrdiff_t i=0; i<npoints; ++i) {
+          CCTK_REAL tmp = scale * varptr[i];
+          for (int n=0; n<nsrcs; ++n) {
+            tmp += facts[n] * srcptrs[n][i];
+          }
+          varptr[i] = tmp;
+        }
+        break;
+      }
+      
+    }
+    break;
+  }
+    
+  case CCTK_VARIABLE_COMPLEX: {
+    // Obtain pointer to variable data
+    // TODO: check that all variable types are CCTK_COMPLEX
+    CCTK_COMPLEX *restrict const varptr = CCTK_VarDataPtrI(cctkGH, 0, var);
+    if (!varptr) error_no_storage(var, 0);
+    CCTK_COMPLEX const *restrict srcptrs[nsrcs];
+    for (int n=0; n<nsrcs; ++n) {
+      srcptrs[n] = CCTK_VarDataPtrI(cctkGH, tls[n], srcs[n]);
+      if (!srcptrs[n]) error_no_storage(srcs[n], tls[n]);
+    }
+    
+    if (scale == 0.0) {
+      // Set (overwrite) target variable
+      // Loop over all grid points
+#pragma omp parallel for
+      for (ptrdiff_t i=0; i<npoints; ++i) {
+        CCTK_COMPLEX tmp = CCTK_Cmplx(0.0, 0.0);
+        for (int n=0; n<nsrcs; ++n) {
+          tmp = CCTK_CmplxAdd(tmp, CCTK_CmplxMul(CCTK_Cmplx(facts[n], 0.0),
+                                                 srcptrs[n][i]));
+        }
+        varptr[i] = tmp;
+      }
+    } else {
+      // Update (add to) target variable
+#pragma omp parallel for
+      for (ptrdiff_t i=0; i<npoints; ++i) {
+        CCTK_COMPLEX tmp = CCTK_CmplxMul(CCTK_Cmplx(scale, 0.0), varptr[i]);
+        for (int n=0; n<nsrcs; ++n) {
+          tmp = CCTK_CmplxAdd(tmp, CCTK_CmplxMul(CCTK_Cmplx(facts[n], 0.0),
+                                                 srcptrs[n][i]));
+        }
+        varptr[i] = tmp;
+      }
+    }
+    break;
+  }
+    
+  default:
+    // Other types (e.g. CCTK_REAL4) could be supported as well
+    CCTK_WARN(CCTK_WARN_ABORT, "Unsupported variable type");
+  }
+  
+  if (CCTK_IsFunctionAliased("Accelerator_NotifyDataModified")) {
+    CCTK_INT const tl = 0;
+    Accelerator_NotifyDataModified(cctkGH, &var, &tl, 1, 0);
+  }
+
+  // Done
+  return 0;
+}
diff --git a/src/Operators.h b/src/Operators.h
new file mode 100644
index 0000000..72618e4
--- /dev/null
+++ b/src/Operators.h
@@ -0,0 +1,15 @@
+#ifndef OPERATORS_H
+#define OPERATORS_H
+
+#include <cctk.h>
+
+CCTK_INT
+MoL_LinearCombination(cGH const *const cctkGH,
+                      CCTK_INT   const var,
+                      CCTK_REAL  const scale,
+                      CCTK_INT   const srcs[],
+                      CCTK_INT   const tls[],
+                      CCTK_REAL  const facts[],
+                      CCTK_INT   const nsrcs);
+
+#endif  // #ifndef OPERATORS_H
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index b80aa1b..96494bf 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -116,6 +116,12 @@ void MoL_ParamCheck(CCTK_ARGUMENTS)
     }
   }
 
+  if ( (CCTK_Equals(ODE_Method, "Euler"))&&(!(MoL_Intermediate_Steps == 1)) )
+  {
+    CCTK_PARAMWARN("When using the Euler evolver the "
+                   "number of intermediate steps must be 1");
+  }
+
   if ( (CCTK_Equals(ODE_Method, "RK2"))&&(!(MoL_Intermediate_Steps == 2)) )
   {
     CCTK_PARAMWARN("When using the efficient RK2 evolver the "
@@ -160,6 +166,12 @@ void MoL_ParamCheck(CCTK_ARGUMENTS)
                    "and the number of scratch levels at least 13.");
   }
   
+  if ( (CCTK_Equals(ODE_Method, "AB"))&&(!(MoL_Intermediate_Steps == 1)) )
+  {
+    CCTK_PARAMWARN("When using the Adams-Bashforth evolver the "
+                   "number of intermediate steps must be 1");
+  }
+
   if ( (CCTK_Equals(ODE_Method, "RK2-MR-2:1"))&&( !((MoL_Intermediate_Steps == 5) || (MoL_Num_Scratch_Levels > 4))) )
   {
     CCTK_PARAMWARN("When using the multirate 2-1 RK2 evolver the "
diff --git a/src/RK2-MR-2_1.c b/src/RK2-MR-2_1.c
index 8eaabd2..9179cdd 100644
--- a/src/RK2-MR-2_1.c
+++ b/src/RK2-MR-2_1.c
@@ -99,7 +99,7 @@ CCTK_WARN(0, "not implemented");
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-     totalsize *= cctk_lsh[arraydim];
+     totalsize *= cctk_ash[arraydim];
   }
 
   if (scratchspace_firstindex == -99)
diff --git a/src/RK2.c b/src/RK2.c
index 58dd03f..ea4f12c 100644
--- a/src/RK2.c
+++ b/src/RK2.c
@@ -16,6 +16,7 @@
 #include "cctk_Parameters.h"
 
 #include "ExternalVariables.h"
+#include "Operators.h"
 
 static const char *rcsid = "$Header$";
 
@@ -68,15 +69,10 @@ void MoL_RK2Add(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
 
-  cGroupDynamicData arraydata;
-  CCTK_INT groupindex, ierr;
-  CCTK_INT arraytotalsize, arraydim;
+  CCTK_INT arraydim;
   
-  CCTK_INT index, var;
+  CCTK_INT var;
   CCTK_INT totalsize;
-  CCTK_REAL const * restrict OldVar;
-  CCTK_REAL       * restrict UpdateVar;
-  CCTK_REAL const * restrict RHSVar;
   
   /* FIXME */
 
@@ -103,7 +99,7 @@ void MoL_RK2Add(CCTK_ARGUMENTS)
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }
 
   switch (*MoL_Intermediate_Step)
@@ -113,45 +109,24 @@ void MoL_RK2Add(CCTK_ARGUMENTS)
       {
         for (var = 0; var < MoLNumEvolvedVariables; var++)
         {
-          UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                   EvolvedVariableIndex[var]);
-          RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                RHSVariableIndex[var]);
-          
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateVar[index] += CCTK_DELTA_TIME * RHSVar[index];
-          }
+          int const nsrcs = 1;
+          CCTK_INT const srcs[] = {RHSVariableIndex[var]};
+          CCTK_INT const tls[] = {0};
+          CCTK_REAL const facts[] = {CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedVariableIndex[var], 1.0,
+                                srcs, tls, facts, nsrcs);
         }
 
         for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
         {
-          UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                   EvolvedArrayVariableIndex[var]);
-          RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                RHSArrayVariableIndex[var]);
-          groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
-          ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
-                                       &arraydata);
-          if (ierr)
-          {
-            CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, 
-                       "The driver does not return group information "
-                       "for group '%s'.", 
-                       CCTK_GroupName(groupindex));
-          }
-          arraytotalsize = 1;
-          for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
-          {
-            arraytotalsize *= arraydata.lsh[arraydim];
-          }
-
-#pragma omp parallel for
-          for (index = 0; index < arraytotalsize; index++)
-          {
-            UpdateVar[index] += CCTK_DELTA_TIME * RHSVar[index];
-          }
+          int const nsrcs = 1;
+          CCTK_INT const srcs[] = {RHSArrayVariableIndex[var]};
+          CCTK_INT const tls[] = {0};
+          CCTK_REAL const facts[] = {CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedArrayVariableIndex[var], 1.0,
+                                srcs, tls, facts, nsrcs);
         }
 
   /* FIXME */
@@ -160,18 +135,13 @@ void MoL_RK2Add(CCTK_ARGUMENTS)
 
         for (var = 0; var < MoLNumEvolvedComplexVariables; var++)
         {
-          UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                             EvolvedComplexVariableIndex[var]);
-          RHSComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                          RHSComplexVariableIndex[var]);
-          
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateComplexVar[index] = CCTK_CmplxAdd(UpdateComplexVar[index],
-                                                    CCTK_CmplxMul(Complex_Delta_Time,
-                                                                  RHSComplexVar[index]));
-          }
+          int const nsrcs = 1;
+          CCTK_INT const srcs[] = {RHSComplexVariableIndex[var]};
+          CCTK_INT const tls[] = {0};
+          CCTK_REAL const facts[] = {CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedComplexVariableIndex[var], 1.0,
+                                srcs, tls, facts, nsrcs);
         }
 
 #endif
@@ -182,52 +152,26 @@ void MoL_RK2Add(CCTK_ARGUMENTS)
       {
         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, 
-                                                RHSVariableIndex[var]);
-          
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateVar[index] = 0.5 * (OldVar[index] + UpdateVar[index]) + 
-                                      CCTK_DELTA_TIME * RHSVar[index];
-          }
+          int const nsrcs = 2;
+          CCTK_INT const srcs[] =
+            {EvolvedVariableIndex[var], RHSVariableIndex[var]};
+          CCTK_INT const tls[] = {1, 0};
+          CCTK_REAL const facts[] = {0.5, CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedVariableIndex[var], 0.5,
+                                srcs, tls, facts, nsrcs);
         }
 
         for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
         {
-          OldVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 1,
-                                                   EvolvedArrayVariableIndex[var]);
-          UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                   EvolvedArrayVariableIndex[var]);
-          RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                RHSArrayVariableIndex[var]);
-
-          groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
-          ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
-                                       &arraydata);
-          if (ierr)
-          {
-            CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, 
-                       "The driver does not return group information "
-                       "for group '%s'.", 
-                       CCTK_GroupName(groupindex));
-          }
-          arraytotalsize = 1;
-          for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
-          {
-            arraytotalsize *= arraydata.lsh[arraydim];
-          }
-          
-#pragma omp parallel for
-          for (index = 0; index < arraytotalsize; index++)
-          {
-            UpdateVar[index] = 0.5 * (OldVar[index] + UpdateVar[index]) + 
-                                      CCTK_DELTA_TIME * RHSVar[index];
-          }
+          int const nsrcs = 2;
+          CCTK_INT const srcs[] =
+            {EvolvedArrayVariableIndex[var], RHSArrayVariableIndex[var]};
+          CCTK_INT const tls[] = {1, 0};
+          CCTK_REAL const facts[] = {0.5, CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedArrayVariableIndex[var], 0.5,
+                                srcs, tls, facts, nsrcs);
         }
 
   /* FIXME */
@@ -236,23 +180,14 @@ void MoL_RK2Add(CCTK_ARGUMENTS)
 
         for (var = 0; var < MoLNumEvolvedComplexVariables; var++)
         {
-          OldComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 1,
-                                                       EvolvedComplexVariableIndex[var]);
-          UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                             EvolvedComplexVariableIndex[var]);
-          RHSComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                          RHSComplexVariableIndex[var]);
-          
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateComplexVar[index] =
-              CCTK_CmplxAdd(CCTK_CmplxMul(Complex_Half,
-                                          (CCTK_CmplxAdd(OldComplexVar[index],
-                                                         UpdateComplexVar[index]))), 
-                            CCTK_CmplxMul(Complex_Delta_Time, 
-                                          RHSComplexVar[index])); 
-          }
+          int const nsrcs = 2;
+          CCTK_INT const srcs[] =
+            {EvolvedComplexVariableIndex[var], RHSComplexVariableIndex[var]};
+          CCTK_INT const tls[] = {1, 0};
+          CCTK_REAL const facts[] = {0.5, CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedComplexVariableIndex[var], 0.5,
+                                srcs, tls, facts, nsrcs);
         }
 
 #endif
diff --git a/src/RK3.c b/src/RK3.c
index f3b5e9c..38a55ab 100644
--- a/src/RK3.c
+++ b/src/RK3.c
@@ -15,6 +15,7 @@
 #include "cctk_Parameters.h"
 
 #include "ExternalVariables.h"
+#include "Operators.h"
 
 static const char *rcsid = "$Header$";
 
@@ -67,15 +68,10 @@ void MoL_RK3Add(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
 
-  cGroupDynamicData arraydata;
-  CCTK_INT groupindex, ierr;
-  CCTK_INT arraytotalsize, arraydim;
+  CCTK_INT arraydim;
   
-  CCTK_INT index, var;
+  CCTK_INT var;
   CCTK_INT totalsize;
-  CCTK_REAL const * restrict OldVar;
-  CCTK_REAL       * restrict UpdateVar;
-  CCTK_REAL const * restrict RHSVar;
   
   /* FIXME */
 
@@ -106,7 +102,7 @@ void MoL_RK3Add(CCTK_ARGUMENTS)
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }
 
   switch (*MoL_Intermediate_Step)
@@ -116,45 +112,24 @@ void MoL_RK3Add(CCTK_ARGUMENTS)
       {
         for (var = 0; var < MoLNumEvolvedVariables; var++)
         {
-          UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                   EvolvedVariableIndex[var]);
-          RHSVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                RHSVariableIndex[var]);
-          
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateVar[index] += CCTK_DELTA_TIME * RHSVar[index];
-          }
+          int const nsrcs = 1;
+          CCTK_INT const srcs[] = {RHSVariableIndex[var]};
+          CCTK_INT const tls[] = {0};
+          CCTK_REAL const facts[] = {CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedVariableIndex[var], 1.0,
+                                srcs, tls, facts, nsrcs);
         }
 
         for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
         {
-          UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                   EvolvedArrayVariableIndex[var]);
-          RHSVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                RHSArrayVariableIndex[var]);
-          groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
-          ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
-                                       &arraydata);
-          if (ierr)
-          {
-            CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, 
-                       "The driver does not return group information "
-                       "for group '%s'.", 
-                       CCTK_GroupName(groupindex));
-          }
-          arraytotalsize = 1;
-          for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
-          {
-            arraytotalsize *= arraydata.lsh[arraydim];
-          }
-
-#pragma omp parallel for
-          for (index = 0; index < arraytotalsize; index++)
-          {
-            UpdateVar[index] += CCTK_DELTA_TIME * RHSVar[index];
-          }
+          int const nsrcs = 1;
+          CCTK_INT const srcs[] = {RHSArrayVariableIndex[var]};
+          CCTK_INT const tls[] = {0};
+          CCTK_REAL const facts[] = {CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedArrayVariableIndex[var], 1.0,
+                                srcs, tls, facts, nsrcs);
         }
 
   /* FIXME */
@@ -163,18 +138,13 @@ void MoL_RK3Add(CCTK_ARGUMENTS)
 
         for (var = 0; var < MoLNumEvolvedComplexVariables; var++)
         {
-          UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                             EvolvedComplexVariableIndex[var]);
-          RHSComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                          RHSComplexVariableIndex[var]);
-          
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateComplexVar[index] = CCTK_CmplxAdd(UpdateComplexVar[index],
-                                                    CCTK_CmplxMul(Complex_Delta_Time,
-                                                                  RHSComplexVar[index]));
-          }
+          int const nsrcs = 1;
+          CCTK_INT const srcs[] = {RHSComplexVariableIndex[var]};
+          CCTK_INT const tls[] = {0};
+          CCTK_REAL const facts[] = {CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedComplexVariableIndex[var], 1.0,
+                                srcs, tls, facts, nsrcs);
         }
 
 #endif
@@ -186,52 +156,26 @@ void MoL_RK3Add(CCTK_ARGUMENTS)
       {
         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, 
-                                                RHSVariableIndex[var]);
-          
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateVar[index] = 0.25 * (3*OldVar[index] +
-                                       UpdateVar[index]) + 
-              CCTK_DELTA_TIME * RHSVar[index]; 
-          }
+          int const nsrcs = 2;
+          CCTK_INT const srcs[] =
+            {EvolvedVariableIndex[var], RHSVariableIndex[var]};
+          CCTK_INT const tls[] = {1, 0};
+          CCTK_REAL const facts[] = {0.75, CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedVariableIndex[var], 0.25,
+                                srcs, tls, facts, nsrcs);
         }
 
         for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
         {
-          OldVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 1,
-                                                EvolvedArrayVariableIndex[var]);
-          UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                   EvolvedArrayVariableIndex[var]);
-          RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                RHSArrayVariableIndex[var]);
-          groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
-          ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
-                                       &arraydata);
-          if (ierr)
-          {
-            CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, 
-                       "The driver does not return group information "
-                       "for group '%s'.", 
-                       CCTK_GroupName(groupindex));
-          }
-          arraytotalsize = 1;
-          for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
-          {
-            arraytotalsize *= arraydata.lsh[arraydim];
-          }
-
-#pragma omp parallel for
-          for (index = 0; index < arraytotalsize; index++)
-          {
-            UpdateVar[index] = 0.25*(3*OldVar[index] + UpdateVar[index]) +
-              CCTK_DELTA_TIME * RHSVar[index]; 
-          }
+          int const nsrcs = 2;
+          CCTK_INT const srcs[] =
+            {EvolvedArrayVariableIndex[var], RHSArrayVariableIndex[var]};
+          CCTK_INT const tls[] = {1, 0};
+          CCTK_REAL const facts[] = {0.75, CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedArrayVariableIndex[var], 0.25,
+                                srcs, tls, facts, nsrcs);
         }
 
   /* FIXME */
@@ -240,24 +184,14 @@ void MoL_RK3Add(CCTK_ARGUMENTS)
 
         for (var = 0; var < MoLNumEvolvedComplexVariables; var++)
         {
-          OldComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 1,
-                                                          EvolvedComplexVariableIndex[var]);
-          UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                             EvolvedComplexVariableIndex[var]);
-          RHSComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                          RHSComplexVariableIndex[var]);
-          
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateComplexVar[index] =
-              CCTK_CmplxAdd(CCTK_CmplxAdd(CCTK_CmplxMul(Complex_ThreeQuarter,
-                                                        OldComplexVar[index]),
-                                          CCTK_CmplxMul(ComplexQuarter,
-                                                        UpdateComplexVar[index])), 
-                            CCTK_CmplxMul(CCTK_CmplxMul(Complex_Delta_Time, 
-                                                        RHSComplexVar[index]))); 
-          }
+          int const nsrcs = 2;
+          CCTK_INT const srcs[] =
+            {EvolvedComplexVariableIndex[var], RHSComplexVariableIndex[var]};
+          CCTK_INT const tls[] = {1, 0};
+          CCTK_REAL const facts[] = {0.75, CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedComplexVariableIndex[var], 0.25,
+                                srcs, tls, facts, nsrcs);
         }
 
 #endif
@@ -266,56 +200,28 @@ void MoL_RK3Add(CCTK_ARGUMENTS)
       }
     case 1:
       {
-        const CCTK_REAL one_third = 1.0 / 3.0;
-
         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, 
-                                                RHSVariableIndex[var]);
-          
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateVar[index] = (OldVar[index] + 2*UpdateVar[index]) * one_third
-              + CCTK_DELTA_TIME * RHSVar[index];
-          }
+          int const nsrcs = 2;
+          CCTK_INT const srcs[] =
+            {EvolvedVariableIndex[var], RHSVariableIndex[var]};
+          CCTK_INT const tls[] = {1, 0};
+          CCTK_REAL const facts[] = {1.0/3.0, CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedVariableIndex[var], 2.0/3.0,
+                                srcs, tls, facts, nsrcs);
         }
 
         for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
         {
-          OldVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 1,
-                                                   EvolvedArrayVariableIndex[var]);
-          UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                   EvolvedArrayVariableIndex[var]);
-          RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                RHSArrayVariableIndex[var]);
-
-          groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
-          ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
-                                       &arraydata);
-          if (ierr)
-          {
-            CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, 
-                       "The driver does not return group information "
-                       "for group '%s'.", 
-                       CCTK_GroupName(groupindex));
-          }
-          arraytotalsize = 1;
-          for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
-          {
-            arraytotalsize *= arraydata.lsh[arraydim];
-          }
-          
-#pragma omp parallel for
-          for (index = 0; index < arraytotalsize; index++)
-          {
-            UpdateVar[index] = (OldVar[index] + 2*UpdateVar[index]) * one_third
-              + CCTK_DELTA_TIME * RHSVar[index];
-          }
+          int const nsrcs = 2;
+          CCTK_INT const srcs[] =
+            {EvolvedArrayVariableIndex[var], RHSArrayVariableIndex[var]};
+          CCTK_INT const tls[] = {1, 0};
+          CCTK_REAL const facts[] = {1.0/3.0, CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedArrayVariableIndex[var], 2.0/3.0,
+                                srcs, tls, facts, nsrcs);
         }
 
   /* FIXME */
@@ -324,24 +230,14 @@ void MoL_RK3Add(CCTK_ARGUMENTS)
 
         for (var = 0; var < MoLNumEvolvedComplexVariables; var++)
         {
-          OldComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 1,
-                                                          EvolvedComplexVariableIndex[var]);
-          UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
-                                                             EvolvedComplexVariableIndex[var]);
-          RHSComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 0, 
-                                                          RHSComplexVariableIndex[var]);
-          
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateComplexVar[index] =
-              CCTK_CmplxAdd(CCTK_CmplxAdd(CCTK_CmplxMul(Complex_Third,
-                                                        OldComplexVar[index]),
-                                          CCTK_CmplxMul(ComplexTwoThird,
-                                                        UpdateComplexVar[index])), 
-                            CCTK_CmplxMul(Complex_Delta_Time, 
-                                          RHSComplexVar[index])); 
-          }
+          int const nsrcs = 2;
+          CCTK_INT const srcs[] =
+            {EvolvedComplexVariableIndex[var], RHSComplexVariableIndex[var]};
+          CCTK_INT const tls[] = {1, 0};
+          CCTK_REAL const facts[] = {1.0/3.0, CCTK_DELTA_TIME};
+          MoL_LinearCombination(cctkGH,
+                                EvolvedComplexVariableIndex[var], 2.0/3.0,
+                                srcs, tls, facts, nsrcs);
         }
 
 #endif
diff --git a/src/RK4-MR-2_1.c b/src/RK4-MR-2_1.c
index 815a41e..0a48cb0 100644
--- a/src/RK4-MR-2_1.c
+++ b/src/RK4-MR-2_1.c
@@ -99,7 +99,7 @@ CCTK_WARN(0, "not implemented");
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-     totalsize *= cctk_lsh[arraydim];
+     totalsize *= cctk_ash[arraydim];
   }
 
   if (scratchspace_firstindex == -99)
diff --git a/src/RK4-RK2.c b/src/RK4-RK2.c
index 6e5b85e..472418e 100644
--- a/src/RK4-RK2.c
+++ b/src/RK4-RK2.c
@@ -54,7 +54,7 @@ void MoL_RK4_RK2_Add(CCTK_ARGUMENTS)
   int const step = MoL_Intermediate_Steps - *MoL_Intermediate_Step;
   
   int totalsize = 1;
-  for (int d=0; d<cctk_dim; ++d) totalsize *= cctk_lsh[d];
+  for (int d=0; d<cctk_dim; ++d) totalsize *= cctk_ash[d];
   
   CCTK_REAL const dt = *Original_Delta_Time / cctkGH->cctk_timefac;
   
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;
diff --git a/src/RK45.c b/src/RK45.c
index ea547fd..2e520da 100644
--- a/src/RK45.c
+++ b/src/RK45.c
@@ -159,7 +159,7 @@ void MoL_RK45Add(CCTK_ARGUMENTS)
   int totalsize = 1;
   for (int arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }  
 
   /* Real GFs */
diff --git a/src/RK65.c b/src/RK65.c
index f333741..e5dc629 100644
--- a/src/RK65.c
+++ b/src/RK65.c
@@ -115,7 +115,7 @@ void MoL_RK65Add(CCTK_ARGUMENTS)
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }  
 
   /* Real GFs */
diff --git a/src/RK87.c b/src/RK87.c
index dee29da..68c7ff2 100644
--- a/src/RK87.c
+++ b/src/RK87.c
@@ -14,6 +14,7 @@
 #include "cctk_Parameters.h"
 
 #include "ExternalVariables.h"
+#include "Operators.h"
 
 static const char *rcsid = "$Header$";
 
@@ -70,17 +71,9 @@ void MoL_RK87Add(CCTK_ARGUMENTS)
   CCTK_INT arraydim;
 
   static CCTK_INT scratchspace_firstindex = -99;
-  CCTK_INT index, var, scratchstep;
+  CCTK_INT var, scratchstep;
   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_REAL beta, gamma, gamma_error;
-
   static const CCTK_REAL beta_array[12][12] = {
     { 1.0/18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
     { 1.0/48.0, 1.0/16.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
@@ -131,7 +124,7 @@ void MoL_RK87Add(CCTK_ARGUMENTS)
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }  
 
   if (scratchspace_firstindex == -99)
@@ -145,122 +138,101 @@ void MoL_RK87Add(CCTK_ARGUMENTS)
 
   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)));
-#pragma omp parallel for
-    for (index = 0; index < totalsize; index++)
-    {
-      ScratchVar[index] = tmp * RHSVar[index];
-    }
-
+    int const step = MoL_Intermediate_Steps - (*MoL_Intermediate_Step);
+    int const scratchvarindex =
+      scratchspace_firstindex + MoL_Num_Evolved_Vars * step + var;
+
+    int const nsrcs = 1;
+    CCTK_INT const srcs[] = {RHSVariableIndex[var]};
+    CCTK_INT const tls[] = {0};
+    CCTK_REAL const facts[] = {(*Original_Delta_Time) / cctkGH->cctk_timefac};
+    MoL_LinearCombination(cctkGH,
+                          scratchvarindex, 1.0,
+                          srcs, tls, facts, nsrcs);
   }
 
-
   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);
+    int const error_firstindex = CCTK_FirstVarIndex("MOL::ERRORESTIMATE");
+    int const step = MoL_Intermediate_Steps - (*MoL_Intermediate_Step);
 
     if (*MoL_Intermediate_Step - 1)
-    { 
-
-#pragma omp parallel for
-      for (index = 0; index < totalsize; index++)
+    {
+      int const nsrcs = step + 2;
+      CCTK_INT srcs[nsrcs];
+      CCTK_INT tls[nsrcs];
+      CCTK_REAL facts[nsrcs];
+      srcs[0]  = EvolvedVariableIndex[var];
+      tls[0]   = 1;
+      facts[0] = 1.0;
+
+      for (scratchstep = 0; scratchstep < step + 1; scratchstep++)
       {
-        UpdateVar[index] = OldVar[index];
+        int const scratchvarindex =
+          scratchspace_firstindex + MoL_Num_Evolved_Vars * scratchstep + var;
+        srcs [scratchstep+1] = scratchvarindex;
+        tls  [scratchstep+1] = 0;
+        facts[scratchstep+1] = beta_array[step][scratchstep];
       }
 
-      for (scratchstep = 0; 
-           scratchstep < MoL_Intermediate_Steps - (*MoL_Intermediate_Step) + 1;
-           scratchstep++)
+      MoL_LinearCombination(cctkGH,
+                            EvolvedVariableIndex[var], 0.0,
+                            srcs, tls, facts, nsrcs);
+    }
+    else
+    {
       {
+        int const nsrcs = 14;
+        CCTK_INT srcs[nsrcs];
+        CCTK_INT tls[nsrcs];
+        CCTK_REAL facts[nsrcs];
+        srcs[0]  = EvolvedVariableIndex[var];
+        tls[0]   = 1;
+        facts[0] = 1.0;
         
-        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 (scratchstep = 0; scratchstep < 13; scratchstep++)
         {
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateVar[index] += beta * ScratchVar[index];
-          }
+          int const scratchvarindex =
+            scratchspace_firstindex + MoL_Num_Evolved_Vars * scratchstep + var;
+          srcs [scratchstep+1] = scratchvarindex;
+          tls  [scratchstep+1] = 0;
+          facts[scratchstep+1] = gamma_array[scratchstep];
         }
-        
-      }
-      
-    }
-    else
-    {
 
-#pragma omp parallel for
-      for (index = 0; index < totalsize; index++)
-      {
-        UpdateVar[index] = OldVar[index];
-        ErrorVar[index] = 0;
+        MoL_LinearCombination(cctkGH,
+                              EvolvedVariableIndex[var], 0.0,
+                              srcs, tls, facts, nsrcs);
       }
-    
-      for (scratchstep = 0; scratchstep < 13; scratchstep++)
+
       {
+        int const nsrcs = 13;
+        CCTK_INT srcs[nsrcs];
+        CCTK_INT tls[nsrcs];
+        CCTK_REAL facts[nsrcs];
         
-        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 (scratchstep = 0; scratchstep < 13; scratchstep++)
         {
-#pragma omp parallel for
-          for (index = 0; index < totalsize; index++)
-          {
-            UpdateVar[index] += gamma * ScratchVar[index];
-            ErrorVar[index] += gamma_error * ScratchVar[index];
-          }
+          int const scratchvarindex =
+            scratchspace_firstindex + MoL_Num_Evolved_Vars * scratchstep + var;
+          srcs [scratchstep] = scratchvarindex;
+          tls  [scratchstep] = 0;
+          facts[scratchstep] =
+            gamma_array[scratchstep] - gammastar_array[scratchstep];
         }
 
+        MoL_LinearCombination(cctkGH,
+                              error_firstindex + var, 0.0,
+                              srcs, tls, facts, nsrcs);
       }
-      
+
     }
-    
+
   }
 
   /* Real arrays */
 
-  for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
-  {
-
-    CCTK_WARN(0, "Peter has been too lazy to write the RK87 routine "
-              "out for array variables. Better send him an email...");
-
-  }
+  CCTK_WARN(0, "Peter has been too lazy to write the RK87 routine "
+            "out for array variables. Better send him an email...");
         
   return;
 }
diff --git a/src/SandR.c b/src/SandR.c
index e644c53..577a8fe 100644
--- a/src/SandR.c
+++ b/src/SandR.c
@@ -91,7 +91,7 @@ void MoL_RestoreSandR(CCTK_ARGUMENTS)
   totalsize = 1;
   for (arraydim = 0; arraydim < cctk_dim; arraydim++)
   {
-    totalsize *= cctk_lsh[arraydim];
+    totalsize *= cctk_ash[arraydim];
   }
 
   for (var = 0; var < MoLNumSandRVariables; var++) 
diff --git a/src/SetTime.c b/src/SetTime.c
index edef89f..b60f3bb 100644
--- a/src/SetTime.c
+++ b/src/SetTime.c
@@ -234,6 +234,13 @@ void MoL_ResetTime(CCTK_ARGUMENTS)
     cctkGH->cctk_time = previous_times[MoL_Intermediate_Steps - 
                                       *MoL_Intermediate_Step];
   }
+  else if (CCTK_EQUALS(ODE_Method,"Euler"))
+  {
+    if (*MoL_Intermediate_Step == 1)
+    {
+      cctkGH->cctk_time = (*Original_Time);
+    }
+  }
   else if (CCTK_EQUALS(ODE_Method,"RK2"))
   {
     if (*MoL_Intermediate_Step == 1)
@@ -433,6 +440,13 @@ void MoL_ResetDeltaTime(CCTK_ARGUMENTS)
                                                 (*MoL_Intermediate_Step)] *
       (*Original_Delta_Time);
   }
+  else if (CCTK_EQUALS(ODE_Method,"Euler"))
+  {
+    if (*MoL_Intermediate_Step == 1)
+    {
+      cctkGH->cctk_delta_time = (*Original_Delta_Time);
+    }
+  }
   else if (CCTK_EQUALS(ODE_Method,"RK2"))
   {
     if (*MoL_Intermediate_Step == 1)
diff --git a/src/make.code.defn b/src/make.code.defn
index 80e1b56..f11c36b 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,14 +2,20 @@
 # $Header$
 
 # Source files in this directory
-SRCS = 	ChangeType.c \
+SRCS = 	AB.c \
+	ChangeType.c \
 	Counter.c \
+	Euler.c \
 	GenericRK.c \
 	ICN.c \
 	IndexArrays.c \
 	InitialCopy.c \
+	Operators.c \
 	ParamCheck.c \
 	RK2.c \
+	RK2-MR-2_1.c \
+	RK4-MR-2_1.c \
+	RK4-RK2.c \
 	RK3.c \
 	RK4.c \
 	RK45.c \
@@ -21,11 +27,8 @@ SRCS = 	ChangeType.c \
 	SandR.c \
 	SetTime.c \
 	Startup.c \
-	StepSize.c \
-	RK2-MR-2_1.c \
-	RK4-MR-2_1.c \
-	RK4-RK2.c
-	
+	StepSize.c
+
 # Subdirectories containing source files
 SUBDIRS =