Merge branch opencl into master

22 files changed, 1715 insertions, 1181 deletions

diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/configuration.ccl b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/configuration.ccl
index 477b219..664b79a 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/configuration.ccl
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/configuration.ccl
@@ -7,3 +7,7 @@ PROVIDES GenericFD
   SCRIPT
   LANG
 }
+
+OPTIONAL Vectors
+{
+}
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
index 401b4e0..ace83a0 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
@@ -43,18 +43,30 @@ extern "C" {
 #include "MathematicaCompat.h"
 
 #ifdef KRANC_C
+
+  /* Grid function access */
+  /* var is a pointer to a grid point, i,j,k are offsets with respect
+     to that point.
+     For example: KRANC_GFINDEX3D_OFFSET(&u[ind3d],-1,-1,0) */
+#ifndef VECTORISE
+  /* standard, thorn Vectors is not used */
+  /* simple implementation */
+  /* #  define KRANC_GFOFFSET3D(var,i,j,k) ((var)[di*(i)+dj*(j)+dk*(k)]) */
+  /* more efficient implementation for some compilers */
+#  define KRANC_GFOFFSET3D(var,i,j,k)                                   \
+  (*(CCTK_REAL const*)&((char const*)(var))[cdi*(i)+cdj*(j)+cdk*(k)])
+#else
+  /* vectorised version */
+#  define KRANC_GFOFFSET3D(var,i,j,k)                                   \
+  vec_loadu_maybe3((i),(j),(k),                                         \
+                   *(CCTK_REAL const*)&                                 \
+                   ((char const*)(var))[cdi*(i)+cdj*(j)+cdk*(k)])
+#endif
+
 int sgn(CCTK_REAL x);
 
 int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH);
 
-#ifdef __cplusplus
-// Define the restrict qualifier
-#  ifdef CCTK_CXX_RESTRICT
-#    undef restrict
-#    define restrict CCTK_CXX_RESTRICT
-#  endif
-#endif
-
 void GenericFD_GetBoundaryInfo(cGH const * restrict cctkGH,
                                int const * restrict cctk_lsh,
                                int const * restrict cctk_lssh,
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h
index 1823cd8..25ec32c 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h
@@ -30,6 +30,8 @@
 #define Cosh(x)       (cosh(x))
 #define Tanh(x)       (tanh(x))
 
+#define Csch(x)       (1./sinh(x))
+
 #ifdef KRANC_C
 #  define Sign(x)     (copysign(1.0,(x)))
 #  define ToReal(x)   ((CCTK_REAL)(x))
@@ -38,38 +40,6 @@
 #  define ToReal(x)   (real((x),kind(khalf)))
 #endif
 
-#if 0
-
-/* TODO: use fma(x,y,z) to implement fmadd and friends?  Note that fma
-   may be unsupported, or may be slow.  */
-
-/* #define fmadd(x,y,z)  ((x)*(y)+(z)) */
-/* #define fmsub(x,y,z)  ((x)*(y)-(z)) */
-/* #define fnmadd(x,y,z) (-(z)-(x)*(y)) */
-/* #define fnmsub(x,y,z) (+(z)-(x)*(y)) */
-
-#define fpos(x)   (+(x))
-#define fneg(x)   (-(x))
-#define fmul(x,y) ((x)*(y))
-#define fdiv(x,y) ((x)/(y))
-#define fadd(x,y) ((x)+(y))
-#define fsub(x,y) ((x)-(y))
-
-#define fmadd(x,y,z)  (fadd(fmul(x,y),z))
-#define fmsub(x,y,z)  (fsub(fmul(x,y),z))
-#define fnmadd(x,y,z) (fsub(fneg(z),fmul(x,y)))
-#define fnmsub(x,y,z) (fsub(z,fmul(x,y)))
-
-#define kexp(x)    (exp(x))
-#define kfabs(x)   (fabs(x))
-#define kfmax(x,y) (fmax(x,y))
-#define kfmin(x,y) (fmin(x,y))
-#define klog(x)    (log(x))
-#define kpow(x,y)  (pow(x,y))
-#define ksqrt(x)   (sqrt(x))
-
-#endif
-
 #ifdef KRANC_C
 #  define E           M_E
 #  define Pi          M_PI
diff --git a/Bin/kranc b/Bin/kranc
index 9f16432..e82e75a 100755
--- a/Bin/kranc
+++ b/Bin/kranc
@@ -1,9 +1,22 @@
 #!/bin/bash
 
-# Assume that this script is called from the Kranc/Bin directory
+set -e
+
+# Assume that this script is called from the Kranc/Bin directory.
+# This will not work if someone creates a symlink to the kranc script
+# somewhere else
 export KRANCDIR=$(dirname $0)/..
+
+if [ ! -r "$KRANCDIR/Tools/CodeGen" ]; then
+    echo "Cannot find Kranc (the kranc script must be run directly from the Kranc/Bin directory - symbolic links are not currently allowed)"
+    exit 1
+fi
+
 export KRANCVERBOSE=no
 
+echo "Using Kranc installation at $KRANCDIR"
+# It would be good to find a portable way to canonicalise KRANCDIR.
+
 while getopts "v" flag
 do
   case $flag in
diff --git a/Examples/Laplace.m b/Examples/Laplace.m
new file mode 100644
index 0000000..ad8a28f
--- /dev/null
+++ b/Examples/Laplace.m
@@ -0,0 +1,68 @@
+
+groups = {{"phi_group", {phi}}};
+
+derivatives =
+{
+  PDstandard[i_] ->
+    StandardCenteredDifferenceOperator[1,fdOrder/2,i],
+  PDstandard[i_, i_] ->
+    StandardCenteredDifferenceOperator[2,fdOrder/2,i],
+  PDstandard[i_, j_] ->
+    StandardCenteredDifferenceOperator[1,fdOrder/2,i] StandardCenteredDifferenceOperator[1,fdOrder/2,j]
+};
+
+PD = PDstandard;
+
+initialCalc = 
+{
+  Name -> "Laplace_initial",
+  Schedule -> {"AT INITIAL"},
+  Where -> Interior,
+  Equations -> 
+  {
+    phi -> phi0 Sum[4/(Pi n) Sin[n Pi x/Lx] Sinh[n Pi y/Lx]/Sinh[n Pi Ly/Lx], {n, 1, 1, 2}]
+  }
+};
+
+initialBoundaryCalc = 
+{
+  Name -> "Laplace_initial_boundary",
+  Schedule -> {"AT INITIAL after Laplace_initial"},
+  Where -> Boundary,
+  Equations -> 
+  {
+    phi -> IfThen[Abs[y-Ly]<10^-10, phi0, 0]
+  }
+};
+
+evolveCalc = 
+{
+  Name -> "Laplace_relax",
+  Schedule -> {"in MoL_CalcRHS", "AT ANALYSIS"},
+  Equations ->
+  {
+    dot[phi] -> mu Euc[ui,uj] PD[phi,li,lj]
+  }
+};
+
+boundaryCalc = 
+{
+  Name -> "Laplace_boundary",
+  Schedule -> {"in MoL_RHSBoundaries", "AT ANALYSIS"},
+  Where -> Boundary,
+  Equations ->
+  {
+    dot[phi] -> 0
+  }
+};
+
+
+CreateKrancThornTT[
+  groups, ".", 
+  "Laplace", 
+  Calculations -> {initialCalc, initialBoundaryCalc, evolveCalc, boundaryCalc},
+  PartialDerivatives -> derivatives,
+  ZeroDimensions -> {3},
+  RealParameters -> {Lx,Ly,phi0,mu},
+  IntParameters -> {{Name -> fdOrder, Default -> 2, AllowedValues -> {2, 4}}},
+  DeclaredGroups -> {"phi_group"}];
diff --git a/Examples/laplace.par b/Examples/laplace.par
new file mode 100644
index 0000000..a63db80
--- /dev/null
+++ b/Examples/laplace.par
@@ -0,0 +1,91 @@
+
+Cactus::cctk_final_time = 1
+Cactus::terminate = "time"
+
+ActiveThorns = "IOUtil Carpet CarpetLib CarpetSlab CoordBase CoordBase SymBase CartGrid3D Slab CarpetIOBasic CarpetIOASCII CarpetIOScalar Laplace Time MoL Boundary GenericFD CarpetReduce LoopControl CarpetIOHDF5 HDF5"
+
+CoordBase::domainsize                   = minmax
+
+CoordBase::boundary_size_x_lower        = 1
+CoordBase::boundary_size_y_lower        = 1
+CoordBase::boundary_size_z_lower        = 0
+CoordBase::boundary_shiftout_x_lower    = 0
+CoordBase::boundary_shiftout_y_lower    = 0
+CoordBase::boundary_shiftout_z_lower    = 1
+
+CoordBase::boundary_size_x_upper        = 1
+CoordBase::boundary_size_y_upper        = 1
+CoordBase::boundary_size_z_upper        = 0
+CoordBase::boundary_shiftout_x_upper    = 0
+CoordBase::boundary_shiftout_y_upper    = 0
+CoordBase::boundary_shiftout_z_upper    = 1
+
+CartGrid3D::type                        = "coordbase"
+CartGrid3D::domain                      = "full"
+CartGrid3D::avoid_origin                = "no"
+
+CoordBase::xmin                         = 0
+CoordBase::ymin                         = 0
+CoordBase::zmin                         = 0
+CoordBase::xmax                         = 1
+CoordBase::ymax                         = 1
+CoordBase::zmax                         = 0
+CoordBase::dx                           = 0.02
+CoordBase::dy                           = 0.02
+CoordBase::dz                           = 0.02
+
+driver::ghost_size_x = 1
+driver::ghost_size_y = 1
+driver::ghost_size_z = 0
+
+Carpet::domain_from_coordbase           = "yes"
+Carpet::poison_new_timelevels           = "yes"
+Carpet::check_for_poison                = "no"
+
+Carpet::max_refinement_levels = 1
+Carpet::prolongation_order_space = 3
+Carpet::prolongation_order_time = 2
+Carpet::init_each_timelevel = no
+Carpet::print_timestats_every = 0
+
+Time::timestep = 0.0001
+Time::timestep_method                   = "given"
+
+MethodOfLines::ode_method             = "generic"
+MethodOfLines::generic_type             = "RK"
+MethodOfLines::MoL_Intermediate_Steps = 1
+MethodOfLines::MoL_Num_Scratch_Levels = 1
+
+# MethodOfLines::ode_method             = "RK2"
+# MethodOfLines::MoL_NaN_Check          = "no"
+# MethodOfLines::initial_data_is_crap   = "yes"
+# MethodOfLines::MoL_Intermediate_Steps = 2
+# MethodOfLines::MoL_Num_Scratch_Levels = 1
+
+IO::out_dir = $parfile
+
+IOBasic::outInfo_every                  = 100
+IOBasic::outInfo_reductions             = "minimum maximum"
+IOBasic::outInfo_vars                   = "
+  Carpet::physical_time_per_hour
+  Laplace::phi
+  Laplace::phirhs
+"
+
+IOScalar::outScalar_every               = 10
+IOScalar::outScalar_reductions          = "norm2"
+IOScalar::outScalar_vars                = "Laplace::phirhs"
+
+Laplace::Lx = 1
+Laplace::Ly = 1
+Laplace::phi0 = 1
+Laplace::mu = 1
+
+IOHDF5::out_every                       = 100
+IOHDF5::out_vars                        = "
+  Laplace::phi
+  Laplace::phirhs
+"
+
+# TimerReport::out_every = 0
+# TimerReport::n_top_timers = 50
diff --git a/Tools/CodeGen/CalculationBoundaries.m b/Tools/CodeGen/CalculationBoundaries.m
index 3d331ab..3a66aba 100644
--- a/Tools/CodeGen/CalculationBoundaries.m
+++ b/Tools/CodeGen/CalculationBoundaries.m
@@ -19,7 +19,7 @@
 *)
 
 BeginPackage["CalculationBoundaries`", {"Errors`", "Helpers`",
-  "Kranc`", "MapLookup`", "KrancGroups`", "CodeGen`"}];
+  "Kranc`", "MapLookup`", "KrancGroups`", "CodeGen`", "CodeGenCactus`", "CodeGenC`"}];
 
 CalculationBoundariesFunction;
 
diff --git a/Tools/CodeGen/CalculationFunction.m b/Tools/CodeGen/CalculationFunction.m
index 35e9625..1fd5894 100644
--- a/Tools/CodeGen/CalculationFunction.m
+++ b/Tools/CodeGen/CalculationFunction.m
@@ -19,9 +19,9 @@
     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *)
 
-BeginPackage["CalculationFunction`", {"CodeGen`",
+BeginPackage["CalculationFunction`", {"CodeGenCactus`", "CodeGenC`", "CodeGen`",
   "MapLookup`", "KrancGroups`", "Differencing`", "Errors`",
-  "Helpers`", "Kranc`", "Optimize`", "Jacobian`"}];
+  "Helpers`", "Kranc`", "Optimize`", "Jacobian`", "Profile`"}];
 
 CreateCalculationFunction::usage = "";
 VerifyCalculation::usage = "";
@@ -69,7 +69,7 @@ VerifyListContent[l_, type_, while_] :=
   Module[{types},
     If[!(Head[l] === List),
       ThrowError["Expecting a list of ", type,
-        " objects, but found the following, which is not a List object: ", l, while]];
+        " objects, but found the following, which is not a List object.  Error occured ", while, l]];
     types = Union[Map[Head, l]];
     If [!(types === {type}) && !(types === {}),
       ThrowError["Expecting a list of ", type ,
@@ -241,7 +241,7 @@ simpCollect[collectList_, eqrhs_, localvar_, debug_] :=
     collectCoeff = Collect[rhs, localCollectList];
     InfoMessage[InfoFull, "ByteCount[terms collected]: ", ByteCount@collectCoeff];
 
-    all = Collect[rhs, localCollectList, Simplify];
+    all = Profile["Collect/Simplify", Collect[rhs, localCollectList, Simplify]];
     InfoMessage[InfoFull, "ByteCount[simplified rhs]: ", ByteCount@all];
 
     all];
@@ -373,11 +373,13 @@ pdCanonicalOrdering[name_[inds___] -> x_] :=
 
 Options[CreateCalculationFunction] = ThornOptions;
 
-CreateCalculationFunction[calcp_, debug_, imp_, opts:OptionsPattern[]] :=
+DefFn[
+  CreateCalculationFunction[calcp_, debug_, imp_, opts:OptionsPattern[]] :=
   Module[{gfs, allSymbols, knownSymbols,
           shorts, eqs, parameters, parameterRules,
           functionName, dsUsed, groups, pddefs, cleancalc, eqLoop, where,
-          addToStencilWidth, pDefs, haveCondTextuals, condTextuals, calc},
+          addToStencilWidth, pDefs, haveCondTextuals, condTextuals, calc,
+          kernelCall},
 
   calc = If[OptionValue[UseJacobian], InsertJacobian[calcp, opts], calcp];
 
@@ -451,7 +453,7 @@ CreateCalculationFunction[calcp_, debug_, imp_, opts:OptionsPattern[]] :=
   (* Check that there are no unknown symbols in the calculation *)
   allSymbols = calculationSymbols[cleancalc];
   knownSymbols = Join[lookupDefault[cleancalc, AllowedSymbols, {}], gfs, shorts, parameters,
-    {dx,dy,dz,idx,idy,idz,t, Pi, E, Symbol["i"], Symbol["j"], Symbol["k"], normal1, normal2,
+    {dx,dy,dz,dt,idx,idy,idz,t, Pi, E, Symbol["i"], Symbol["j"], Symbol["k"], normal1, normal2,
     normal3, tangentA1, tangentA2, tangentA3, tangentB1, tangentB2, tangentB3},
     If[OptionValue[UseJacobian], JacobianSymbols[], {}]];
 
@@ -461,26 +463,41 @@ CreateCalculationFunction[calcp_, debug_, imp_, opts:OptionsPattern[]] :=
      ThrowError["Unknown symbols in calculation.  Symbols are:", unknownSymbols,
        "Calculation is:", cleancalc]];
 
+  kernelCall = Switch[where,
+    Everywhere,
+      "GenericFD_LoopOverEverything(cctkGH, &" <> bodyFunctionName <> ");\n",
+    Interior,
+      "GenericFD_LoopOverInterior(cctkGH, &" <> bodyFunctionName <> ");\n",
+    InteriorNoSync,
+      "GenericFD_LoopOverInterior(cctkGH, &" <> bodyFunctionName <> ");\n",
+    Boundary,
+      "GenericFD_LoopOverBoundary(cctkGH, &" <> bodyFunctionName <> ");\n",
+    BoundaryWithGhosts,
+      "GenericFD_LoopOverBoundaryWithGhosts(cctkGH, &" <> bodyFunctionName <> ");\n",
+    _,
+      ThrowError["Unknown 'Where' entry in calculation " <>
+        functionName <> ": " <> ToString[where]]];
+
   {
-  DefineFunction[bodyFunctionName, "static void", "cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]",
+  DefineFunction[bodyFunctionName, "static void", "cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]",
   {
     "DECLARE_CCTK_ARGUMENTS;\n",
     "DECLARE_CCTK_PARAMETERS;\n\n",
-    If[!OptionValue[UseLoopControl], DeclareGridLoopVariables[], {}],
-    DeclareFDVariables[],
-
-    ConditionalOnParameterTextual["verbose > 1",
-      "CCTK_VInfo(CCTK_THORNSTRING,\"Entering " <> bodyFunctionName <> "\");\n"],
 
-    ConditionalOnParameterTextual["cctk_iteration % " <> functionName <> "_calc_every != " <>
-      functionName <> "_calc_offset", "return;\n"],
+    (* OpenCL kernel prologue *)
+    (* We could (or probably should) write this into a source file of its own *)
+    If[OptionValue[UseOpenCL],
+       {
+         "char const * const source =\n"
+       },
+       {
+       }],
 
-    CheckGroupStorage[groupsInCalculation[cleancalc, imp], functionName],
+    If[OptionValue[UseOpenCL], Stringify, Identity][{
 
-    "\n",
-    CheckStencil[pddefs, eqs, functionName, lookup[{opts}, IntParameters, {}]],
+    If[!OptionValue[UseLoopControl], DeclareGridLoopVariables[], {}],
 
-    If[haveCondTextuals, Map[ConditionalOnParameterTextual["!(" <> # <> ")", "return;\n"] &,condTextuals], {}],
+    DeclareFDVariables[],
 
     CommentedBlock["Include user-supplied include files",
       Map[IncludeFile, lookupDefault[cleancalc, DeclarationIncludes, {}]]],
@@ -499,38 +516,65 @@ CreateCalculationFunction[calcp_, debug_, imp_, opts:OptionsPattern[]] :=
       {
 	eqLoop = equationLoop[eqs, cleancalc, gfs, shorts, {}, groups,
           pddefs, where, addToStencilWidth, opts]},
-
-       ConditionalOnParameterTextual["verbose > 1",
-         "CCTK_VInfo(CCTK_THORNSTRING,\"Leaving " <> bodyFunctionName <> "\");\n"],
       {}]
+
     }],
-  Switch[where,
-    Everywhere,
-      DefineCCTKSubroutine[functionName,
-        "GenericFD_LoopOverEverything(cctkGH, &" <> bodyFunctionName <> ");\n"],
-    Interior,
-      DefineCCTKSubroutine[functionName,
-        "GenericFD_LoopOverInterior(cctkGH, &" <> bodyFunctionName <> ");\n"],
-    InteriorNoSync,
-      DefineCCTKSubroutine[functionName,
-        "GenericFD_LoopOverInterior(cctkGH, &" <> bodyFunctionName <> ");\n"],
-    Boundary,
-      DefineCCTKSubroutine[functionName,
-        "GenericFD_LoopOverBoundary(cctkGH, &" <> bodyFunctionName <> ");\n"],
-    BoundaryWithGhosts,
-      DefineCCTKSubroutine[functionName,
-        "GenericFD_LoopOverBoundaryWithGhosts(cctkGH, &" <> bodyFunctionName <> ");\n"],
-    PenaltyPrim2Char,
-      DefineFunction[functionName, "CCTK_INT",
-               "CCTK_POINTER_TO_CONST const cctkGH_, CCTK_INT const dir, CCTK_INT const face, CCTK_REAL const * restrict const base_, CCTK_INT const * restrict const lbnd, CCTK_INT const * restrict const lsh, CCTK_INT const * restrict const from, CCTK_INT const * restrict const to, CCTK_INT const rhs_flag, CCTK_INT const num_modes, CCTK_POINTER const * restrict const modes, CCTK_POINTER const * restrict const speeds",
-        "GenericFD_PenaltyPrim2Char(cctkGH, &" <> bodyFunctionName <> ");\n"],
-    _,
-      ThrowError["Unknown 'Where' entry in calculation " <>
-        functionName <> ": " <> ToString[where]]]}];
+
+    (* OpenCL kernel epologue *)
+    If[OptionValue[UseOpenCL],
+       {
+         ";\n\n",
+         Module[
+           {ignoreGroups, groupsNames, groupNameList},
+           ignoreGroups = {"TmunuBase::stress_energy_scalar",
+                           "TmunuBase::stress_energy_vector",
+                           "TmunuBase::stress_energy_tensor"};
+           groupNames = groupsInCalculation[cleancalc, imp];
+           groupNames = Select[groupNames, !MemberQ[ignoreGroups, #] &];
+           {
+             "char const * const groups[] = {",
+               Riffle[Join[Map[Quote, groupNames], {"NULL"}], ","],
+               "};\n\n"
+           }
+         ],
+         "static struct OpenCLKernel * kernel = NULL;\n",
+         "char const * const sources[] = {differencing, source, NULL};\n",
+         "OpenCLRunTime_CallKernel (cctkGH, CCTK_THORNSTRING, \"" <> functionName <> "\",\n",
+         "                          sources, groups, NULL, NULL, NULL, -1,\n",
+         "                          imin, imax, &kernel);\n\n"
+       },
+       {
+       }]
+    }],
+  
+    DefineCCTKSubroutine[functionName,
+      FlattenBlock[{
+        ConditionalOnParameterTextual["verbose > 1",
+          "CCTK_VInfo(CCTK_THORNSTRING,\"Entering " <> bodyFunctionName <> "\");\n"],
+
+        ConditionalOnParameterTextual["cctk_iteration % " <> functionName <> "_calc_every != " <>
+          functionName <> "_calc_offset", "return;\n"],
+  
+        CheckGroupStorage[groupsInCalculation[cleancalc, imp], functionName],
+        "\n",
+
+        CheckStencil[pddefs, eqs, functionName, OptionValue[ZeroDimensions],
+                     lookup[{opts}, IntParameters, {}]],
+        "\n",
+  
+        If[haveCondTextuals, Map[ConditionalOnParameterTextual["!(" <> # <> ")", "return;\n"] &,condTextuals], {}],
+
+        kernelCall,
+
+        ConditionalOnParameterTextual["verbose > 1",
+          "CCTK_VInfo(CCTK_THORNSTRING,\"Leaving " <> bodyFunctionName <> "\");\n"]
+      }]]
+  }]];
 
 Options[equationLoop] = ThornOptions;
 
-equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_,
+DefFn[
+  equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_,
              where_, addToStencilWidth_,
              opts:OptionsPattern[]] :=
   Module[{rhss, lhss, gfsInRHS, gfsInLHS, gfsOnlyInRHS, localGFs,
@@ -563,10 +607,12 @@ equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_,
     localMap = Map[# -> localName[#] &, gfs];
 
     derivSwitch =
-      GridFunctionDerivativesInExpression[pddefs, eqsOrdered] != {};
+      GridFunctionDerivativesInExpression[pddefs, eqsOrdered,
+                                          OptionValue[ZeroDimensions]] != {};
 
     gfsDifferentiated = Map[First,
-      GridFunctionDerivativesInExpression[pddefs, eqsOrdered]];
+      GridFunctionDerivativesInExpression[pddefs, eqsOrdered,
+                                          OptionValue[ZeroDimensions]]];
 
     gfsDifferentiatedAndOnLHS = Intersection[gfsDifferentiated, gfsInLHS];
 
@@ -579,8 +625,9 @@ equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_,
 
     (* Replace the partial derivatives *)
     {defsWithoutShorts, defsWithShorts} = splitPDDefsWithShorthands[pddefs, shorts];
-    eqs2 = ReplaceDerivatives[defsWithoutShorts, eqsOrdered, True];
-    eqs2 = ReplaceDerivatives[defsWithShorts, eqs2, False];
+    eqs2 = ReplaceDerivatives[defsWithoutShorts, eqsOrdered, True,
+                              OptionValue[ZeroDimensions]];
+    eqs2 = ReplaceDerivatives[defsWithShorts, eqs2, False, OptionValue[ZeroDimensions]];
 
     checkEquationAssignmentOrder[eqs2, shorts];
     functionName = ToString@lookup[cleancalc, Name];
@@ -673,7 +720,9 @@ equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_,
         Map[IncludeFile, incs]],
 
       CommentedBlock["Precompute derivatives",
-        PrecomputeDerivatives[defsWithoutShorts, eqsOrdered, lookup[{opts}, IntParameters, {}]]],
+        PrecomputeDerivatives[defsWithoutShorts, eqsOrdered, 
+                              lookup[{opts}, IntParameters, {}],
+                              OptionValue[ZeroDimensions]]],
 
       CommentedBlock["Calculate temporaries and grid functions", calcCode],
 
@@ -681,37 +730,40 @@ equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_,
         Map[InfoVariable[#[[1]]] &, (eqs2 /. localMap)],
         ""],
 
-      If[OptionValue[UseVectors], {
-        CommentedBlock["If necessary, store only partial vectors after the first iteration",
-          ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 2 && CCTK_BUILTIN_EXPECT(i < lc_imin && i+CCTK_REAL_VEC_SIZE > lc_imax, 0)",
-            {
-              DeclareAssignVariable["ptrdiff_t", "elt_count_lo", "lc_imin-i"],
-              DeclareAssignVariable["ptrdiff_t", "elt_count_hi", "lc_imax-i"],
-              Map[StoreMiddlePartialVariableInLoop[GridName[#], localName[#], "elt_count_lo", "elt_count_hi"] &,
-                  gfsInLHS],
-              "break;\n"
-            }]],
-        CommentedBlock["If necessary, store only partial vectors after the first iteration",
-          ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 1 && CCTK_BUILTIN_EXPECT(i < lc_imin, 0)",
-            {
-              DeclareAssignVariable["ptrdiff_t", "elt_count", "lc_imin-i"],
-              Map[StoreHighPartialVariableInLoop[GridName[#], localName[#], "elt_count"] &,
-                  gfsInLHS],
-              "continue;\n"
-            }]],
-        CommentedBlock["If necessary, store only partial vectors after the last iteration",
-          ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 1 && CCTK_BUILTIN_EXPECT(i+CCTK_REAL_VEC_SIZE > lc_imax, 0)",
-            {
-              DeclareAssignVariable["ptrdiff_t", "elt_count", "lc_imax-i"],
-              Map[StoreLowPartialVariableInLoop[GridName[#], localName[#], "elt_count"] &,
-                  gfsInLHS],
-              "break;\n"
-            }]]}, {}],
-      CommentedBlock["Copy local copies back to grid functions",
-        Map[(If[OptionValue[UseVectors], StoreVariableInLoop, AssignVariableInLoop][GridName[#], localName[#]]) &,
-            gfsInLHS]],
-
-      If[debugInLoop, Map[InfoVariable[GridName[#]] &, gfsInLHS], ""]}, opts]];
+      Which[OptionValue[UseOpenCL],
+              CommentedBlock["Copy local copies back to grid functions",
+                Map[StorePartialVariableInLoop[GridName[#], localName[#]] &, gfsInLHS]],
+            OptionValue[UseVectors],
+              {
+                CommentedBlock["If necessary, store only partial vectors after the first iteration",
+                  ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 2 && CCTK_BUILTIN_EXPECT(i < lc_imin && i+CCTK_REAL_VEC_SIZE > lc_imax, 0)",
+                  {
+                    DeclareAssignVariable["ptrdiff_t", "elt_count_lo", "lc_imin-i"],
+                    DeclareAssignVariable["ptrdiff_t", "elt_count_hi", "lc_imax-i"],
+                    Map[StoreMiddlePartialVariableInLoop[GridName[#], localName[#], "elt_count_lo", "elt_count_hi"] &, gfsInLHS],
+                    "break;\n"
+                  }]],
+                CommentedBlock["If necessary, store only partial vectors after the first iteration",
+                  ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 1 && CCTK_BUILTIN_EXPECT(i < lc_imin, 0)",
+                  {
+                    DeclareAssignVariable["ptrdiff_t", "elt_count", "lc_imin-i"],
+                    Map[StoreHighPartialVariableInLoop[GridName[#], localName[#], "elt_count"] &, gfsInLHS],
+                    "continue;\n"
+                  }]],
+                CommentedBlock["If necessary, store only partial vectors after the last iteration",
+                  ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 1 && CCTK_BUILTIN_EXPECT(i+CCTK_REAL_VEC_SIZE > lc_imax, 0)",
+                  {
+                    DeclareAssignVariable["ptrdiff_t", "elt_count", "lc_imax-i"],
+                    Map[StoreLowPartialVariableInLoop[GridName[#], localName[#], "elt_count"] &, gfsInLHS],
+                    "break;\n"
+                  }]],
+                Map[StoreVariableInLoop[GridName[#], localName[#]] &, gfsInLHS]
+              },
+            True,
+              CommentedBlock["Copy local copies back to grid functions",
+                Map[AssignVariableInLoop[GridName[#], localName[#]] &, gfsInLHS]]],
+
+      If[debugInLoop, Map[InfoVariable[GridName[#]] &, gfsInLHS], ""]}, opts]]];
 
 End[];
 
diff --git a/Tools/CodeGen/CodeGen.m b/Tools/CodeGen/CodeGen.m
index 127ee81..c2cb890 100644
--- a/Tools/CodeGen/CodeGen.m
+++ b/Tools/CodeGen/CodeGen.m
@@ -22,13 +22,6 @@
 
 BeginPackage["CodeGen`", {"Errors`", "Kranc`"}];
 
-SOURCELANGUAGE::usage = "global variable == \"C\" or \"Fortran\" determines language
-                        for code generation";
-SOURCESUFFIX::usage = "file suffix for source files";
-
-EOL::usage = "the end of line termination string";
-SetSourceLanguage::usage = "set the source language to  \"C\" or \"Fortran\"";
-
 FlattenBlock::usage = "FlattenBlock[block] converts 'block' to a string.";
 SeparatedBlock::usage = "SeparatedBlock[block] returns a version of 'block' with " <>
   "a newline before it.";
@@ -36,961 +29,133 @@ GenerateFile::usage = "GenerateFile[name, block] writes 'block' to a file of the
   "specified 'name'.";
 AddToFile::usage = "AddToFile[name, block] appends 'block' to a file of the " <>
     "specified 'name'.";
-IncludeFile::usage = "IncludeFile[name] returns a block of code" <>
-  "that includes a header file (i.e '#include \"name\"').";
-IncludeSystemFile::usage = "IncludeFile[name] returns a block of code" <>
-  "that includes a system header file (i.e '#include <name>').";
-DeclareVariable::usage = "DeclareVariable[name, type] returns a block of code " <>
-  "that declares a variable of given name and type.  'name' and 'type' should be " <>
-  "strings.";
-DeclareVariableNoInit::usage = "DeclareVariableNoInit[name, type] returns a block of code " <>
-  "that declares a variable of given name and type without initialising it.  'name' and 'type' should be " <>
-  "strings.";
-DeclareVariables::usage = "DeclareVariables[names, type] returns a block of code " <>
-  "that declares a list of variables of given name and type.  'names' should be a list" <>
-  " of strings and 'type' should be a string string.";
-DeclarePointer::usage = "DeclarePointer[name, type] returns a block of code " <>
-  "that declares a pointer of given name and type.  'name' and 'type' should be " <>
-  "strings.";
-DeclarePointers::usage = "DeclarePointers[names, type] returns a block of code " <>
-  "that declares a list of pointers of given name and type.  'names' should be a list" <>
-  " of strings and 'type' should be a string string.";
-DefineVariable::usage = "DefineVariable[name, type, value] returns a block of " <>
-  "code that declares and initialised a variable 'name' of type 'type' to value 'value'.";
-AssignVariable::usage = "AssignVariable[dest_, src_] returns a block of code " <>
-  "that assigns 'src' to 'dest'.";
-DeclareAssignVariable::usage = "DeclareAssignVariable[type_, dest_, src_] returns a block of code " <>
-  "that declares and sets a constant variable of given name and type.";
-AssignVariableInLoop::usage = "AssignVariableInLoop[dest_, src_] returns a block of code " <>
-  "that assigns 'src' to 'dest'.";
-StoreVariableInLoop::usage = "StoreVariableInLoop[dest_, src_] returns a block of code " <>
-  "that assigns 'src' to 'dest'.";
-StoreLowPartialVariableInLoop::usage = "StoreLowPartialVariableInLoop[dest_, src_, count_] returns a block of code " <>
-  "that assigns 'src' to 'dest'.";
-StoreHighPartialVariableInLoop::usage = "StoreHighPartialVariableInLoop[dest_, src_, count_] returns a block of code " <>
-  "that assigns 'src' to 'dest'.";
-StoreMiddlePartialVariableInLoop::usage = "StoreMiddlePartialVariableInLoop[dest_, src_, countLow_, countHigh_] returns a block of code " <>
-  "that assigns 'src' to 'dest'.";
-DeclareAssignVariableInLoop::usage = "DeclareAssignVariableInLoop[type_, dest_, src_] returns a block of code " <>
-  "that assigns 'src' to 'dest'.";
-MaybeAssignVariableInLoop::usage = "MaybeAssignVariableInLoop[dest_, src_, cond_] returns a block of code " <>
-  "that assigns 'src' to 'dest'.";
-DeclareMaybeAssignVariableInLoop::usage = "DeclareMaybeAssignVariableInLoop[type_, dest_, src_, cond_] returns a block of code " <>
-  "that assigns 'src' to 'dest'.";
-TestForNaN::usage = "TestForNaN[expr_] returns a block of code " <>
-  "that tests 'expr' for nan.";
-CommentedBlock::usage = "CommentedBlock[comment, block] returns a block consisting " <>
-  "of 'comment' followed by 'block'.";
-DefineCCTKFunction::usage = "DefineCCTKFunction[name, type, block] returns a block " <>
-  "of code that defines a CCTK function of name 'name' returning type 'type' with " <>
-  "body 'block'.";
-DefineCCTKSubroutine::usage = "DefineCCTKSubroutine[name, block] returns a block " <>
-  "of code that defines a CCTK Fortran subroutine of name 'name' with body 'block'.";
-GridName::usage = "GridName[variable] returns the name needed to access variable " <>
-  "assuming it is a grid variable when inside a grid loop.";
-DeclareGridLoopVariables::usage = "DeclareGridLoopVariables[] returns a block " <>
-  "that defines the variables needed during a grid loop.";
-InitialiseGridLoopVariables::usage = "InitialiseGridLoopVariables[] returns a block " <>
-  "that initialises variables needed by a grid loop.";
-InitialiseGridLoopVariablesWithStencil::usage = "InitialiseGridLoopVariables[] returns a block " <>
-  "that initialises variables needed by a grid loop using the evolution stencils.";
-ConditionalOnParameter::usage = "ConditionalOnParameter[name, value, block] returns " <>
-  "a block that introduces a conditional expression whereby 'block' is only executed " <>
-  "if the Cactus parameter 'name' has value 'value'.";
-(*
-GridLoop::usage = "GridLoop[block] returns a block that is looped over for every " <>
-  "grid point.  Must have previously set up the grid loop variables (see " <>
-  "DeclareGridLoopVariables and InitialiseGridLoopVariables.";
-*)
-GridLoop::usage = "GridLoop[block] returns a block that is looped over for every " <>
-  "grid point.  Must have previously set up the grid loop variables (see " <>
-  "InitialiseGridLoopVariables.";
-
-DeclareArray::usage = "";
-
-DefineFunction::usage = "";
-ConditionalOnParameterTextual::usage = "";
-
-SpaceSeparated::usage = "";   (* This should not really be in CodeGen *)
-NewlineSeparated::usage = ""; (* This should not really be in CodeGen *)
-
-CBlock::usage = "";
-SuffixedCBlock::usage = "";
+SpaceSeparated::usage = "";
+NewlineSeparated::usage = "";
 InfoVariable::usage = "";
-
-DeclareFDVariables::usage = "";
-
-InitialiseFDVariables::usage = "";
-
-CommaNewlineSeparated::usage = ""; (* This should not really be in CodeGen *)
+CommaNewlineSeparated::usage = "";
 CommaSeparated::usage = "";
-ReplacePowers::usage = "";
-CFormHideStrings::usage = "";
-BoundaryLoop::usage = "";
-BoundaryWithGhostsLoop::usage = "";
-GenericGridLoop::usage = "";
-
-NameRoot::usage = "";
-PartitionVarList::usage = "";
+Stringify::usage = "";
 Quote::usage = "Quote[x] returns x surrounded by quotes";
-DataType::usage = "DataType[] returns a string for the grid function data type (e.g. CCTK_REAL)";
-SetDataType::usage = "SetDataType[type] sets a string for the grid function data type (e.g. CCTK_REAL)";
-Conditional;
-SwitchStatement;
-
-Begin["`Private`"];
-
-SOURCELANGUAGE =  "C";
-SOURCESUFFIX   = ".cc";
-
-setSourceSuffix[lang_] :=
-If[ (lang == "C"),
-  SOURCESUFFIX = ".cc";
-,
-  SOURCESUFFIX = ".F90";
-];
-
-
-SetSourceLanguage[lang_] :=
-If[ (lang == "C" || lang == "Fortran"),
-  SOURCELANGUAGE = lang;
-  setSourceSuffix[lang];
-  InfoMessage[Terse, "User set source language to " <> lang],
+IndentBlock::usage = "";
+CheckBlock::usage = "";
 
-  SOURCELANGUAGE = "C";
-  setSourceSuffix[".cc"];
-  InfoMessage[Terse, "Setting Source Language to C"];
-];
+CodeGenBlock := _String | _?AtomQ | List[(_?(MatchQ[#, CodeGenBlock] &)) ...];
+Boolean = (True | False);
 
-SetDataType[type_String] :=
-  dataType = type;
-
-DataType[] :=
-  If[dataType === Symbol["datatype"],
-    Throw["DataType: Have not set a data type"],
-    dataType];
+Begin["`Private`"];
 
 (* Code generation utilities; not specific to any language *)
 
-FlattenBlock[b_] := Apply[StringJoin,Map[ToString,If[! AtomQ[b], Flatten[b, Infinity], b]]];
-
-indentBlock[block_] :=
-  StringDrop["  " <> StringReplace[FlattenBlock[block], {"\n" -> "\n  "}],-2];
-
-SeparatedBlock[block_] := {"\n", block};
-
-GenerateFile[filename_, contents_] :=
-  Module[{fp = OpenWrite[filename]},
-    WriteString[fp, FlattenBlock[contents]];
-    Close[fp]];
-
-AddToFile[filename_, contents_] :=
-  Module[{fp = OpenAppend[filename]},
-      WriteString[fp, FlattenBlock[contents]];
-          Close[fp]];
-
-intersperse[l_, x_] :=
-  If[l == {},
-    {},
-    If[Rest[l] == {},
-      {l[[1]]},
-      Join[{l[[1]]}, {x}, intersperse[Rest[l],x]]]];
-
-CommaNewlineSeparated[l_] := intersperse[l, ",\n"];
-
-SpaceSeparated[l_] := 
-  Module[{},
-  If[!ListQ[l],
-    ThrowError["SpaceSeparated: Expecting a list, but was given", l]];
-  intersperse[l, " "]];
-
-CommaSeparated[l_] := 
-  intersperse[l, ", "];
-
-NewlineSeparated[l_] := 
-  intersperse[l, "\n"];
-
-CommaInitSeparated[l_] :=
-  intersperse[Map[{#," = INITVALUE"} &, l], ", "];
-(*  intersperse[l, " = INITVALUE, "];*)
-
-
-
-NameRoot[name_] := Module[{dropNumberRule, root},
+DefFn[
+  CheckBlock[s_String] := s];
 
-      dropNumberRule = {"1" -> "", "2" -> "", "3" -> "", "4" -> "", "5" -> "",
-                        "6" -> "", "7" -> "", "8" -> "", "9" -> "", "0" -> "", "rhs" -> ""};
+DefFn[
+  CheckBlock[a_?AtomQ] := a];
 
-      root = StringReplace[ToString@name, dropNumberRule]
-      ];
+DefFn[
+  CheckBlock[l_List] := Map[CheckBlock, l]];
 
-PartitionVarList[list_]:= Module[{partition, split},
+DefFn[
+  FlattenBlock[b_] :=
+  Module[
+    {flattenBlock},
+    flattenBlock[x_String] := x;
+    flattenBlock[l_List] := StringJoin@@Map[FlattenBlock, l];
+    flattenBlock[a_?AtomQ] := ToString[a];
 
-partition[locallist_] := Module[{cutoff},
-  cutoff = 6;
-  If[Length@locallist > cutoff, Partition[locallist, cutoff, cutoff, {1,1}, {}], {locallist}]
-];
+    CheckBlock[b];
+    flattenBlock[b]]];
 
-split = Split[list, NameRoot[#1] == NameRoot[#2] &];
-split = Flatten[Map[partition, split], 1];
+DefFn[
+  IndentBlock[block:CodeGenBlock] :=
+  StringDrop["  " <> StringReplace[FlattenBlock[block], {"\n" -> "\n  "}],-2]];
 
-split
-];
+DefFn[
+  SeparatedBlock[block:CodeGenBlock] := {"\n", block}];
+ErrorDefinition[SeparatedBlock];
 
+DefFn[
+  GenerateFile[filename_String, contents_] :=
+  Module[
+    {fp = OpenWrite[filename]},
+    CheckBlock[contents];
+    WriteString[fp, FlattenBlock[contents]];
+    Close[fp]]];
 
-(* Code generation for generic C and C-preprocessed Fortran *)
-
-EOL[dummy___] := If[SOURCELANGUAGE == "C" || SOURCELANGUAGE == "C++", ";\n", "\n"];
-
-IncludeFile[filename_] :=
-  {"#include \"", filename, "\"\n"};
-
-IncludeSystemFile[filename_] :=
-  {"#include <", filename, ">\n"};
-
-DeclareVariable[name_, type_] :=
-If[SOURCELANGUAGE == "C",
-  {type, " ",    name, " = INITVALUE" <> EOL[]},
-  {type, " :: ", name, EOL[]} (* no value init here to avoid implicit SAVE attribute *)
-  ];
-
-DeclareVariableNoInit[name_, type_] :=
-If[SOURCELANGUAGE == "C",
-  {type, " ",    name, EOL[]},
-  {type, " :: ", name, EOL[]} (* no value init here to avoid implicit SAVE attribute *)
-  ];
-
-
-DeclareVariables[names_?ListQ, type_] := 
-If[SOURCELANGUAGE == "C",
-   {type, " ",    CommaSeparated@names, EOL[]},
-   {type, " :: ", CommaSeparated@names,     EOL[]} (* no value init avoids implicit SAVE attribute *)
-   ];
-
-DeclarePointer[name_, type_] :=
-If[SOURCELANGUAGE == "C",
-  {type, " *",    name, EOL[]},
-  {type, ", target :: ", name, EOL[]}
-  ];
-
-DeclarePointers[names_?ListQ, type_] :=
-If[SOURCELANGUAGE == "C",
-   {type, " *",           CommaInitSeparated@names, EOL[]},
-   {type, ", target :: ", CommaSeparated@names,     EOL[]} 
-   ];
-
-DeclareArray[name_, dim_, type_] :=
-  If[SOURCELANGUAGE == "C",
-     DeclareArrayC[name, dim, type],
-     DeclareArrayFortran[name, dim, type]];
-
-DeclareArrayC[name_, dim_, type_] :=
-  {type, " ", name, "[", dim, "];","\n"};
-
-DeclareArrayFortran[name_, dim_, type_] :=
-  {type, " :: ", name, "(", dim, ")","\n"};
-
-DefineVariable[name_, type_, value_] :=
-  {type, " ", name, " = ", value, EOL[]};
-
-AssignVariable[dest_, src_] :=
-  {dest, " = ", src, EOL[]};
-
-DeclareAssignVariable[type_, dest_, src_] :=
-  {type, " const ", dest, " = ", src, EOL[]};
-
-AssignVariableInLoop[dest_, src_, vectorise_:False] :=
-  Module[{loader},
-    loader[x_] := If[vectorise, {"vec_load(", x, ")"}, x];
-    {dest, " = ", loader[src], EOL[]}];
-(*
-  {dest, " = ", src, EOL[],
-   TestForNaN[dest]};
-*)
-
-StoreVariableInLoop[dest_, src_] :=
-  {"vec_store_nta(", dest, ",", src, ")", EOL[]};
-
-StoreLowPartialVariableInLoop[dest_, src_, count_] :=
-  {"vec_store_nta_partial_lo(", dest, ",", src, ",", count, ")", EOL[]};
-
-StoreHighPartialVariableInLoop[dest_, src_, count_] :=
-  {"vec_store_nta_partial_hi(", dest, ",", src, ",", count, ")", EOL[]};
-
-StoreMiddlePartialVariableInLoop[dest_, src_, countLow_, countHigh_] :=
-  {"vec_store_nta_partial_mid(", dest, ",", src, ",", countLow, ",", countHigh, ")", EOL[]};
-
-DeclareAssignVariableInLoop[type_, dest_, src_] :=
-  {type, " const ", dest, " = vec_load(", src, ")", EOL[]};
-
-MaybeAssignVariableInLoop[dest_, src_, cond_] :=
-  If [cond,
-      {dest, " = useMatter ? vec_load(", src, ") : ToReal(0.0)", EOL[]},
-      {dest, " = vec_load(", src, ")", EOL[]}];
-
-DeclareMaybeAssignVariableInLoop[type_, dest_, src_, mmaCond_, codeCond_, vectorise_:False] :=
-  Module[{loader},
-    loader[x_] := If[vectorise, {"vec_load(", x, ")"}, x];
-    If [mmaCond,
-        {type, " ", dest, " = (", codeCond, ") ? ", loader[src], " : ToReal(0.0)", EOL[]},
-        {type, " ", dest, " = ", loader[src], EOL[]}]];
-
-TestForNaN[expr_] :=
-  {"if (isnan(", expr, ")) {\n",
-   "  CCTK_VInfo(CCTK_THORNSTRING, \"NaN found\");\n",
-   "  CCTK_VInfo(CCTK_THORNSTRING, \"ipos: %d %d %d\", i, j, k);\n",
-   "  CCTK_VInfo(CCTK_THORNSTRING, \"lbnd: %d %d %d\", cctk_lbnd[0], cctk_lbnd[1], cctk_lbnd[2]);\n",
-   "  CCTK_VInfo(CCTK_THORNSTRING, \"lsh: %d %d %d\", cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]);\n",
-   "  CCTK_VInfo(CCTK_THORNSTRING, \"LSSH: %d %d %d\", CCTK_LSSH(0,0), CCTK_LSSH(0,1), CCTK_LSSH(0,2));\n",
-   "  CCTK_VInfo(CCTK_THORNSTRING, \"", expr, ": %.17g\", (double)", expr, ");\n",
-   "}\n"};
-
-(* comments are always done C-style because they are killed by cpp anyway *) 
-insertComment[text_] := {"/* ", text, " */\n"};
-
-CBlock[block_] :=
- {"{\n",
-  indentBlock[block],
-  "}\n"};
-
-SuffixedCBlock[block_, suffix_] :=
- {"{\n",
-  indentBlock[block],
-  "} ", suffix, "\n"};
-
-
-loopOverInteger[name_, start_, endplusone_, block_] :=
-If[SOURCELANGUAGE == "C" || SOURCELANGUAGE == "C++",
-
-  {"for (", name, " = ", start, "; ", name, " < ", endplusone, "; ", name, "++)\n",
-   "{\n",
-   indentBlock[block],
-   "}\n"},
-
-  {"Do ", name, " = ", start, ", ", endplusone, "\n",
-   "\n",
-   indentBlock[block],
-   "End Do\n"}
-];
-
-
-CommentedBlock[comment_, block_] :=
-  SeparatedBlock[{insertComment[comment],
-                  block}];
-
-(* FUNCTIONS *)
-
-defineFunctionC[name_, type_, args_, contents_] :=
-  SeparatedBlock[
-    {type, " ", name, "(", args, ")\n",
-      CBlock[contents]}];
-     
-defineFunctionF[name_, args_, contents_] :=
-  SeparatedBlock[
-    {"FUNCTION", " ", name, "(", args, ")\n",
-      indentBlock[contents]}];
-
-DefineFunction[name_, type_, args_, contents_] :=
-  If[SOURCELANGUAGE == "C",
-    defineFunctionC[name, type, args, contents],
-    defineFunctionF[name, args, contents]];
-
-(* SUBROUTINES *)
-
-defineSubroutine[name_, args_, contents_] :=
-  If[SOURCELANGUAGE == "C",
-    defineSubroutineC[name, args, contents],
-    defineSubroutineF[name, args, contents]];
-
-defineSubroutineC[name_, args_, contents_] :=
-  SeparatedBlock[
-    {"extern \"C\" void ", name, "(", args, ")", "\n",
-     CBlock[contents]}];
-
-defineSubroutineF[name_, args_, contents_] :=
-  SeparatedBlock[
-    {"subroutine ", name, "(", args, ")", "\n",
-     "\nimplicit none\n\n",
-     contents,
-     "end subroutine\n"}];
-
-
-
-
-(********* Code generation for Cactus C or Fortran code **********)
-
-
-
-
-(* This is a Cactus-callable function *)
-DefineCCTKFunction[name_, type_, contents_] :=
-  DefineFunction[name, "extern \"C\" " <> type, "CCTK_ARGUMENTS", 
-    {
-      "DECLARE_CCTK_ARGUMENTS;\n",
-      "DECLARE_CCTK_PARAMETERS;\n\n",
-      contents
-    }];
-
-(* This is a Cactus-callable subroutine *)
-DefineCCTKSubroutine[name_, contents_] :=
-  defineSubroutine[
-    name, "CCTK_ARGUMENTS", 
-    {
-      "DECLARE_CCTK_ARGUMENTS;\n",
-      "DECLARE_CCTK_PARAMETERS;\n\n",
-      contents
-    }];
-
-DeclareFDVariables[] := 
-(*
-  CommentedBlock["Declare finite differencing variables",
-    {Map[DeclareVariables[#, "CCTK_REAL"] &, {{"dx", "dy", "dz"}, 
-                                              {"dxi", "dyi", "dzi"},
-                                              {khalf,kthird,ktwothird,kfourthird,keightthird},
-                                              {"hdxi", "hdyi", "hdzi"}}],
-     "\n"},
-    {Map[DeclareVariables[#, "ptrdiff_t"] &, {{"di", "dj", "dk"}}],
-     "\n"}];
-*)
-  CommentedBlock["Declare finite differencing variables", {}];
-
-InitialiseFDSpacingVariablesC[] :=
-  {
-    (* DeclareAssignVariable["ptrdiff_t", "di", "CCTK_GFINDEX3D(cctkGH,1,0,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"], *)
-    DeclareAssignVariable["ptrdiff_t", "di", "1"],
-    DeclareAssignVariable["ptrdiff_t", "dj", "CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"],
-    DeclareAssignVariable["ptrdiff_t", "dk", "CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0)"],
-    DeclareAssignVariable["ptrdiff_t", "cdi", "sizeof(CCTK_REAL) * di"],
-    DeclareAssignVariable["ptrdiff_t", "cdj", "sizeof(CCTK_REAL) * dj"],
-    DeclareAssignVariable["ptrdiff_t", "cdk", "sizeof(CCTK_REAL) * dk"],
-    DeclareAssignVariable[DataType[], "dx", "ToReal(CCTK_DELTA_SPACE(0))"],
-    DeclareAssignVariable[DataType[], "dy", "ToReal(CCTK_DELTA_SPACE(1))"],
-    DeclareAssignVariable[DataType[], "dz", "ToReal(CCTK_DELTA_SPACE(2))"],
-    DeclareAssignVariable[DataType[], "dt", "ToReal(CCTK_DELTA_TIME)"]
-  };
-
-InitialiseFDSpacingVariablesFortran[] := 
-  {
-    AssignVariable["dt", "CCTK_DELTA_TIME"],
-    AssignVariable["dx", "CCTK_DELTA_SPACE(1)"],
-    AssignVariable["dy", "CCTK_DELTA_SPACE(2)"],
-    AssignVariable["dz", "CCTK_DELTA_SPACE(3)"]
-  }
-
-
-InitialiseFDVariables[vectorise_] :=
-  CommentedBlock["Initialise finite differencing variables",
-  { If[SOURCELANGUAGE == "Fortran",
-       InitialiseFDSpacingVariablesFortran[],
-       InitialiseFDSpacingVariablesC[]],
-    
-    DeclareAssignVariable[DataType[], "dxi", "INV(dx)"],
-    DeclareAssignVariable[DataType[], "dyi", "INV(dy)"],
-    DeclareAssignVariable[DataType[], "dzi", "INV(dz)"],
-    If[vectorise,
-     {DeclareAssignVariable[DataType[], "khalf", "ToReal(0.5)"],
-      DeclareAssignVariable[DataType[], "kthird", "ToReal(1.0/3.0)"],
-      DeclareAssignVariable[DataType[], "ktwothird", "ToReal(2.0/3.0)"],
-      DeclareAssignVariable[DataType[], "kfourthird", "ToReal(4.0/3.0)"],
-      DeclareAssignVariable[DataType[], "keightthird", "ToReal(8.0/3.0)"],
-      DeclareAssignVariable[DataType[], "hdxi", "kmul(ToReal(0.5), dxi)"],
-      DeclareAssignVariable[DataType[], "hdyi", "kmul(ToReal(0.5), dyi)"],
-      DeclareAssignVariable[DataType[], "hdzi", "kmul(ToReal(0.5), dzi)"]},
-     {DeclareAssignVariable[DataType[], "khalf", "0.5"],
-      DeclareAssignVariable[DataType[], "kthird", "1/3.0"],
-      DeclareAssignVariable[DataType[], "ktwothird", "2.0/3.0"],
-      DeclareAssignVariable[DataType[], "kfourthird", "4.0/3.0"],
-      DeclareAssignVariable[DataType[], "keightthird", "8.0/3.0"],
-      DeclareAssignVariable[DataType[], "hdxi", "0.5 * dxi"],
-      DeclareAssignVariable[DataType[], "hdyi", "0.5 * dyi"],
-      DeclareAssignVariable[DataType[], "hdzi", "0.5 * dzi"]}]}];
-
-GridName[x_] := If[SOURCELANGUAGE == "C",
-                   ToExpression[ToString[x] <> "[index]"],
-                   ToString[x] <> "(i,j,k)"
-                ];
-
-DeclareGridLoopVariables[] :=
-  SeparatedBlock[
-    {insertComment["Declare the variables used for looping over grid points"],
-     Map[DeclareVariables[#, "CCTK_INT"] &, 
-         {{"i", "j", "k"}(*, {"istart", "jstart", "kstart"}, 
-          {"iend", "jend", "kend"},
-          {"index_offset_x", "index_offset_y", "index_offset_z", "dir", "face"} *)}] (*,
-     Map[DeclareArray[#, 6, "CCTK_INT"] &, {"is_symbnd", "is_physbnd", "is_ipbnd"}],
-     Map[DeclareArray[#, 3, "CCTK_INT"] &, {"imin", "imax", "bmin", "bmax"}] *), 
-
-     If[SOURCELANGUAGE == "C", DeclareVariable["index", "// CCTK_INT"], "\n"]
-  }];
-
-(* Access an element of an array; syntax is different between C and
-   Fortran.  Always give this function a C-style array index. *)
-arrayElement[var_, i_] :=
-  If[SOURCELANGUAGE == "C",
-     {var, "[", arrayIndex[i], "]"},
-     {var, "(", arrayIndex[i], ")"}];
-
-(* Given a C-style variable index, return the corresponding index for
-   the language currently in use.  The idea is that the caller does not
-   need to know what language is being used. *)
-arrayIndex[i_] :=
-  If[SOURCELANGUAGE == "C",
-     i,
-     If[NumberQ[i], i+1, {i, " + 1"}]];
-
-max[]:= If[SOURCELANGUAGE == "C", "IMAX", "max"];
-
-InitialiseGridLoopVariables[derivativesUsedSwitch_, addToStencilWidth_] :=
-  CommentedBlock["Set up variables used in the grid loop for the physical grid points",
-
-  If[ (derivativesUsedSwitch),
-  {
-  AssignVariable["index_offset_x", max[] <>"(stencil_width, stencil_width_x) + " <> ToString[addToStencilWidth]],
-  AssignVariable["index_offset_y", max[] <>"(stencil_width, stencil_width_y) + " <> ToString[addToStencilWidth]],
-  AssignVariable["index_offset_z", max[] <>"(stencil_width, stencil_width_z) + " <> ToString[addToStencilWidth]],
-
-  "\n",
-  AssignVariable["istart", arrayIndex["index_offset_x"]],
-  AssignVariable["jstart", arrayIndex["index_offset_y"]],
-  AssignVariable["kstart", arrayIndex["index_offset_z"]],
-
-  "\n",
-  AssignVariable["iend", {"CCTK_LSSH(0,0) - index_offset_x"}],
-  AssignVariable["jend", {"CCTK_LSSH(0,1) - index_offset_y"}],
-  AssignVariable["kend", {"CCTK_LSSH(0,2) - index_offset_z"}]
-  },
-
-  {
-  AssignVariable["istart", arrayIndex[0]],
-  AssignVariable["jstart", arrayIndex[0]],
-  AssignVariable["kstart", arrayIndex[0]],
-
-  "\n",
-  AssignVariable["iend", "CCTK_LSSH(0,0)"],
-  AssignVariable["jend", "CCTK_LSSH(0,1)"],
-  AssignVariable["kend", "CCTK_LSSH(0,2)"]
-  }]
-];
-
-
-ConditionalOnParameter[name_, value_, block_] :=
-  SeparatedBlock[
-  {"if (CCTK_EQUALS(", name, ", \"", value, "\"))\n",
-   "{\n",
-   indentBlock[block],
-   "}\n"}];
-
-ConditionalOnParameterTextual[text_, block_] :=
-  SeparatedBlock[
-  {"if (", text, ")\n",
-   "{\n",
-   indentBlock[block],
-   "}\n"}];
-
-(*
-GridLoop[block_] :=
-  CommentedBlock["Loop over the grid points",
-   loopOverInteger["k", "kstart", "kend",
-     loopOverInteger["j", "jstart", "jend",
-       loopOverInteger["i", "istart", "iend",
-
-       { If[SOURCELANGUAGE == "C",  
-            AssignVariable["index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"],
-            ""],
-	 block
-       }
-        ]]]];
-
-*)
-
-(*
-GridLoop[block_] :=
-  If[SOURCELANGUAGE == "C",  
-    CommentedBlock["Loop over the grid points",
-      {
-        "#pragma omp parallel\n",
-        "LC_LOOP3 (unnamed,\n",
-        "          i,j,k, istart,jstart,kstart, iend,jend,kend,\n",
-        "          cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])\n",
-        "{\n",
-        indentBlock[
-          {
-             DeclareVariable["index", "int"],
-             AssignVariable["index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"],
-             block
-          }
-        ],
-        "}\n",
-        "LC_ENDLOOP3 (unnamed);\n"
-      }
-    ],
-    CommentedBlock["Loop over the grid points",
-      {
-        "#pragma omp parallel\n",
-        "LC_LOOP3 (unnamed,\n",
-        "          i,j,k, istart,jstart,kstart, iend,jend,kend,\n",
-        "          cctk_lsh(1),cctk_lsh(2),cctk_lsh(3))\n",
-        indentBlock[block],
-        "LC_ENDLOOP3 (unnamed)\n"
-      }
-    ]
-  ];
-*)
-
-Options[GenericGridLoop] = ThornOptions;
-
-GenericGridLoop[functionName_, block_, opts:OptionsPattern[]] :=
-  If[OptionValue[UseLoopControl],
-    GenericGridLoopUsingLoopControl[functionName, block, OptionValue[UseVectors]],
-    GenericGridLoopTraditional[block]];
-
-GenericGridLoopTraditional[block_] :=
-  CommentedBlock["Loop over the grid points",
-   loopOverInteger["k", "min[2]", "max[2]",
-     loopOverInteger["j", "min[1]", "max[1]",
-       loopOverInteger["i", "min[0]", "max[0]",
-
-       { If[SOURCELANGUAGE == "C",  
-            {
-              DeclareAssignVariable["int", "index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"]
-            }
-            ""],
-	 block
-       }
-        ]]]];
-
-GenericGridLoopUsingLoopControl[functionName_, block_, vectorise_] :=
-  If[SOURCELANGUAGE == "C",  
-    CommentedBlock["Loop over the grid points",
-      {
-        "#pragma omp parallel\n",
-        If[vectorise, "LC_LOOP3VEC", "LC_LOOP3"] <> " (", functionName, ",\n",
-        "  i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],\n",
-        "  cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]", If[vectorise, {",\n",
-        "  CCTK_REAL_VEC_SIZE"},""] <> ")\n",
-        "{\n",
-        indentBlock[
-          {
-             (* DeclareVariable["index", "// int"], *)
-             (* DeclareAssignVariable["int", "index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"], *)
-             DeclareAssignVariable["ptrdiff_t", "index", "di*i + dj*j + dk*k"],
-             block
-          }
-        ],
-        "}\n",
-        If[vectorise, "LC_ENDLOOP3VEC", "LC_ENDLOOP3"] <> " (", functionName, ");\n"
-      }
-    ],
-    ""
-  ];
-
-switchOptions[{value_, block_}] :=
-{
-  "case ", value, ":\n", indentBlock[{block,"break;\n"}]
-}
-
-SwitchStatement[var_, pairs__] :=
-{
-  "switch(", var, ")\n",
-  CBlock[{Riffle[Map[switchOptions, {pairs}],"\n"]}]
-}
-
-
-
-BoundaryLoop[block_] :=
-{
-  "\nGenericFD_GetBoundaryInfo(cctkGH, cctk_lsh, cctk_lssh, cctk_bbox, cctk_nghostzones, imin, imax, is_symbnd, is_physbnd, is_ipbnd);\n",
-
-  CommentedBlock["Start by looping over the whole grid, minus the NON-PHYSICAL boundary points, which are set by synchronization.  ", {
-  AssignVariable[arrayElement["bmin", 0], "is_physbnd[0*2+0] ? 0 : imin[0]"],
-  AssignVariable[arrayElement["bmin", 1], "is_physbnd[1*2+0] ? 0 : imin[1]"],
-  AssignVariable[arrayElement["bmin", 2], "is_physbnd[2*2+0] ? 0 : imin[2]"],
-  AssignVariable[arrayElement["bmax", 0], "is_physbnd[0*2+1] ? CCTK_LSSH(0,0) : imax[0]"],
-  AssignVariable[arrayElement["bmax", 1], "is_physbnd[1*2+1] ? CCTK_LSSH(0,1) : imax[1]"],
-  AssignVariable[arrayElement["bmax", 2], "is_physbnd[2*2+1] ? CCTK_LSSH(0,2) : imax[2]"]}], 
-
-  CommentedBlock["Loop over all faces",
-   loopOverInteger["dir", "0", "3",
-     loopOverInteger["face", "0", "2",
-     {
-      CommentedBlock["Now restrict to only the boundary points on the current face",
-       SwitchStatement["face", 
-        {0,  {AssignVariable[arrayElement["bmax", "dir"], {arrayElement["imin", "dir"], ""}], 
-              AssignVariable[arrayElement["bmin", "dir"], {0, ""}]}},
-        {1,  {AssignVariable[arrayElement["bmin", "dir"], {arrayElement["imax", "dir"], "" }],
-              AssignVariable[arrayElement["bmax", "dir"], {"CCTK_LSSH(0,dir)", ""}]}}]],
-       conditional[arrayElement["is_physbnd", "dir * 2 + face"],
-         loopOverInteger["k", arrayElement["bmin",2], arrayElement["bmax",2],
-           loopOverInteger["j", arrayElement["bmin",1], arrayElement["bmax",1],
-             loopOverInteger["i", arrayElement["bmin",0], arrayElement["bmax",0],
-
-         { If[SOURCELANGUAGE == "C",  
-              AssignVariable["index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"],
-              ""],
-	   block
-         }
-      
-      ]]]
-      ]}
-     ]]]};
-
-BoundaryWithGhostsLoop[block_] :=
-{
-  "\nGenericFD_GetBoundaryInfo(cctkGH, cctk_lsh, cctk_lssh, cctk_bbox, cctk_nghostzones, imin, imax, is_symbnd, is_physbnd, is_ipbnd);\n",
-
-  CommentedBlock["Start by looping over the whole grid, including the NON-PHYSICAL boundary points.  ", {
-  AssignVariable[arrayElement["bmin", 0], "0"],
-  AssignVariable[arrayElement["bmin", 1], "0"],
-  AssignVariable[arrayElement["bmin", 2], "0"],
-  AssignVariable[arrayElement["bmax", 0], "CCTK_LSSH(0,0)"],
-  AssignVariable[arrayElement["bmax", 1], "CCTK_LSSH(0,1)"],
-  AssignVariable[arrayElement["bmax", 2], "CCTK_LSSH(0,2)"]}], 
-
-  CommentedBlock["Loop over all faces",
-   loopOverInteger["dir", "0", "3",
-     loopOverInteger["face", "0", "2",
-     {
-      CommentedBlock["Now restrict to only the boundary points on the current face",
-       SwitchStatement["face", 
-        {0,  {AssignVariable[arrayElement["bmax", "dir"], {arrayElement["imin", "dir"], ""}], 
-              AssignVariable[arrayElement["bmin", "dir"], {0, ""}]}},
-        {1,  {AssignVariable[arrayElement["bmin", "dir"], {arrayElement["imax", "dir"], "" }],
-              AssignVariable[arrayElement["bmax", "dir"], {"CCTK_LSSH(0,dir)]", ""}]}}]],
-       conditional[arrayElement["is_physbnd", "dir * 2 + face"],
-         loopOverInteger["k", arrayElement["bmin",2], arrayElement["bmax",2],
-           loopOverInteger["j", arrayElement["bmin",1], arrayElement["bmax",1],
-             loopOverInteger["i", arrayElement["bmin",0], arrayElement["bmax",0],
-
-         { If[SOURCELANGUAGE == "C",  
-              AssignVariable["index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"],
-              ""],
-	   block
-         }
-      
-      ]]]
-      ]}
-     ]]]};
-
-Conditional[condition_, block_] :=
- {"if (", condition, ")\n",
-  CBlock[block]};
-
-Conditional[condition_, block1_, block2_] :=
- {"if (", condition, ")\n",
-  CBlock[block1], "else\n", CBlock[block2]};
-
-onceInGridLoop[block_] :=
-  Conditional["i == 5 && j == 5 && k == 5",
-              block];
-
-InfoVariable[name_] :=
- onceInGridLoop[
-  { "char buffer[255];\n",
-     "sprintf(buffer,\"" , name , " == %f\", " , name , ");\n",
-     "CCTK_INFO(buffer);\n"}];
-
-(* Code generation for Cactus .par files *)
-
-activeThorns[list_] :=
-  {"ActiveThorns = \"", SpaceSeparated[list], "\"\n"};
-
-setParameter[thorn_, par_, value_] :=
-  {thorn, " = ", If[NumberQ[value], ToString[value], "\"" <> value <> "\""], "\n"};
-
-setParametersForThorn[thorn_, map_] :=
-  Map[setParameter[thorn, #[[1]], #[[2]]] &, map];
-
-insertFile[name_] := 
-  Module[{istream_, contents_},
+DefFn[
+  AddToFile[filename_String, contents:CodeGenBlock] :=
+  Module[
+    {fp = OpenAppend[filename]},
+    WriteString[fp, FlattenBlock[contents]];
+    Close[fp]]];
+
+DefFn[
+  CommaNewlineSeparated[l_List] :=
+  Riffle[l, ",\n"]];
+
+DefFn[
+  SpaceSeparated[l_List] :=
+  Riffle[l, " "]];
+
+DefFn[
+  CommaSeparated[l_List] :=
+  Riffle[l, ", "]];
+
+DefFn[
+  NewlineSeparated[l_List] :=
+  Riffle[l, "\n"]];
+
+DefFn[
+  CommaInitSeparated[l_List] :=
+  Riffle[Map[{#," = INITVALUE"} &, l], ", "]];
+
+(* Turn a section of code into a string:
+   1. quote all quotes (replace all quotes with backslash-quote)
+   2. break the string into lines to make it readable (replace all newlines
+      with quote-newline-quote)
+   3. surround the result with quotes *)
+DefFn[
+  Stringify[x:CodeGenBlock] :=
+  "\"" <> StringReplace[StringReplace[FlattenBlock[x], "\"" -> "\\\""],
+                        "\n" -> "\\n\"\n\""] <> "\"\n"];
+
+DefFn[
+  PartitionVarList[list_List] :=
+  Module[
+    {partition, split},
+
+    partition[locallist_] :=
+    Module[
+      {cutoff},
+      cutoff = 6;
+      If[Length@locallist > cutoff, Partition[locallist, cutoff, cutoff, {1,1}, {}],
+         {locallist}]];
+
+    split = Split[list, NameRoot[#1] == NameRoot[#2] &];
+    split = Flatten[Map[partition, split], 1];
+
+    split]];
+
+DefFn[
+  insertFile[name_String] :=
+  Module[
+    {istream_, contents_},
     istream = OpenRead[name];
     contents = ReadList[istream, String];
     Close[istream];
-    contents];
-
-vectoriseExpression[exprp_] :=
-  Module[{isNotMinusOneQ, isNotTimesMinusOneQ, fmaRules, isNotKneg, arithRules, undoRules, expr},
-    expr = exprp;
-
-    (* Constants *)
-    (* expr = expr /. xx_Integer/; xx!=-1 :> ToReal[xx]; *)
-    expr = expr /. xx_Integer -> ToReal[xx];
-    expr = expr /. xx_Real -> ToReal[xx];
-    removeToRealRules = {
-      - ToReal[xx_] -> ToReal[- xx],
-      ToReal[xx_] + ToReal[yy_] -> ToReal[xx + yy],
-      ToReal[xx_] - ToReal[yy_] -> ToReal[xx - yy],
-      ToReal[xx_] * ToReal[yy_] -> ToReal[xx * yy],
-      ToReal[xx_] == ToReal[yy_] -> ToReal[xx == yy],
-      ToReal[xx_] != ToReal[yy_] -> ToReal[xx != yy],
-      pow[xx_, ToReal[power_]] -> pow[xx, power]
-    };
-    expr = expr //. removeToRealRules;
-
-    (* Replace all operators and functions *)
-    (* kneg, kadd, ksub, kmul, kdiv *)
-    isNotKneg[n_] := ! MatchQ[n,kneg[_]];
-    arithRules = {
-      - xx_ -> kneg[xx],
-      xx_ * yy_ -> kmul[xx,yy],
-      xx_ / yy_ -> kdiv[xx,yy],
-      xx_ + yy_ -> kadd[xx,yy],
-      xx_ - yy_ -> ksub[xx,yy],
-      kmul[-1,xx_]            -> kneg[xx],
-      kmul[xx_,-1]            -> kneg[xx],
-      kmul[ToReal[-1],xx_]    -> kneg[xx],
-      kmul[xx_,ToReal[-1]]    -> kneg[xx],
-      ToReal[- xx_]           -> kneg[ToReal[xx]],
-      (* kmul[xx_,INV[yy_]]      -> kdiv[xx,yy], *)
-      (* kmul[INV[xx_],yy_]      -> kdiv[yy,xx], *)
-      kdiv[xx_,kdiv[yy_,zz_]] -> kdiv[kmul[xx,zz],yy],
-      kdiv[kdiv[xx_,yy_],zz_] -> kdiv[xx,kmul[yy,zz]],
-      kmul[kneg[xx_],yy_]     -> kneg[kmul[xx,yy]],
-      kmul[xx_,kneg[yy_]]     -> kneg[kmul[xx,yy]],
-      kdiv[kneg[xx_],yy_]     -> kneg[kdiv[xx,yy]],
-      kdiv[xx_,kneg[yy_]]     -> kneg[kdiv[xx,yy]],
-      kadd[kneg[xx_],yy_]     -> ksub[yy,xx],
-      ksub[kneg[xx_],yy_]     -> kneg[kadd[xx,yy]],
-      kadd[xx_,kneg[yy_]]     -> ksub[xx,yy],
-      ksub[xx_,kneg[yy_]]     -> kadd[xx,yy],
-      kneg[ksub[xx_,yy_]]     -> ksub[yy,xx],
-      Abs[xx_]      -> kfabs[xx],
-      Log[xx_]      -> klog[xx],
-      fabs[xx_]     -> kfabs[xx],
-      fmax[xx_,yy_] -> kfmax[xx,yy],
-      fmin[xx_,yy_] -> kfmin[xx,yy],
-      sqrt[xx_]     -> ksqrt[xx],
-      exp[xx_]      -> kexp[xx],
-      log[xx_]      -> klog[xx],
-      pow[xx_,yy_]  -> kpow[xx,yy],
-      kneg[kneg[xx_]]   -> xx,
-      kfabs[kneg[xx_]]  -> kfabs[xx],
-      kfnabs[kneg[xx_]] -> kfnabs[xx],
-      kneg[kfabs[xx_]]  -> kfnabs[xx],
-      kneg[kfnabs[xx_]] -> kfabs[xx]
-    };
-    expr = expr //. arithRules;
-
-    (* Undo some transformations *)
-    undoRules = {
-      IfThen[_, aa_, aa_] -> aa,
-      IfThen[xx_?IntegerQ, aa_, bb_] /; xx!=0 :> aa,
-      IfThen[xx_?IntegerQ, aa_, bb_] /; xx==0 :> bb,
-      IfThen[kmul[xx_,yy_], aa_, bb_] -> IfThen[xx*yy, aa, bb],
-      IfThen[kmul[xx_,yy_] != zz_, aa_, bb_] -> IfThen[xx*yy!=zz, aa, bb],
-      IfThen[ToReal[xx_], aa_, bb_] -> IfThen[xx, aa, bb],
-      Scalar[kmul[xx_,yy_]] -> Scalar[xx*yy],
-      Scalar[kmul[xx_,yy_] != zz_] -> Scalar[xx*yy!=zz],
-      Scalar[ToReal[xx_]] -> Scalar[xx],
-      Scalar[xx_ != ToReal[yy_]] -> Scalar[xx != yy],
-      ToReal[kneg[xx_]] -> ToReal[-xx],
-      ToReal[kadd[xx_,yy_]] -> ToReal[xx+yy],
-      ToReal[ksub[xx_,yy_]] -> ToReal[xx-yy],
-      ToReal[kmul[xx_,yy_]] -> ToReal[xx*yy],
-      ToReal[xx_*kadd[yy_,zz_]] -> ToReal[xx*(yy+zz)],
-      kpow[xx_, kneg[power_]] -> kpow[xx, -power]
-    };
-    expr = expr //. undoRules;
-    
-    (* FMA (fused multiply-add) instructions *)
-    (* kmadd (x,y,z) =   xy+z
-       kmsub (x,y,z) =   xy-z
-       knmadd(x,y,z) = -(xy+z)
-       knmsub(x,y,z) = -(xy-z) *)
-    fmaRules = {
-      kadd[kmul[xx_,yy_],zz_] -> kmadd[xx,yy,zz],
-      kadd[zz_,kmul[xx_,yy_]] -> kmadd[xx,yy,zz],
-      ksub[kmul[xx_,yy_],zz_] -> kmsub[xx,yy,zz],
-      ksub[zz_,kmul[xx_,yy_]] -> knmsub[xx,yy,zz],
-      kneg[kmadd [xx_,yy_,zz_]] -> knmadd[xx,yy,zz],
-      kneg[kmsub [xx_,yy_,zz_]] -> knmsub[xx,yy,zz],
-      kneg[knmadd[xx_,yy_,zz_]] -> kmadd [xx,yy,zz],
-      kneg[knmsub[xx_,yy_,zz_]] -> kmsub [xx,yy,zz]
-      (* we could match this and similar patterns
-      kmul[xx_, kadd[yy_, ToReal[+1]]] -> kmadd[xx, yy, xx],
-      kmul[xx_, kadd[yy_, ToReal[-1]]] -> kmsub[xx, yy, xx],
-      *)
-    };
-    expr = expr //. fmaRules;
-    
-    Return[expr]];
-
-(* Take an expression x and replace occurrences of Powers with the C
-macros SQR, CUB, QAD *)
-ReplacePowers[expr_, vectorise_] :=
-  Module[{rhs},
-    rhs = expr /. Power[xx_, -1] -> INV[xx];
-
-    If[SOURCELANGUAGE == "C",
-           Module[{},
-	     rhs = rhs /. Power[xx_,  2  ] -> SQR[xx];
-	     rhs = rhs /. Power[xx_,  3  ] -> CUB[xx];
-	     rhs = rhs /. Power[xx_,  4  ] -> QAD[xx];
-	     rhs = rhs /. Power[xx_, -2  ] -> INV[SQR[xx]];
-             rhs = rhs /. Power[xx_,  1/2] -> sqrt[xx];
-             rhs = rhs /. Power[xx_, -1/2] -> INV[sqrt[xx]];
-             rhs = rhs /. Power[xx_,  0.5] -> sqrt[xx];
-             rhs = rhs /. Power[xx_, -0.5] -> INV[sqrt[xx]];
-
-             (*
-	     rhs = rhs /.  1/2 ->  khalf
-	     rhs = rhs /. -1/2 -> -khalf;
-
-	     rhs = rhs /.  1/3 ->  kthird;
-	     rhs = rhs /. -1/3 -> -kthird;
-
-	     rhs = rhs /.  2/3 ->  ktwothird;
-	     rhs = rhs /. -2/3 -> -ktwothird;
-
-	     rhs = rhs /.  4/3 ->  kfourthird;
-	     rhs = rhs /. -4/3 -> -kfourthird;
-
-	     rhs = rhs /.  8/3 ->  keightthird;
-	     rhs = rhs /. -8/3 -> -keightthird;
-             *)
-
-             (* Avoid rational numbers *)
-	     rhs = rhs /. Rational[xx_,yy_] :> N[xx/yy, 30];
-
-             rhs = rhs //.      IfThen[cond1_,xx1_,yy1_] +      IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[    xx1 +     xx2], Simplify[    yy1 +     yy2]];
-             rhs = rhs //. ff1_ IfThen[cond1_,xx1_,yy1_] +      IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[ff1 xx1 +     xx2], Simplify[ff1 yy1 +     yy2]];
-             rhs = rhs //.      IfThen[cond1_,xx1_,yy1_] + ff2_ IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[    xx1 + ff2 xx2], Simplify[    yy1 + ff2 yy2]];
-             rhs = rhs //. ff1_ IfThen[cond1_,xx1_,yy1_] + ff2_ IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[ff1 xx1 + ff2 xx2], Simplify[ff1 yy1 + ff2 yy2]];
-
-             (* Is this still a good idea when FMA instructions are used? *)
-	     rhs = rhs //. xx_ yy_ + xx_ zz_ -> xx (yy+zz);
-	     rhs = rhs //. xx_ yy_ - xx_ zz_ -> xx (yy-zz);
-
-             rhs = rhs /. Power[E, power_] -> exp[power];
-
-             (* there have been some problems doing the Max/Min
-                replacement via the preprocessor for C, so we do it
-                here *)
-             (* Note: Mathematica simplifies Max[xx_] -> xx automatically *)
-             rhs = rhs //. Max[xx_, yy__] -> fmax[xx, Max[yy]];
-             rhs = rhs //. Min[xx_, yy__] -> fmin[xx, Min[yy]];
-
-             rhs = rhs /. Power[xx_, power_] -> pow[xx, power];
-
-             If[vectorise === True,
-              rhs = vectoriseExpression[rhs]];
-
-             (* Remove Scalar[] after vectorising *)
-             rhs = rhs /. Scalar[xx_] -> xx;
-           ],
-
-           rhs = rhs /. Power[xx_, power_] -> xx^power
-       ];
-(*       Print[rhs//FullForm];*)
-    rhs
-    ];
-
-(* Convert an expression to CForm, but remove the quotes from any
-   strings present *)
-CFormHideStrings[x_, opts___] := StringReplace[ToString[CForm[x,opts]], "\"" -> ""];
-
-
-
-Quote[x_] := {"\"", x, "\""};
+    contents]];
+
+DefFn[
+  NameRoot[name_Symbol] :=
+  Module[
+    {dropNumberRule, root},
+    dropNumberRule = {"1" -> "", "2" -> "", "3" -> "", "4" -> "", "5" -> "",
+                      "6" -> "", "7" -> "", "8" -> "", "9" -> "", "0" -> "", "rhs" -> ""};
+    root = StringReplace[ToString@name, dropNumberRule]]];
+
+DefFn[
+  Quote[x:CodeGenBlock] :=
+  {"\"", x, "\""}];
 
 End[];
 
diff --git a/Tools/CodeGen/CodeGenC.m b/Tools/CodeGen/CodeGenC.m
new file mode 100644
index 0000000..95d1591
--- /dev/null
+++ b/Tools/CodeGen/CodeGenC.m
@@ -0,0 +1,245 @@
+
+(* $Id$ *)
+
+(*  Copyright 2004 Sascha Husa, Ian Hinder, Christiane Lechner
+
+    This file is part of Kranc.
+
+    Kranc is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    Kranc is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Kranc; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*)
+
+BeginPackage["CodeGenC`", {"Errors`", "Kranc`", "CodeGen`"}];
+
+SOURCELANGUAGE::usage = "global variable == \"C\" or \"Fortran\" determines language
+                        for code generation";
+SOURCESUFFIX::usage = "file suffix for source files";
+SetSourceLanguage::usage = "set the source language to  \"C\" or \"Fortran\"";
+
+EOL::usage = "the end of line termination string";
+
+CommentedBlock::usage = "CommentedBlock[comment, block] returns a block consisting " <>
+  "of 'comment' followed by 'block'.";
+CBlock::usage = "";
+SuffixedCBlock::usage = "";
+IncludeFile::usage = "IncludeFile[name] returns a block of code" <>
+  "that includes a header file (i.e '#include \"name\"').";
+IncludeSystemFile::usage = "IncludeFile[name] returns a block of code" <>
+  "that includes a system header file (i.e '#include <name>').";
+DeclareVariable::usage = "DeclareVariable[name, type] returns a block of code " <>
+  "that declares a variable of given name and type.  'name' and 'type' should be " <>
+  "strings.";
+DeclareVariableNoInit::usage = "DeclareVariableNoInit[name, type] returns a block of code " <>
+  "that declares a variable of given name and type without initialising it.  'name' and 'type' should be " <>
+  "strings.";
+DeclareVariables::usage = "DeclareVariables[names, type] returns a block of code " <>
+  "that declares a list of variables of given name and type.  'names' should be a list" <>
+  " of strings and 'type' should be a string string.";
+DeclarePointer::usage = "DeclarePointer[name, type] returns a block of code " <>
+  "that declares a pointer of given name and type.  'name' and 'type' should be " <>
+  "strings.";
+DeclarePointers::usage = "DeclarePointers[names, type] returns a block of code " <>
+  "that declares a list of pointers of given name and type.  'names' should be a list" <>
+  " of strings and 'type' should be a string string.";
+DefineVariable::usage = "DefineVariable[name, type, value] returns a block of " <>
+  "code that declares and initialised a variable 'name' of type 'type' to value 'value'.";
+AssignVariable::usage = "AssignVariable[dest_, src_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+DeclareAssignVariable::usage = "DeclareAssignVariable[type_, dest_, src_] returns a block of code " <>
+  "that declares and sets a constant variable of given name and type.";
+DeclareArray::usage = "";
+DefineFunction::usage = "";
+DefineSubroutine::usage = "";
+Conditional::usage = "";
+SwitchStatement::usage = "";
+CFormHideStrings::usage = "";
+InsertComment::usage = "";
+
+Begin["`Private`"];
+
+SOURCELANGUAGE =  "C";
+SOURCESUFFIX   = ".cc";
+
+setSourceSuffix[lang_] :=
+  If[lang == "C", SOURCESUFFIX = ".cc", SOURCESUFFIX = ".F90"];
+
+SetSourceLanguage[lang_] :=
+  If[lang == "C" || lang == "Fortran",
+     SOURCELANGUAGE = lang;
+     setSourceSuffix[lang];
+     InfoMessage[Terse, "User set source language to " <> lang],
+     (* else *)
+     SOURCELANGUAGE = "C";
+     setSourceSuffix[".cc"];
+     InfoMessage[Terse, "Setting Source Language to C"]];
+
+EOL[dummy___] :=
+  If[SOURCELANGUAGE == "C" || SOURCELANGUAGE == "C++", ";\n", "\n"];
+
+DefFn[
+  IncludeFile[filename_String] :=
+  {"#include \"", filename, "\"\n"}];
+
+DefFn[
+  IncludeSystemFile[filename_String] :=
+  {"#include <", filename, ">\n"}];
+
+DefFn[
+  DeclareVariable[name:(_String|_Symbol), type_String] :=
+  If[SOURCELANGUAGE == "C",
+     {type, " ",    name, " = INITVALUE" <> EOL[]},
+     {type, " :: ", name, EOL[]} (* no value init here to avoid implicit SAVE attribute *)]];
+
+DefFn[
+  DeclareVariableNoInit[name:(_String|_Symbol), type_String] :=
+  If[SOURCELANGUAGE == "C",
+     {type, " ",    name, EOL[]},
+     {type, " :: ", name, EOL[]} (* no value init here to avoid implicit SAVE attribute *)]];
+
+DefFn[
+  DeclareVariables[names_?ListQ, type_String] := 
+  If[SOURCELANGUAGE == "C",
+     {type, " ",    CommaSeparated@names, EOL[]},
+     {type, " :: ", CommaSeparated@names,     EOL[]} (* no value init avoids implicit SAVE attribute *)]];
+
+DefFn[
+  DeclarePointer[name:(_String|_Symbol), type_String] :=
+  If[SOURCELANGUAGE == "C",
+     {type, " *",    name, EOL[]},
+     {type, ", target :: ", name, EOL[]}]];
+
+DefFn[
+  DeclarePointers[names_?ListQ, type_String] :=
+  If[SOURCELANGUAGE == "C",
+     {type, " *",           CommaInitSeparated@names, EOL[]},
+     {type, ", target :: ", CommaSeparated@names,     EOL[]}]];
+
+DefFn[
+  DeclareArray[name:(_String|_Symbol), dim_Integer, type_String] :=
+  If[SOURCELANGUAGE == "C",
+     DeclareArrayC[name, dim, type],
+     DeclareArrayFortran[name, dim, type]]];
+
+DefFn[
+  DeclareArrayC[name:(_String|_Symbol), dim_Integer, type_String] :=
+  {type, " ", name, "[", dim, "];","\n"}];
+
+DefFn[
+  DeclareArrayFortran[name:(_String|_Symbol), dim_Integer, type_String] :=
+  {type, " :: ", name, "(", dim, ")","\n"}];
+
+DefFn[
+  DefineVariable[name:(_String|_Symbol), type_String, value:CodeGenBlock] :=
+  {type, " ", name, " = ", value, EOL[]}];
+
+DefFn[
+  AssignVariable[dest:(_String|_Symbol), src:CodeGenBlock] :=
+  {dest, " = ", src, EOL[]}];
+
+DefFn[
+  DeclareAssignVariable[type_String, dest:(_String|_Symbol), src:CodeGenBlock] :=
+  {type, " const ", dest, " = ", src, EOL[]}];
+
+(* comments are always done C-style because they are killed by cpp anyway *) 
+DefFn[
+  InsertComment[text:CodeGenBlock] := {"/* ", text, " */\n"}];
+
+DefFn[
+  CBlock[block:CodeGenBlock] :=
+  {"{\n",
+   IndentBlock[block],
+   "}\n"}];
+
+DefFn[
+  SuffixedCBlock[block:CodeGenBlock, suffix_] :=
+  {"{\n",
+   IndentBlock[block],
+   "} ", suffix, "\n"}];
+
+DefFn[
+  CommentedBlock[comment:CodeGenBlock, block:CodeGenBlock] :=
+  SeparatedBlock[{InsertComment[comment],
+                  block}]];
+
+(* FUNCTIONS *)
+
+DefFn[
+  defineFunctionC[name_String, type_String, args:CodeGenBlock, contents:CodeGenBlock] :=
+  SeparatedBlock[
+    {type, " ", name, "(", args, ")\n",
+     CBlock[contents]}]];
+     
+DefFn[
+  defineFunctionF[name_String, args:CodeGenBlock, contents:CodeGenBlock] :=
+  SeparatedBlock[
+    {"FUNCTION", " ", name, "(", args, ")\n",
+     IndentBlock[contents]}]];
+
+DefFn[
+  DefineFunction[name_String, type_String, args:CodeGenBlock, contents:CodeGenBlock] :=
+  If[SOURCELANGUAGE == "C",
+     defineFunctionC[name, type, args, contents],
+     defineFunctionF[name, args, contents]]];
+
+(* SUBROUTINES *)
+
+DefFn[
+  DefineSubroutine[name_String, args:CodeGenBlock, contents:CodeGenBlock] :=
+  If[SOURCELANGUAGE == "C",
+     DefineSubroutineC[name, args, contents],
+     DefineSubroutineF[name, args, contents]]];
+
+DefFn[
+  DefineSubroutineC[name_String, args:CodeGenBlock, contents:CodeGenBlock] :=
+  SeparatedBlock[
+    {"extern \"C\" void ", name, "(", args, ")", "\n",
+     CBlock[contents]}]];
+
+DefFn[
+  DefineSubroutineF[name_String, args:CodeGenBlock, contents:CodeGenBlock] :=
+  SeparatedBlock[
+    {"subroutine ", name, "(", args, ")", "\n",
+     "\nimplicit none\n\n",
+     contents,
+     "end subroutine\n"}]];
+
+DefFn[
+  switchOption[{value:(_String|_Symbol|_?NumberQ), block:CodeGenBlock}] :=
+  {"case ", value, ":\n", IndentBlock[{block,"break;\n"}]}];
+(* Outer list unnecessary? *)
+
+DefFn[
+  SwitchStatement[var:(_String|_Symbol), pairs__] :=
+  {"switch(", var, ")\n",
+   CBlock[{Riffle[Map[switchOption, {pairs}],"\n"]}]}];
+
+DefFn[
+  Conditional[condition:CodeGenBlock, block:CodeGenBlock] :=
+  {"if (", condition, ")\n",
+   CBlock[block]}];
+
+DefFn[
+  Conditional[condition:CodeGenBlock, block1:CodeGenBlock, block2:CodeGenBlock] :=
+  {"if (", condition, ")\n",
+   CBlock[block1], "else\n", CBlock[block2]}];
+
+(* Convert an expression to CForm, but remove the quotes from any
+   strings present *)
+DefFn[
+  CFormHideStrings[x_, opts___] :=
+  StringReplace[ToString[CForm[x,opts]], "\"" -> ""]];
+
+End[];
+
+EndPackage[];
diff --git a/Tools/CodeGen/CodeGenCactus.m b/Tools/CodeGen/CodeGenCactus.m
new file mode 100644
index 0000000..092b00f
--- /dev/null
+++ b/Tools/CodeGen/CodeGenCactus.m
@@ -0,0 +1,733 @@
+
+(* $Id$ *)
+
+(*  Copyright 2004 Sascha Husa, Ian Hinder, Christiane Lechner
+
+    This file is part of Kranc.
+
+    Kranc is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    Kranc is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Kranc; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*)
+
+BeginPackage["CodeGenCactus`", {"Errors`", "Kranc`", "CodeGenC`", "CodeGen`"}];
+
+AssignVariableInLoop::usage = "AssignVariableInLoop[dest_, src_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+StoreVariableInLoop::usage = "StoreVariableInLoop[dest_, src_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+StoreLowPartialVariableInLoop::usage = "StoreLowPartialVariableInLoop[dest_, src_, count_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+StoreHighPartialVariableInLoop::usage = "StoreHighPartialVariableInLoop[dest_, src_, count_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+StoreMiddlePartialVariableInLoop::usage = "StoreMiddlePartialVariableInLoop[dest_, src_, countLow_, countHigh_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+StorePartialVariableInLoop::usage = "StorePartialVariableInLoop[dest_, src_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+DeclareAssignVariableInLoop::usage = "DeclareAssignVariableInLoop[type_, dest_, src_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+MaybeAssignVariableInLoop::usage = "MaybeAssignVariableInLoop[dest_, src_, cond_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+DeclareMaybeAssignVariableInLoop::usage = "DeclareMaybeAssignVariableInLoop[type_, dest_, src_, cond_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+TestForNaN::usage = "TestForNaN[expr_] returns a block of code " <>
+  "that tests 'expr' for nan.";
+DefineCCTKFunction::usage = "DefineCCTKFunction[name, type, block] returns a block " <>
+  "of code that defines a CCTK function of name 'name' returning type 'type' with " <>
+  "body 'block'.";
+DefineCCTKSubroutine::usage = "DefineCCTKSubroutine[name, block] returns a block " <>
+  "of code that defines a CCTK Fortran subroutine of name 'name' with body 'block'.";
+GridName::usage = "GridName[variable] returns the name needed to access variable " <>
+  "assuming it is a grid variable when inside a grid loop.";
+DeclareGridLoopVariables::usage = "DeclareGridLoopVariables[] returns a block " <>
+  "that defines the variables needed during a grid loop.";
+InitialiseGridLoopVariables::usage = "InitialiseGridLoopVariables[] returns a block " <>
+  "that initialises variables needed by a grid loop.";
+InitialiseGridLoopVariablesWithStencil::usage = "InitialiseGridLoopVariables[] returns a block " <>
+  "that initialises variables needed by a grid loop using the evolution stencils.";
+ConditionalOnParameter::usage = "ConditionalOnParameter[name, value, block] returns " <>
+  "a block that introduces a conditional expression whereby 'block' is only executed " <>
+  "if the Cactus parameter 'name' has value 'value'.";
+(*
+GridLoop::usage = "GridLoop[block] returns a block that is looped over for every " <>
+  "grid point.  Must have previously set up the grid loop variables (see " <>
+  "DeclareGridLoopVariables and InitialiseGridLoopVariables.";
+*)
+GridLoop::usage = "GridLoop[block] returns a block that is looped over for every " <>
+  "grid point.  Must have previously set up the grid loop variables (see " <>
+  "InitialiseGridLoopVariables.";
+ConditionalOnParameterTextual::usage = "";
+DeclareFDVariables::usage = "";
+InitialiseFDVariables::usage = "";
+ReplacePowers::usage = "";
+BoundaryLoop::usage = "";
+BoundaryWithGhostsLoop::usage = "";
+GenericGridLoop::usage = "";
+NameRoot::usage = "";
+PartitionVarList::usage = "";
+DataType::usage = "DataType[] returns a string for the grid function data type (e.g. CCTK_REAL)";
+SetDataType::usage = "SetDataType[type] sets a string for the grid function data type (e.g. CCTK_REAL)";
+
+Begin["`Private`"];
+
+DefFn[
+  SetDataType[type_String] :=
+  dataType = type];
+
+DefFn[
+  DataType[] :=
+  If[dataType === Symbol["datatype"],
+     Throw["DataType: Have not set a data type"],
+     dataType]];
+
+DefFn[
+  AssignVariableInLoop[dest:(_String|_Symbol), src:CodeGenBlock, vectorise_:False] :=
+  Module[
+    {loader},
+    loader[x_] := If[vectorise, {"vec_load(", x, ")"}, x];
+    {dest, " = ", loader[src], EOL[]}]];
+
+DefFn[
+  StoreVariableInLoop[dest:(_String|_Symbol), src:(_String|_Symbol)] :=
+  {"vec_store_nta(", dest, ",", src, ")", EOL[]}];
+
+DefFn[
+  StoreLowPartialVariableInLoop[dest:(_String|_Symbol), src:(_String|_Symbol), count_String] :=
+  {"vec_store_nta_partial_lo(", dest, ",", src, ",", count, ")", EOL[]}];
+
+DefFn[
+  StoreHighPartialVariableInLoop[dest:(_String|_Symbol), src:(_String|_Symbol), count_String] :=
+  {"vec_store_nta_partial_hi(", dest, ",", src, ",", count, ")", EOL[]}];
+
+DefFn[
+  StoreMiddlePartialVariableInLoop[dest:(_String|_Symbol), src:(_String|_Symbol), countLow_String, countHigh_String] :=
+  {"vec_store_nta_partial_mid(", dest, ",", src, ",", countLow, ",", countHigh, ")", EOL[]}];
+
+DefFn[
+  StorePartialVariableInLoop[dest:(_String|_Symbol), src:(_String|_Symbol)] :=
+  {"vec_store_nta_partial(", dest, ",", src, ")", EOL[]}];
+
+DefFn[
+  DeclareAssignVariableInLoop[type_String, dest:(_String|_Symbol), src:(_String|_Symbol)] :=
+  {type, " const ", dest, " = vec_load(", src, ")", EOL[]}];
+
+DefFn[
+  MaybeAssignVariableInLoop[dest:(_String|_Symbol), src:(_String|_Symbol), cond:Boolean] :=
+  If[cond,
+     {dest, " = useMatter ? vec_load(", src, ") : ToReal(0.0)", EOL[]},
+     {dest, " = vec_load(", src, ")", EOL[]}]];
+
+DefFn[
+  DeclareMaybeAssignVariableInLoop[type_String, dest:(_String|_Symbol), src:(_String|_Symbol),
+                                 mmaCond:Boolean, codeCond:CodeGenBlock,
+                                 vectorise:Boolean:False] :=
+  Module[
+    {loader},
+    loader[x_] := If[vectorise, {"vec_load(", x, ")"}, x];
+    If[mmaCond,
+       {type, " ", dest, " = (", codeCond, ") ? ", loader[src], " : ToReal(0.0)", EOL[]},
+       {type, " ", dest, " = ", loader[src], EOL[]}]]];
+
+DefFn[
+  TestForNaN[expr:CodeGenBlock] :=
+  {"if (isnan(", expr, ")) {\n",
+   "  CCTK_VInfo(CCTK_THORNSTRING, \"NaN found\");\n",
+   "  CCTK_VInfo(CCTK_THORNSTRING, \"ipos: %d %d %d\", i, j, k);\n",
+   "  CCTK_VInfo(CCTK_THORNSTRING, \"lbnd: %d %d %d\", cctk_lbnd[0], cctk_lbnd[1], cctk_lbnd[2]);\n",
+   "  CCTK_VInfo(CCTK_THORNSTRING, \"lsh: %d %d %d\", cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]);\n",
+   "  CCTK_VInfo(CCTK_THORNSTRING, \"LSSH: %d %d %d\", CCTK_LSSH(0,0), CCTK_LSSH(0,1), CCTK_LSSH(0,2));\n",
+   "  CCTK_VInfo(CCTK_THORNSTRING, \"", expr, ": %.17g\", (double)", expr, ");\n",
+   "}\n"}];
+
+DefFn[
+  loopOverInteger[name_String, start_String, endplusone_String, block:CodeGenBlock] :=
+  If[SOURCELANGUAGE == "C" || SOURCELANGUAGE == "C++",
+     {"for (", name, " = ", start, "; ", name, " < ", endplusone, "; ", name, "++)\n",
+      "{\n",
+      IndentBlock[block],
+      "}\n"},
+     
+     {"Do ", name, " = ", start, ", ", endplusone, "\n",
+      "\n",
+      IndentBlock[block],
+      "End Do\n"}
+    ]];
+
+(* This is a Cactus-callable function *)
+DefFn[
+  DefineCCTKFunction[name_String, type_String, contents:CodeGenBlock] :=
+  DefineFunction[
+    name, "extern \"C\" " <> type, "CCTK_ARGUMENTS", 
+    {
+      "DECLARE_CCTK_ARGUMENTS;\n",
+      "DECLARE_CCTK_PARAMETERS;\n\n",
+      contents
+    }]];
+
+(* This is a Cactus-callable subroutine *)
+DefFn[
+  DefineCCTKSubroutine[name_String, contents:CodeGenBlock] :=
+  DefineSubroutine[
+    name, "CCTK_ARGUMENTS", 
+    {
+      "DECLARE_CCTK_ARGUMENTS;\n",
+      "DECLARE_CCTK_PARAMETERS;\n\n",
+      contents
+    }]];
+
+DefFn[
+  DeclareFDVariables[] := 
+(*
+  CommentedBlock["Declare finite differencing variables",
+    {Map[DeclareVariables[#, "CCTK_REAL"] &, {{"dx", "dy", "dz"}, 
+                                              {"dxi", "dyi", "dzi"},
+                                              {khalf,kthird,ktwothird,kfourthird,keightthird},
+                                              {"hdxi", "hdyi", "hdzi"}}],
+     "\n"},
+    {Map[DeclareVariables[#, "ptrdiff_t"] &, {{"di", "dj", "dk"}}],
+     "\n"}];
+*)
+  CommentedBlock["Declare finite differencing variables", {}]];
+
+DefFn[
+  InitialiseFDSpacingVariablesC[] :=
+  {
+    (* DeclareAssignVariable["ptrdiff_t", "di", "CCTK_GFINDEX3D(cctkGH,1,0,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"], *)
+    DeclareAssignVariable["ptrdiff_t", "di", "1"],
+    DeclareAssignVariable["ptrdiff_t", "dj", "CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"],
+    DeclareAssignVariable["ptrdiff_t", "dk", "CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0)"],
+    DeclareAssignVariable["ptrdiff_t", "cdi", "sizeof(CCTK_REAL) * di"],
+    DeclareAssignVariable["ptrdiff_t", "cdj", "sizeof(CCTK_REAL) * dj"],
+    DeclareAssignVariable["ptrdiff_t", "cdk", "sizeof(CCTK_REAL) * dk"],
+    DeclareAssignVariable[DataType[], "dx", "ToReal(CCTK_DELTA_SPACE(0))"],
+    DeclareAssignVariable[DataType[], "dy", "ToReal(CCTK_DELTA_SPACE(1))"],
+    DeclareAssignVariable[DataType[], "dz", "ToReal(CCTK_DELTA_SPACE(2))"],
+    DeclareAssignVariable[DataType[], "dt", "ToReal(CCTK_DELTA_TIME)"],
+    DeclareAssignVariable[DataType[], "t", "ToReal(cctk_time)"]}];
+
+DefFn[
+  InitialiseFDSpacingVariablesFortran[] := 
+  {
+    AssignVariable["dt", "CCTK_DELTA_TIME"],
+    AssignVariable["dx", "CCTK_DELTA_SPACE(1)"],
+    AssignVariable["dy", "CCTK_DELTA_SPACE(2)"],
+    AssignVariable["dz", "CCTK_DELTA_SPACE(3)"]
+  }];
+
+DefFn[
+  InitialiseFDVariables[vectorise:Boolean] :=
+  CommentedBlock[
+    "Initialise finite differencing variables",
+    {
+      If[SOURCELANGUAGE == "Fortran",
+         InitialiseFDSpacingVariablesFortran[],
+         InitialiseFDSpacingVariablesC[]],
+      
+      DeclareAssignVariable[DataType[], "dxi", "INV(dx)"],
+      DeclareAssignVariable[DataType[], "dyi", "INV(dy)"],
+      DeclareAssignVariable[DataType[], "dzi", "INV(dz)"],
+      If[vectorise,
+         {DeclareAssignVariable[DataType[], "khalf", "ToReal(0.5)"],
+          DeclareAssignVariable[DataType[], "kthird", "ToReal(1.0/3.0)"],
+          DeclareAssignVariable[DataType[], "ktwothird", "ToReal(2.0/3.0)"],
+          DeclareAssignVariable[DataType[], "kfourthird", "ToReal(4.0/3.0)"],
+          DeclareAssignVariable[DataType[], "keightthird", "ToReal(8.0/3.0)"],
+          DeclareAssignVariable[DataType[], "hdxi", "kmul(ToReal(0.5), dxi)"],
+          DeclareAssignVariable[DataType[], "hdyi", "kmul(ToReal(0.5), dyi)"],
+          DeclareAssignVariable[DataType[], "hdzi", "kmul(ToReal(0.5), dzi)"]},
+         (* else *)
+         {DeclareAssignVariable[DataType[], "khalf", "0.5"],
+          DeclareAssignVariable[DataType[], "kthird", "1/3.0"],
+          DeclareAssignVariable[DataType[], "ktwothird", "2.0/3.0"],
+          DeclareAssignVariable[DataType[], "kfourthird", "4.0/3.0"],
+          DeclareAssignVariable[DataType[], "keightthird", "8.0/3.0"],
+          DeclareAssignVariable[DataType[], "hdxi", "0.5 * dxi"],
+          DeclareAssignVariable[DataType[], "hdyi", "0.5 * dyi"],
+          DeclareAssignVariable[DataType[], "hdzi", "0.5 * dzi"]}]}]];
+
+DefFn[
+  GridName[x_Symbol] :=
+  If[SOURCELANGUAGE == "C",
+     ToString[x] <> "[index]",
+     ToString[x] <> "(i,j,k)"]];
+
+DefFn[
+  DeclareGridLoopVariables[] :=
+  SeparatedBlock[
+    {InsertComment["Declare the variables used for looping over grid points"],
+     Map[DeclareVariables[#, "CCTK_INT"] &, 
+         {{"i", "j", "k"}
+          (*, {"istart", "jstart", "kstart"}, 
+            {"iend", "jend", "kend"},
+            {"index_offset_x", "index_offset_y", "index_offset_z", "dir", "face"} *)}]
+     (*, Map[DeclareArray[#, 6, "CCTK_INT"] &, {"is_symbnd", "is_physbnd", "is_ipbnd"}],
+       Map[DeclareArray[#, 3, "CCTK_INT"] &, {"imin", "imax", "bmin", "bmax"}] *), 
+     If[SOURCELANGUAGE == "C", DeclareVariable["index", "// CCTK_INT"], "\n"]}]];
+
+(* Access an element of an array; syntax is different between C and
+   Fortran.  Always give this function a C-style array index. *)
+DefFn[
+  arrayElement[var_String, i_String] :=
+  If[SOURCELANGUAGE == "C",
+     {var, "[", arrayIndex[i], "]"},
+     {var, "(", arrayIndex[i], ")"}]];
+
+(* Given a C-style variable index, return the corresponding index for
+   the language currently in use.  The idea is that the caller does not
+   need to know what language is being used. *)
+DefFn[
+  arrayIndex[i:(_Integer|_String|_Symbol)] :=
+  If[SOURCELANGUAGE == "C",
+     i,
+     If[NumberQ[i], i+1, {i, " + 1"}]]];
+
+DefFn[
+  max[]:= If[SOURCELANGUAGE == "C", "IMAX", "max"]];
+
+DefFn[
+  InitialiseGridLoopVariables[derivativesUsedSwitch:Boolean, addToStencilWidth_Integer] :=
+  CommentedBlock[
+    "Set up variables used in the grid loop for the physical grid points",
+
+    If[derivativesUsedSwitch,
+       {AssignVariable["index_offset_x", max[] <>"(stencil_width, stencil_width_x) + " <> ToString[addToStencilWidth]],
+        AssignVariable["index_offset_y", max[] <>"(stencil_width, stencil_width_y) + " <> ToString[addToStencilWidth]],
+        AssignVariable["index_offset_z", max[] <>"(stencil_width, stencil_width_z) + " <> ToString[addToStencilWidth]],
+        "\n",
+        AssignVariable["istart", arrayIndex["index_offset_x"]],
+        AssignVariable["jstart", arrayIndex["index_offset_y"]],
+        AssignVariable["kstart", arrayIndex["index_offset_z"]],
+        "\n",
+        AssignVariable["iend", {"CCTK_LSSH(0,0) - index_offset_x"}],
+        AssignVariable["jend", {"CCTK_LSSH(0,1) - index_offset_y"}],
+        AssignVariable["kend", {"CCTK_LSSH(0,2) - index_offset_z"}]},
+       (* else *)
+       {AssignVariable["istart", arrayIndex[0]],
+        AssignVariable["jstart", arrayIndex[0]],
+        AssignVariable["kstart", arrayIndex[0]],
+        "\n",
+        AssignVariable["iend", "CCTK_LSSH(0,0)"],
+        AssignVariable["jend", "CCTK_LSSH(0,1)"],
+        AssignVariable["kend", "CCTK_LSSH(0,2)"]}]]];
+
+DefFn[
+  ConditionalOnParameter[name_String, value_String, block:CodeGenBlock] :=
+  SeparatedBlock[
+    {"if (CCTK_EQUALS(", name, ", \"", value, "\"))\n",
+     "{\n",
+     IndentBlock[block],
+     "}\n"}]];
+
+DefFn[
+  ConditionalOnParameterTextual[text:CodeGenBlock, block:CodeGenBlock] :=
+  SeparatedBlock[
+    {"if (", text, ")\n",
+     "{\n",
+     IndentBlock[block],
+     "}\n"}]];
+
+(*
+GridLoop[block_] :=
+  CommentedBlock["Loop over the grid points",
+   loopOverInteger["k", "kstart", "kend",
+     loopOverInteger["j", "jstart", "jend",
+       loopOverInteger["i", "istart", "iend",
+
+       { If[SOURCELANGUAGE == "C",  
+            AssignVariable["index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"],
+            ""],
+	 block
+       }
+        ]]]];
+
+*)
+
+(*
+GridLoop[block_] :=
+  If[SOURCELANGUAGE == "C",  
+    CommentedBlock["Loop over the grid points",
+      {
+        "#pragma omp parallel\n",
+        "LC_LOOP3 (unnamed,\n",
+        "          i,j,k, istart,jstart,kstart, iend,jend,kend,\n",
+        "          cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])\n",
+        "{\n",
+        IndentBlock[
+          {
+             DeclareVariable["index", "int"],
+             AssignVariable["index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"],
+             block
+          }
+        ],
+        "}\n",
+        "LC_ENDLOOP3 (unnamed);\n"
+      }
+    ],
+    CommentedBlock["Loop over the grid points",
+      {
+        "#pragma omp parallel\n",
+        "LC_LOOP3 (unnamed,\n",
+        "          i,j,k, istart,jstart,kstart, iend,jend,kend,\n",
+        "          cctk_lsh(1),cctk_lsh(2),cctk_lsh(3))\n",
+        IndentBlock[block],
+        "LC_ENDLOOP3 (unnamed)\n"
+      }
+    ]
+  ];
+*)
+
+Options[GenericGridLoop] = ThornOptions;
+
+DefFn[
+  GenericGridLoop[functionName_String, block:CodeGenBlock, opts:OptionsPattern[]] :=
+  If[OptionValue[UseLoopControl],
+     GenericGridLoopUsingLoopControl[functionName, block, OptionValue[UseVectors]],
+     GenericGridLoopTraditional[block]]];
+
+DefFn[
+  GenericGridLoopTraditional[block:CodeGenBlock] :=
+  CommentedBlock[
+    "Loop over the grid points",
+    loopOverInteger[
+      "k", "imin[2]", "imax[2]",
+      loopOverInteger[
+        "j", "imin[1]", "imax[1]",
+        loopOverInteger[
+          "i", "imin[0]", "imax[0]",
+          {If[SOURCELANGUAGE == "C",
+              DeclareAssignVariable["int", "index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"],
+              ""],
+            block}]]]]];
+
+DefFn[
+  GenericGridLoopUsingLoopControl[functionName_String, block:CodeGenBlock, vectorise:Boolean] :=
+  If[SOURCELANGUAGE == "C",  
+     CommentedBlock[
+       "Loop over the grid points",
+       {"#pragma omp parallel\n",
+        If[vectorise, "LC_LOOP3VEC", "LC_LOOP3"],
+        " (", functionName, ",\n",
+        "  i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],\n",
+        "  cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]",
+        If[vectorise, {",\n", "  CCTK_REAL_VEC_SIZE"}, ""],
+        ")\n",
+        "{\n",
+        IndentBlock[
+          {(* DeclareVariable["index", "// int"], *)
+           (* DeclareAssignVariable["int", "index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"], *)
+           DeclareAssignVariable["ptrdiff_t", "index", "di*i + dj*j + dk*k"],
+           block}], "}\n",
+        If[vectorise, "LC_ENDLOOP3VEC", "LC_ENDLOOP3"] <> " (", functionName, ");\n"}],
+     (* else *)
+     ""]];
+
+DefFn[
+  BoundaryLoop[block:CodeGenBlock] :=
+  {
+    "\nGenericFD_GetBoundaryInfo(cctkGH, cctk_lsh, cctk_lssh, cctk_bbox, cctk_nghostzones, imin, imax, is_symbnd, is_physbnd, is_ipbnd);\n",
+
+    CommentedBlock[
+      "Start by looping over the whole grid, minus the NON-PHYSICAL boundary points, which are set by synchronization.  ", {
+        AssignVariable[arrayElement["bmin", 0], "is_physbnd[0*2+0] ? 0 : imin[0]"],
+        AssignVariable[arrayElement["bmin", 1], "is_physbnd[1*2+0] ? 0 : imin[1]"],
+        AssignVariable[arrayElement["bmin", 2], "is_physbnd[2*2+0] ? 0 : imin[2]"],
+        AssignVariable[arrayElement["bmax", 0], "is_physbnd[0*2+1] ? CCTK_LSSH(0,0) : imax[0]"],
+        AssignVariable[arrayElement["bmax", 1], "is_physbnd[1*2+1] ? CCTK_LSSH(0,1) : imax[1]"],
+        AssignVariable[arrayElement["bmax", 2], "is_physbnd[2*2+1] ? CCTK_LSSH(0,2) : imax[2]"]}], 
+
+    CommentedBlock[
+      "Loop over all faces",
+      loopOverInteger[
+        "dir", "0", "3",
+        loopOverInteger[
+          "face", "0", "2",
+          {CommentedBlock[
+            "Now restrict to only the boundary points on the current face",
+            SwitchStatement[
+              "face", 
+              {0,  {AssignVariable[arrayElement["bmax", "dir"], {arrayElement["imin", "dir"], ""}], 
+                    AssignVariable[arrayElement["bmin", "dir"], {0, ""}]}},
+              {1,  {AssignVariable[arrayElement["bmin", "dir"], {arrayElement["imax", "dir"], "" }],
+                    AssignVariable[arrayElement["bmax", "dir"], {"CCTK_LSSH(0,dir)", ""}]}}]],
+           conditional[
+             arrayElement["is_physbnd", "dir * 2 + face"],
+             loopOverInteger[
+               "k", arrayElement["bmin",2], arrayElement["bmax",2],
+               loopOverInteger[
+                 "j", arrayElement["bmin",1], arrayElement["bmax",1],
+                 loopOverInteger[
+                   "i", arrayElement["bmin",0], arrayElement["bmax",0],
+                   {If[SOURCELANGUAGE == "C", AssignVariable["index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"], ""],
+                    block}]]]]}]]]}];
+
+DefFn[
+  BoundaryWithGhostsLoop[block:CodeGenBlock] :=
+  {
+  "\nGenericFD_GetBoundaryInfo(cctkGH, cctk_lsh, cctk_lssh, cctk_bbox, cctk_nghostzones, imin, imax, is_symbnd, is_physbnd, is_ipbnd);\n",
+
+  CommentedBlock[
+    "Start by looping over the whole grid, including the NON-PHYSICAL boundary points.  ", {
+  AssignVariable[arrayElement["bmin", 0], "0"],
+  AssignVariable[arrayElement["bmin", 1], "0"],
+  AssignVariable[arrayElement["bmin", 2], "0"],
+  AssignVariable[arrayElement["bmax", 0], "CCTK_LSSH(0,0)"],
+  AssignVariable[arrayElement["bmax", 1], "CCTK_LSSH(0,1)"],
+  AssignVariable[arrayElement["bmax", 2], "CCTK_LSSH(0,2)"]}], 
+
+  CommentedBlock[
+    "Loop over all faces",
+   loopOverInteger[
+     "dir", "0", "3",
+     loopOverInteger[
+       "face", "0", "2",
+     {
+      CommentedBlock[
+        "Now restrict to only the boundary points on the current face",
+       SwitchStatement[
+         "face", 
+        {0,  {AssignVariable[arrayElement["bmax", "dir"], {arrayElement["imin", "dir"], ""}], 
+              AssignVariable[arrayElement["bmin", "dir"], {0, ""}]}},
+        {1,  {AssignVariable[arrayElement["bmin", "dir"], {arrayElement["imax", "dir"], "" }],
+              AssignVariable[arrayElement["bmax", "dir"], {"CCTK_LSSH(0,dir)]", ""}]}}]],
+       conditional[arrayElement["is_physbnd", "dir * 2 + face"],
+         loopOverInteger[
+           "k", arrayElement["bmin",2], arrayElement["bmax",2],
+           loopOverInteger[
+             "j", arrayElement["bmin",1], arrayElement["bmax",1],
+             loopOverInteger[
+               "i", arrayElement["bmin",0], arrayElement["bmax",0],
+
+         {If[SOURCELANGUAGE == "C", AssignVariable["index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"], ""],
+	   block}]]]]}]]]}];
+
+DefFn[
+  onceInGridLoop[block:CodeGenBlock] :=
+  Conditional[
+    "i == 5 && j == 5 && k == 5",
+    block]];
+
+DefFn[
+  InfoVariable[name_String] :=
+  onceInGridLoop[
+    {"char buffer[255];\n",
+     "sprintf(buffer,\"" , name , " == %f\", " , name , ");\n",
+     "CCTK_INFO(buffer);\n"}]];
+
+(* Code generation for Cactus .par files *)
+
+DefFn[
+  activeThorns[list:{_String...}] :=
+  {"ActiveThorns = \"", SpaceSeparated[list], "\"\n"}];
+
+DefFn[
+  setParameter[thorn_String, par_String, value_] :=
+  {thorn, " = ", If[NumberQ[value], ToString[value], "\"" <> value <> "\""], "\n"}];
+
+DefFn[
+  setParametersForThorn[thorn_String, map_List] :=
+  Map[setParameter[thorn, #[[1]], #[[2]]] &, map]];
+
+DefFn[
+  vectoriseExpression[exprp_] :=
+  Module[
+    {isNotMinusOneQ, isNotTimesMinusOneQ, fmaRules, isNotKneg, arithRules, undoRules, expr},
+    expr = exprp;
+
+    (* Constants *)
+    (* expr = expr /. xx_Integer/; xx!=-1 :> ToReal[xx]; *)
+    expr = expr /. xx_Integer -> ToReal[xx];
+    expr = expr /. xx_Real -> ToReal[xx];
+
+    removeToRealRules =
+    {-ToReal[xx_] -> ToReal[- xx],
+     ToReal[xx_] + ToReal[yy_] -> ToReal[xx + yy],
+     ToReal[xx_] - ToReal[yy_] -> ToReal[xx - yy],
+     ToReal[xx_] * ToReal[yy_] -> ToReal[xx * yy],
+     ToReal[xx_] == ToReal[yy_] -> ToReal[xx == yy],
+     ToReal[xx_] != ToReal[yy_] -> ToReal[xx != yy],
+     pow[xx_, ToReal[power_]] -> pow[xx, power]};
+
+    expr = expr //. removeToRealRules;
+
+    (* Replace all operators and functions *)
+    (* kneg, kadd, ksub, kmul, kdiv *)
+    isNotKneg[n_] := ! MatchQ[n,kneg[_]];
+
+    arithRules =
+    { - xx_ -> kneg[xx],
+      xx_ * yy_ -> kmul[xx,yy],
+      xx_ / yy_ -> kdiv[xx,yy],
+      xx_ + yy_ -> kadd[xx,yy],
+      xx_ - yy_ -> ksub[xx,yy],
+      kmul[-1,xx_]            -> kneg[xx],
+      kmul[xx_,-1]            -> kneg[xx],
+      kmul[ToReal[-1],xx_]    -> kneg[xx],
+      kmul[xx_,ToReal[-1]]    -> kneg[xx],
+      ToReal[- xx_]           -> kneg[ToReal[xx]],
+      (* kmul[xx_,INV[yy_]]      -> kdiv[xx,yy], *)
+      (* kmul[INV[xx_],yy_]      -> kdiv[yy,xx], *)
+      kdiv[xx_,kdiv[yy_,zz_]] -> kdiv[kmul[xx,zz],yy],
+      kdiv[kdiv[xx_,yy_],zz_] -> kdiv[xx,kmul[yy,zz]],
+      kmul[kneg[xx_],yy_]     -> kneg[kmul[xx,yy]],
+      kmul[xx_,kneg[yy_]]     -> kneg[kmul[xx,yy]],
+      kdiv[kneg[xx_],yy_]     -> kneg[kdiv[xx,yy]],
+      kdiv[xx_,kneg[yy_]]     -> kneg[kdiv[xx,yy]],
+      kadd[kneg[xx_],yy_]     -> ksub[yy,xx],
+      ksub[kneg[xx_],yy_]     -> kneg[kadd[xx,yy]],
+      kadd[xx_,kneg[yy_]]     -> ksub[xx,yy],
+      ksub[xx_,kneg[yy_]]     -> kadd[xx,yy],
+      kneg[ksub[xx_,yy_]]     -> ksub[yy,xx],
+      Abs[xx_]      -> kfabs[xx],
+      Cos[xx_]      -> kcos[xx],
+      Log[xx_]      -> klog[xx],
+      Sin[xx_]      -> ksin[xx],
+      Tan[xx_]      -> ktan[xx],
+      exp[xx_]      -> kexp[xx],
+      fabs[xx_]     -> kfabs[xx],
+      fmax[xx_,yy_] -> kfmax[xx,yy],
+      fmin[xx_,yy_] -> kfmin[xx,yy],
+      log[xx_]      -> klog[xx],
+      pow[xx_,yy_]  -> kpow[xx,yy],
+      sqrt[xx_]     -> ksqrt[xx],
+      kcos[kneg[xx_]]   -> kcos[xx],
+      kfabs[kneg[xx_]]  -> kfabs[xx],
+      kfnabs[kneg[xx_]] -> kfnabs[xx],
+      kneg[kfabs[xx_]]  -> kfnabs[xx],
+      kneg[kfnabs[xx_]] -> kfabs[xx],
+      kneg[kneg[xx_]]   -> xx,
+      ksin[kneg[xx_]]   -> kneg[ksin[xx]],
+      ktan[kneg[xx_]]   -> kneg[ktan[xx]]};
+    expr = expr //. arithRules;
+
+    (* Undo some transformations *)
+    undoRules = 
+    { IfThen[_, aa_, aa_] -> aa,
+      IfThen[xx_?IntegerQ, aa_, bb_] /; xx!=0 :> aa,
+      IfThen[xx_?IntegerQ, aa_, bb_] /; xx==0 :> bb,
+      IfThen[kmul[xx_,yy_], aa_, bb_] -> IfThen[xx*yy, aa, bb],
+      IfThen[kmul[xx_,yy_] != zz_, aa_, bb_] -> IfThen[xx*yy!=zz, aa, bb],
+      IfThen[ToReal[xx_], aa_, bb_] -> IfThen[xx, aa, bb],
+      Scalar[kmul[xx_,yy_]] -> Scalar[xx*yy],
+      Scalar[kmul[xx_,yy_] != zz_] -> Scalar[xx*yy!=zz],
+      Scalar[ToReal[xx_]] -> Scalar[xx],
+      Scalar[xx_ != ToReal[yy_]] -> Scalar[xx != yy],
+      ToReal[kneg[xx_]] -> ToReal[-xx],
+      ToReal[kadd[xx_,yy_]] -> ToReal[xx+yy],
+      ToReal[ksub[xx_,yy_]] -> ToReal[xx-yy],
+      ToReal[kmul[xx_,yy_]] -> ToReal[xx*yy],
+      ToReal[xx_*kadd[yy_,zz_]] -> ToReal[xx*(yy+zz)],
+      kpow[xx_, kneg[power_]] -> kpow[xx, -power]};
+    expr = expr //. undoRules;
+    
+    (* FMA (fused multiply-add) instructions *)
+    (* kmadd (x,y,z) =   xy+z
+       kmsub (x,y,z) =   xy-z
+       knmadd(x,y,z) = -(xy+z)
+       knmsub(x,y,z) = -(xy-z) *)
+    fmaRules =
+    { kadd[kmul[xx_,yy_],zz_] -> kmadd[xx,yy,zz],
+      kadd[zz_,kmul[xx_,yy_]] -> kmadd[xx,yy,zz],
+      ksub[kmul[xx_,yy_],zz_] -> kmsub[xx,yy,zz],
+      ksub[zz_,kmul[xx_,yy_]] -> knmsub[xx,yy,zz],
+      kneg[kmadd [xx_,yy_,zz_]] -> knmadd[xx,yy,zz],
+      kneg[kmsub [xx_,yy_,zz_]] -> knmsub[xx,yy,zz],
+      kneg[knmadd[xx_,yy_,zz_]] -> kmadd [xx,yy,zz],
+      kneg[knmsub[xx_,yy_,zz_]] -> kmsub [xx,yy,zz]
+      (* we could match this and similar patterns
+         kmul[xx_, kadd[yy_, ToReal[+1]]] -> kmadd[xx, yy, xx],
+         kmul[xx_, kadd[yy_, ToReal[-1]]] -> kmsub[xx, yy, xx],
+         *)};
+    expr = expr //. fmaRules;
+    
+    Return[expr]]];
+
+(* Take an expression x and replace occurrences of Powers with the C
+  macros SQR, CUB, QAD *)
+DefFn[
+  ReplacePowers[expr_, vectorise:Boolean] :=
+  Module[
+    {rhs},
+    rhs = expr /. Power[xx_, -1] -> INV[xx];
+    If[SOURCELANGUAGE == "C",
+       {rhs = rhs /. Power[xx_,  2  ] -> SQR[xx];
+        rhs = rhs /. Power[xx_,  3  ] -> CUB[xx];
+        rhs = rhs /. Power[xx_,  4  ] -> QAD[xx];
+        rhs = rhs /. Power[xx_, -2  ] -> INV[SQR[xx]];
+        rhs = rhs /. Power[xx_,  1/2] -> sqrt[xx];
+        rhs = rhs /. Power[xx_, -1/2] -> INV[sqrt[xx]];
+        rhs = rhs /. Power[xx_,  0.5] -> sqrt[xx];
+        rhs = rhs /. Power[xx_, -0.5] -> INV[sqrt[xx]];
+        
+        (*
+           rhs = rhs /.  1/2 ->  khalf
+           rhs = rhs /. -1/2 -> -khalf;
+           
+           rhs = rhs /.  1/3 ->  kthird;
+           rhs = rhs /. -1/3 -> -kthird;
+           
+           rhs = rhs /.  2/3 ->  ktwothird;
+           rhs = rhs /. -2/3 -> -ktwothird;
+           
+           rhs = rhs /.  4/3 ->  kfourthird;
+           rhs = rhs /. -4/3 -> -kfourthird;
+           
+           rhs = rhs /.  8/3 ->  keightthird;
+           rhs = rhs /. -8/3 -> -keightthird;
+           *)
+        
+        (* Remove parentheses *)
+        rhs = rhs //. Parenthesis[xx_] -> xx;
+
+        (* Avoid rational numbers *)
+        rhs = rhs /. Rational[xx_,yy_] :> N[xx/yy, 30];
+        
+        rhs = rhs //. IfThen[cond1_,xx1_,yy1_] + IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> 
+        IfThen[cond1, Simplify[    xx1 +     xx2], Simplify[    yy1 +     yy2]];
+
+        rhs = rhs //. ff1_ IfThen[cond1_,xx1_,yy1_] + IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :>
+        IfThen[cond1, Simplify[ff1 xx1 +     xx2], Simplify[ff1 yy1 +     yy2]];
+
+        rhs = rhs //.      IfThen[cond1_,xx1_,yy1_] + ff2_ IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> 
+        IfThen[cond1, Simplify[    xx1 + ff2 xx2], Simplify[    yy1 + ff2 yy2]];
+
+        rhs = rhs //. ff1_ IfThen[cond1_,xx1_,yy1_] + ff2_ IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> 
+        IfThen[cond1, Simplify[ff1 xx1 + ff2 xx2], Simplify[ff1 yy1 + ff2 yy2]];
+
+        (* Is this still a good idea when FMA instructions are used? *)
+        rhs = rhs //. xx_ yy_ + xx_ zz_ -> xx (yy+zz);
+        rhs = rhs //. xx_ yy_ - xx_ zz_ -> xx (yy-zz);
+
+        rhs = rhs /. Power[E, power_] -> exp[power];
+
+        (* there have been some problems doing the Max/Min
+           replacement via the preprocessor for C, so we do it
+           here *)
+        (* Note: Mathematica simplifies Max[xx_] -> xx automatically *)
+        rhs = rhs //. Max[xx_, yy__] -> fmax[xx, Max[yy]];
+        rhs = rhs //. Min[xx_, yy__] -> fmin[xx, Min[yy]];
+
+        rhs = rhs /. Power[xx_, power_] -> pow[xx, power];
+
+        If[vectorise === True,
+           rhs = vectoriseExpression[rhs]];
+
+        (* Remove Scalar[] after vectorising *)
+        rhs = rhs /. Scalar[xx_] -> xx},
+       (* else *)
+       rhs = rhs /. Power[xx_, power_] -> xx^power];
+    (*       Print[rhs//FullForm];*)
+    rhs]];
+
+End[];
+
+EndPackage[];
diff --git a/Tools/CodeGen/Differencing.m b/Tools/CodeGen/Differencing.m
index 7b4e021..60e4b80 100644
--- a/Tools/CodeGen/Differencing.m
+++ b/Tools/CodeGen/Differencing.m
@@ -125,7 +125,7 @@ point. Should be checked by someone competent!
 
 *)
 
-BeginPackage["Differencing`", {"CodeGen`", "Kranc`", "MapLookup`", 
+BeginPackage["Differencing`", {"CodeGen`", "CodeGenC`", "CodeGenCactus`", "Kranc`", "MapLookup`", 
              (* "LinearAlgebra`MatrixManipulation`", *) "Errors`"}];
 
 CreateDifferencingHeader::usage = "";
@@ -153,26 +153,28 @@ DZero[n_] := (DPlus[n] + DMinus[n])/2;
 (* User API *)
 (*************************************************************)
 
-CreateDifferencingHeader[derivOps1_, zeroDims_, vectorise_, intParams_] :=
+DefFn[
+  CreateDifferencingHeader[derivOps1_, zeroDims_, vectorise_, intParams_] :=
   Module[{componentDerivOps, dupsRemoved, expressions, componentDerivOps2, zeroDimRules,
           pDefs, derivOps},
     derivOps = Flatten[Map[expandDerivOpOverParameters[#, intParams] &, derivOps1],1];
 
     Map[DerivativeOperatorVerify, derivOps];
-    zeroDimRules = Map[shift[#] -> 1 &, zeroDims];
 
-    componentDerivOps = Flatten[Map[DerivativeOperatorToComponents, derivOps]];
+    componentDerivOps = Flatten[Map[DerivativeOperatorToComponents[#,zeroDims] &,
+                                    derivOps]];
 
-    componentDerivOps2 = componentDerivOps /. zeroDimRules;
+    componentDerivOps2 = componentDerivOps;
 
     dupsRemoved = Block[{$RecursionLimit = Infinity}, RemoveDuplicateRules[componentDerivOps2]];
 
     mDefPairs = Map[ComponentDerivativeOperatorMacroDefinition[#, vectorise] &, dupsRemoved];
 
     pDefs = Union[Flatten[Map[First, mDefPairs]]];
+
     expressions = Flatten[Map[#[[2]]&, mDefPairs]];
 
-    {pDefs, Map[{#, "\n"} &, expressions]}];
+    {pDefs, Map[{#, "\n"} &, expressions]}]];
 
 ordergfds[_[v1_,___], _[v2_,___]] := 
   Order[v1,v2] != -1;
@@ -180,10 +182,11 @@ ordergfds[_[v1_,___], _[v2_,___]] :=
 getParamName[p_List] := lookup[p,Name];
 getParamName[p_] := p;
 
-PrecomputeDerivatives[derivOps_, expr_, intParams_] :=
+DefFn[
+  PrecomputeDerivatives[derivOps_, expr_, intParams_, zeroDims_] :=
   Module[{componentDerivOps, gfds, sortedgfds, opNames, intParamNames, paramsInOps,
           paramName, opsWithParam, opNamesWithParam, replace, param},
-    gfds = GridFunctionDerivativesInExpression[derivOps, expr];
+    gfds = GridFunctionDerivativesInExpression[derivOps, expr, zeroDims];
     sortedgfds = Sort[gfds, ordergfds];
 
     opNames = Union[Map[Head, sortedgfds]];
@@ -209,72 +212,78 @@ PrecomputeDerivatives[derivOps_, expr_, intParams_] :=
        "\n",
        SwitchStatement[paramName,
          Sequence@@Table[{value, Map[PrecomputeDerivative[# /. replace[value],#] &, sortedgfds]}, 
-                         {value, lookup[param, AllowedValues]}]]}]];
+                         {value, lookup[param, AllowedValues]}]]}]]];
 
-DeclareDerivatives[derivOps_, expr_] :=
+DefFn[
+  DeclareDerivatives[derivOps_, expr_] :=
   Module[{componentDerivOps, gfds, sortedgfds},
     Map[DerivativeOperatorVerify, derivOps];
     gfds = GridFunctionDerivativesInExpression[derivOps, expr];
     sortedgfds = Sort[gfds, ordergfds];
     {"/* Declare derivatives */\n",
-     Map[DeclareDerivative, sortedgfds]}];
+     Map[DeclareDerivative, sortedgfds]}]];
 
-ReplaceDerivatives[derivOps_, expr_, precompute_] :=
+DefFn[
+  ReplaceDerivatives[derivOps_, expr_, precompute_, zeroDims_] :=
   Module[{componentDerivOps, gfds},
     Map[DerivativeOperatorVerify, derivOps];
-    componentDerivOps = Flatten[Map[DerivativeOperatorToComponents, derivOps]];
-    gfds = GridFunctionDerivativesInExpression[derivOps, expr];
+    componentDerivOps = Flatten[Map[DerivativeOperatorToComponents[#,zeroDims] &, derivOps]];
+    gfds = GridFunctionDerivativesInExpression[derivOps, expr, zeroDims];
 
     If[precompute,
       rules = Map[# :> GridFunctionDerivativeName[#] &, gfds],
       rules = Map[# :> evaluateDerivative[#] &, gfds]];
-    expr /. rules];
+    expr /. rules]];
 
 (* Generate code to ensure that there are sufficient ghost and
    boundary points for the passed derivative operators used in eqs *)
-CheckStencil[derivOps_, eqs_, name_] :=
+DefFn[
+  CheckStencil[derivOps_, eqs_, name_, zeroDims_] :=
   Module[{gfds, rgzList, rgz},
-    gfds = Map[GridFunctionDerivativesInExpression[{#}, eqs] &, derivOps];
-    rgzList = MapThread[If[Length[#2] > 0, DerivativeOperatorStencilWidth[#1], {0,0,0}] &, {derivOps, gfds}];
+    gfds = Map[GridFunctionDerivativesInExpression[{#}, eqs, zeroDims] &, derivOps];
+    rgzList = MapThread[If[Length[#2] > 0, DerivativeOperatorStencilWidth[#1,zeroDims], {0,0,0}] &, {derivOps, gfds}];
     If[Length[rgzList] === 0, Return[{}]];
     rgz = Map[Max, Transpose[rgzList]];
     If[Max[rgz] == 0, {},
-    {"GenericFD_EnsureStencilFits(cctkGH, ", Quote@name, ", ", Riffle[rgz,", "], ");\n"}]];
+    {"GenericFD_EnsureStencilFits(cctkGH, ", Quote@name, ", ", Riffle[rgz,", "], ");\n"}]]];
 
 parametersUsedInOps[derivOps_, intParams_] :=
   Union@Flatten[Map[Cases[derivOps, getParamName[#] -> #, Infinity] &,
                     intParams], 1];
 
-CheckStencil[derivOps_, eqs_, name_, intParams_] :=
+DefFn[
+  CheckStencil[derivOps_, eqs_, name_, zeroDims_, intParams_] :=
   Module[{psUsed, p},
     psUsed = parametersUsedInOps[derivOps, intParams];
     If[Length[psUsed] > 1, Throw["Too many parameters in partial derivatives"]];
     If[psUsed === {},
-      CheckStencil[derivOps,eqs,name],
+      CheckStencil[derivOps,eqs,name,zeroDims],
       p = psUsed[[1]];
       SwitchStatement[getParamName[p],
         Sequence@@Table[{value,
                          CheckStencil[derivOps/.getParamName[p]->value, eqs,
-                                      name]},
-                        {value, lookup[p, AllowedValues]}]]]];
+                                      name, zeroDims]},
+                        {value, lookup[p, AllowedValues]}]]]]];
 
 (*************************************************************)
 (* Misc *)
 (*************************************************************)
 
-PrecomputeDerivative[d:pd_[gf_, inds___], vargfd_:Automatic] :=
+DefFn[
+  PrecomputeDerivative[d:pd_[gf_, inds___], vargfd_:Automatic] :=
   Module[{},
     If[vargfd === Automatic,
       DeclareAssignVariable[DataType[], GridFunctionDerivativeName[d], evaluateDerivative[d]],
-      AssignVariable[GridFunctionDerivativeName[vargfd], evaluateDerivative[d]]]];
+      AssignVariable[GridFunctionDerivativeName[vargfd], evaluateDerivative[d]]]]];
 
-evaluateDerivative[d:pd_[gf_, inds___]] :=
+DefFn[
+  evaluateDerivative[d:pd_[gf_, inds___]] :=
   Module[{macroname},
     macroName = ComponentDerivativeOperatorMacroName[pd[inds] -> expr];
     (* Return[ToString[macroName] <> "(" <> ToString[gf] <> ", i, j, k)"] *)
     (* Return[ToString[macroName] <> "(" <> ToString[gf] <> ")"] *)
     Return[ToString[macroName] <> "(&" <> ToString[gf] <> "[index])"]
-  ];
+  ]];
 
 DeclareDerivative[d:pd_[gf_, inds___]] :=
   DeclareVariable[GridFunctionDerivativeName[d], "// CCTK_REAL_VEC"];
@@ -284,37 +293,41 @@ DeclareDerivative[d:pd_[gf_, inds___]] :=
 (* GridFunctionDerivative *)
 (*************************************************************)
 
-GridFunctionDerivativeName[pd_[gf_, inds___]] :=
+DefFn[
+  GridFunctionDerivativeName[pd_[gf_, inds___]] :=
   Module[{},
     stringName = StringJoin[Map[ToString, Join[{pd}, {inds}, {gf}]]];
-    Symbol["Global`" <> stringName]];
+    Symbol["Global`" <> stringName]]];
 
 
-GridFunctionDerivativesInExpression[derivOps_, expr_] := 
+DefFn[
+  GridFunctionDerivativesInExpression[derivOps_, expr_, zeroDims_] := 
   Module[{componentDerivOps, derivs, patterns, dupsRemoved},
-    componentDerivOps = Flatten[Map[DerivativeOperatorToComponents, derivOps]];
+    componentDerivOps = Flatten[Map[DerivativeOperatorToComponents[#,zeroDims]&, derivOps]];
     dupsRemoved = RemoveDuplicateRules[componentDerivOps];
     derivs = Map[First, dupsRemoved];
     patterns = Map[# /. x_[inds___] -> x[y_, inds] &, derivs];
-    Flatten[Map[Union[Cases[{expr}, #, Infinity]] &, patterns]]];
+    Flatten[Map[Union[Cases[{expr}, #, Infinity]] &, patterns]]]];
 
 (* Return the definition associated with a grid function derivative *)
-GridFunctionDerivativeToDef[pd_[gf_, inds___], derivOps_] :=
+GridFunctionDerivativeToDef[pd_[gf_, inds___], derivOps_, zeroDims] :=
   Module[{componentDerivOps},
-    componentDerivOps = Flatten[Map[DerivativeOperatorToComponents, derivOps]];
+    componentDerivOps = Flatten[Map[DerivativeOperatorToComponents[#,zeroDims]&, derivOps]];
     pd[inds] /. componentDerivOps];
 
 (*************************************************************)
 (* DerivativeOperator *)
 (*************************************************************)
 
-sbpMacroDefinition[macroName_, d_] :=
+DefFn[
+  sbpMacroDefinition[macroName_, d_] :=
   Module[{ds = Switch[d, 1, "x", 2, "y", 3, "z"],
           l = Switch[d, 1, "i", 2, "j", 3, "k"]},
     FlattenBlock[{"#define ", macroName, "(u,i,j,k) (sbp_deriv_" <> ds
-    <> "(i,j,k,sbp_" <> l <> "min,sbp_" <> l <> "max,d" <> ds <> ",u,q" <> ds <> ",cctkGH))"}]    ];
+    <> "(i,j,k,sbp_" <> l <> "min,sbp_" <> l <> "max,d" <> ds <> ",u,q" <> ds <> ",cctkGH))"}]    ]];
 
-ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> expr_), vectorise_] :=
+DefFn[
+  ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> expr_), vectorise_] :=
   Module[{macroName, rhs, i = "i", j = "j", k = "k", spacings, spacings2, pat, ss, num, den, newnum, signModifier, quotient, liName, finalDef},
   
     macroName = ComponentDerivativeOperatorMacroName[componentDerivOp];
@@ -369,15 +382,16 @@ ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> e
        Print["Sequenced: ", Apply[SequenceForm,Simplify[1/(ss /. spacings2)],{0,Infinity}]];*)
 
        liName = "p" <> signModifier <> quotient <> ToString[Apply[SequenceForm,Simplify[1/(ss /. spacings2)],{0,Infinity}]];
-(*       Print["liName == ", liName];*)
+       pDefs = {{liName -> CFormHideStrings[ReplacePowers[num / den ss /. spacings2, vectorise]]}};
 
        rhs = rhs /. pat -> Times[liName, rest],
 (*       Print["!!!!!!!!DOES NOT MATCH!!!!!!!!!"];*)
-       rhs = rhs];
+       rhs = rhs;
+       pDefs = {};
+       liName = rhs];
 
 (*    Print["rhs3 == ", FullForm[rhs]];*)
 
-    pDefs = {{liName -> CFormHideStrings[ReplacePowers[num / den ss /. spacings2, vectorise]]}};
 
 (*    rhs = Factor[rhs];*)
     rhs = rhs //. (x_ a_ + x_ b_) -> x (a+b);
@@ -442,17 +456,19 @@ ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> e
     }]}];
 
     finalDef
-];
+]];
 
 ComponentDerivativeOperatorMacroName[componentDerivOp:(name_[inds___] -> expr_)] :=
   Module[{stringName},
     stringName = StringJoin[Map[ToString, Join[{name}, {inds}]]];
     stringName];
 
-DerivativeOperatorStencilWidth[derivOp_] :=
-  Map[Max, Transpose[Map[ComponentDerivativeOperatorStencilWidth, DerivativeOperatorToComponents[derivOp]]]];
+DerivativeOperatorStencilWidth[derivOp_,zeroDims_] :=
+  Map[Max, Transpose[Map[ComponentDerivativeOperatorStencilWidth, 
+                         DerivativeOperatorToComponents[derivOp,zeroDims]]]];
 
-ComponentDerivativeOperatorStencilWidth[componentDerivOp:(name_[inds___] -> expr_)] :=
+DefFn[
+  ComponentDerivativeOperatorStencilWidth[componentDerivOp:(name_[inds___] -> expr_)] :=
   Module[{cases, nx, ny, nz, result},
     result = Table[
     cases = Union[Flatten[Cases[{expr}, shift[d] | Power[shift[d],_], Infinity]]];
@@ -468,7 +484,7 @@ ComponentDerivativeOperatorStencilWidth[componentDerivOp:(name_[inds___] -> expr
 
     If[!And@@Map[NumericQ, result],
       Throw["Stencil width is not numeric in "<>ToString[componentDerivOp]]];
-    result];
+    result]];
 
 (* Farm out each term of a difference operator *)
 DifferenceGF[op_, i_, j_, k_, vectorise_] :=
@@ -477,7 +493,7 @@ DifferenceGF[op_, i_, j_, k_, vectorise_] :=
     
     If[Head[expanded] === Plus,
       Apply[Plus, Map[DifferenceGFTerm[#, i, j, k, vectorise] &, expanded]],
-      DifferenceGFTerm[expanded, i, j, k]]];
+      DifferenceGFTerm[expanded, i, j, k, vectorise]]];
 
 
 (* Return the fragment of a macro definition for defining a derivative
@@ -500,23 +516,18 @@ DifferenceGFTerm[op_, i_, j_, k_, vectorise_] :=
     If[Cases[{remaining}, shift[_], Infinity] != {},
       ThrowError["Could not parse difference operator", op]];
     
-    If[CodeGen`SOURCELANGUAGE == "C",
+    If[CodeGenC`SOURCELANGUAGE == "C",
 
   If[vectorise,
-    remaining "vec_loadu_maybe3" <>
-      "(" <> ToString[CFormHideStrings[nx /. {dir1->1, dir2->1, dir3->1}]] <> "," <>
-             ToString[CFormHideStrings[ny /. {dir1->1, dir2->1, dir3->1}]] <> "," <>
-             ToString[CFormHideStrings[nz /. {dir1->1, dir2->1, dir3->1}]] <> "," <>
-      "*(CCTK_REAL const*)&((char const*)(u))" <>
-        "[cdi*(" <> ToString[CFormHideStrings[nx]] <> ")" <>
-        "+cdj*(" <> ToString[CFormHideStrings[ny]] <> ")" <>
-        "+cdk*(" <> ToString[CFormHideStrings[nz]] <> ")])",
-
-    remaining
-      "(*(CCTK_REAL const*)&((char const*)(u))" <>
-        "[cdi*(" <> ToString[CFormHideStrings[nx]] <> ")" <>
-        "+cdj*(" <> ToString[CFormHideStrings[ny]] <> ")" <>
-        "+cdk*(" <> ToString[CFormHideStrings[nz]] <> ")])"],
+    remaining "KRANC_GFOFFSET3D(u," <>
+      ToString[CFormHideStrings[nx /. {dir1->1, dir2->1, dir3->1}]] <> "," <>
+      ToString[CFormHideStrings[ny /. {dir1->1, dir2->1, dir3->1}]] <> "," <>
+      ToString[CFormHideStrings[nz /. {dir1->1, dir2->1, dir3->1}]] <> ")",
+
+    remaining "KRANC_GFOFFSET3D(u," <>
+      ToString[CFormHideStrings[nx]] <> "," <>
+      ToString[CFormHideStrings[ny]] <> "," <>
+      ToString[CFormHideStrings[nz]] <> ")"],
 
     remaining "u(" <> ToString[FortranForm[i+nx]] <> "," <> 
       ToString[FortranForm[j+ny]] <> "," <> ToString[FortranForm[k+nz]] <> ")"] ];
@@ -524,24 +535,27 @@ DifferenceGFTerm[op_, i_, j_, k_, vectorise_] :=
 
 DerivativeOperatorGFDs[gf_];
 
-DerivativeOperatorToComponents[name_[indPatterns___] -> expr_] :=
-  Module[{ips, symbols, symbolRanges, symbolLHS, table},
+DefFn[
+  DerivativeOperatorToComponents[name_[indPatterns___] -> expr_, zeroDims_] :=
+  Module[{ips, symbols, symbolRanges, symbolLHS, table, zeroDimRules},
     ips = {indPatterns};
 
+    zeroDimRules = Map[shift[#] -> 1 &, zeroDims];
+
     If[MatchQ[ips, List[ (_Pattern) ...]],
 
       symbols = Map[First, ips];
       symbolRanges = Map[{#, 1, 3} &, Union[symbols]];
       symbolLHS = name[Apply[Sequence, symbols]];
       table = Apply[Table, Join[{symbolLHS -> expr}, symbolRanges]];
-      Return[Flatten[table]]];
+      Return[Flatten[table/.zeroDimRules]]];
 
 
     If[MatchQ[ips, List[ (_ ? NumberQ) ...]],
-      Return[{name[indPatterns] -> expr}]];
+      Return[{name[indPatterns] -> expr/.zeroDimRules}]];
 
     Throw["DerivativeOperatorToComponents: Expecting indices which are symbolic patterns or numbers"];
-];
+]];
 
 DerivativeOperatorVerify[derivOp_] :=
   If[!MatchQ[derivOp, pd_[_Pattern ...] -> expr_?DerivativeOperatorRHSVerify] && 
@@ -555,7 +569,8 @@ DerivativeOperatorRHSVerify[expr_] :=
     True];
 
 
-RemoveDuplicates[l_] :=
+DefFn[
+  RemoveDuplicates[l_] :=
   Module[{this,next,rest,positions},
     If[l === {},
       Return[{}]];
@@ -566,7 +581,7 @@ RemoveDuplicates[l_] :=
 
        positions = Position[rest, this];
        next = Delete[rest, positions];
-       Prepend[RemoveDuplicates[next], this]]];
+       Prepend[RemoveDuplicates[next], this]]]];
 
 RemoveDuplicateRules[l_] :=
   Module[{lhs,lhs2,rhs2,result},
diff --git a/Tools/CodeGen/Interface.m b/Tools/CodeGen/Interface.m
index 5f8f613..c547b58 100644
--- a/Tools/CodeGen/Interface.m
+++ b/Tools/CodeGen/Interface.m
@@ -135,6 +135,7 @@ CreateKrancInterface[nonevolvedGroups_, evolvedGroups_, rhsGroups_, groups_,
     interface = CreateInterface[implementation, inheritedImplementations, 
       Join[includeFiles, {CactusBoundary`GetIncludeFiles[]},
            If[OptionValue[UseLoopControl], {"loopcontrol.h"}, {}],
+           If[OptionValue[UseOpenCL], {"OpenCLRunTime.h"}, {}],
            If[OptionValue[UseVectors], {"vectors.h"}, {}]],
       groupStructures,
       UsesFunctions ->
diff --git a/Tools/CodeGen/Jacobian.m b/Tools/CodeGen/Jacobian.m
index 20f24df..bd72845 100644
--- a/Tools/CodeGen/Jacobian.m
+++ b/Tools/CodeGen/Jacobian.m
@@ -18,7 +18,8 @@
     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *)
 
-BeginPackage["Jacobian`", {"Errors`", "Helpers`", "Kranc`", "Differencing`", "MapLookup`", "CodeGen`", "KrancGroups`"}];
+BeginPackage["Jacobian`", {"Errors`", "Helpers`", "Kranc`", "Differencing`", "MapLookup`",
+                           "CodeGen`", "CodeGenC`", "KrancGroups`"}];
 
 JacobianQ;
 InsertJacobian;
@@ -69,11 +70,11 @@ shorthandsInDerivDef[def_, shorthands_] :=
 
 (* Insert a Jacobian shorthand definition into a list of equations
    after all the shorthands that it uses *)
-insertDerivInEqs[deriv_, defs_, eqs_, shorthands_] :=
+insertDerivInEqs[deriv_, defs_, eqs_, shorthands_, zeroDims_] :=
   Module[{shortsUsed, lhss, positions, positions2, position, newShort, derivsInShort, defsUsed},
     newShort = jacobianShorthand[deriv]; (* Definition expression for new shorthand *)
-    derivsInShort = GridFunctionDerivativesInExpression[defs, newShort]; (* Derivatives needed *)
-    defsUsed = Map[GridFunctionDerivativeToDef[#, defs] &, derivsInShort]; (* Definitions of those derivatives *)
+    derivsInShort = GridFunctionDerivativesInExpression[defs, newShort, zeroDims]; (* Derivatives needed *)
+    defsUsed = Map[GridFunctionDerivativeToDef[#, defs, zeroDims] &, derivsInShort]; (* Definitions of those derivatives *)
     shortsUsed = Union[Flatten[Map[shorthandsInDerivDef[#, shorthands] &, defsUsed],1]];
     lhss = Map[First, eqs];
     positions = Map[Position[lhss, #, {1}, 1] &, shortsUsed];
@@ -94,7 +95,8 @@ InsertJacobian[calc_List, opts:OptionsPattern[]] :=
           shorts},
     shorts = lookupDefault[calc, Shorthands, {}];
     pdDefs = OptionValue[PartialDerivatives];
-    derivs = GridFunctionDerivativesInExpression[pdDefs, lookup[calc, Equations]];
+    derivs = GridFunctionDerivativesInExpression[pdDefs, lookup[calc, Equations],
+                                                 OptionValue[ZeroDimensions]];
     If[Length[derivs] === 0, Return[calc]];
     newShortDefs = Map[jacobianShorthand, derivs];
     newShorts = Map[First, newShortDefs];
@@ -104,7 +106,8 @@ InsertJacobian[calc_List, opts:OptionsPattern[]] :=
     combinedEqs = newEqs;
 
     Do[
-      combinedEqs = insertDerivInEqs[derivs[[i]], pdDefs, combinedEqs, shorts],
+      combinedEqs = insertDerivInEqs[derivs[[i]], pdDefs, combinedEqs, shorts,
+                                     OptionValue[ZeroDimensions]],
       {i, Length[derivs], 1, -1}];
 
     combinedCalc = mapReplace[mapReplace[mapEnsureKey[calc, Shorthands, {}], Shorthands, combinedShorts], Equations, combinedEqs];
diff --git a/Tools/CodeGen/Kranc.m b/Tools/CodeGen/Kranc.m
index 3606a74..4b47c52 100644
--- a/Tools/CodeGen/Kranc.m
+++ b/Tools/CodeGen/Kranc.m
@@ -22,10 +22,11 @@ BeginPackage["Kranc`"];
 
 (* CodeGen.m *)
 
-{INV, SQR, CUB, QAD, IfThen, Scalar, ToReal, sqrt, exp, pow, fmax, fmin,
+{INV, SQR, CUB, QAD, IfThen, Parenthesis, Scalar, ToReal,
+ sqrt, exp, pow, fmax, fmin,
  kmadd, kmsub, knmadd, knmsub, kpos, kneg, kadd, ksub, kmul, kdiv,
- kfabs, kfmax, kfmin, ksqrt, kexp, klog, kpow,
- dir1, dir2, dir3, dt, dx, dy, dz,
+ kcos, kfabs, kfmax, kfmin, ksqrt, kexp, klog, kpow, ksin, ktan,
+ dir1, dir2, dir3, dt, dx, dy, dz, t,
  khalf, kthird, ktwothird, kfourthird, keightthird};
 
 (* Helpers.m *)
@@ -38,7 +39,7 @@ dummy;
 LoopPreIncludes, GroupImplementations, PartialDerivatives, NoSimplify,
 Boundary, Interior, InteriorNoSync, Where, AddToStencilWidth,
 Everywhere, normal1, normal2, normal3, INV, SQR, CUB, QAD, dot, pow,
-exp, dt, dx, dy, dz, idx, idy, idz}
+exp, dt, dx, dy, dz, idx, idy, idz, t}
 
 {ConditionalOnKeyword, ConditionalOnKeywords, CollectList, Interior,
 InteriorNoSync, Boundary, BoundaryWithGhosts, Where, PreDefinitions,
@@ -70,6 +71,7 @@ ThornOptions =
   ReflectionSymmetries -> {},
   ZeroDimensions -> {},
   UseLoopControl -> False,
+  UseOpenCL -> False,
   UseVectors -> False,
   ProhibitAssignmentToGridFunctionsRead -> False,
   IncludeFiles -> {},
diff --git a/Tools/CodeGen/KrancThorn.m b/Tools/CodeGen/KrancThorn.m
index 44ee16f..26ca6ef 100644
--- a/Tools/CodeGen/KrancThorn.m
+++ b/Tools/CodeGen/KrancThorn.m
@@ -208,16 +208,22 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[
     InfoMessage[Terse, "Creating differencing header file"];
     {pDefs, diffHeader} = CreateDifferencingHeader[partialDerivs, OptionValue[ZeroDimensions], OptionValue[UseVectors], OptionValue[IntParameters]];
     diffHeader = Join[
-        If[OptionValue[UseVectors], {"#include <assert.h>\n",
-                                     "#include \"vectors.h\"\n",
-                                     "\n"},
-                                    {}],
+        If[OptionValue[UseVectors] && ! OptionValue[UseOpenCL],
+           {"#include <assert.h>\n",
+            "#include \"vectors.h\"\n",
+            "\n"},
+           {}],
         diffHeader];
+    diffHeader = If[OptionValue[UseOpenCL],
+                    "static char const * const differencing =\n" <>
+                    Stringify[diffHeader] <>
+                    ";\n",
+                    diffHeader];
 
     (* Add the predefinitions into the calcs *)
     calcs = Map[Join[#, {PreDefinitions -> pDefs}] &, calcs];
 
-    ext = CodeGen`SOURCESUFFIX;
+    ext = CodeGenC`SOURCESUFFIX;
 
     (* Construct a source file for each calculation *)
     allParams = Join[Map[ParamName, realParamDefs],
@@ -249,7 +255,7 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[
                   {Filename -> "Startup.cc", Contents -> startup}, 
                   {Filename -> "RegisterMoL.cc", Contents -> molregister},
                   {Filename -> "RegisterSymmetries.cc", Contents -> symregister},
-                 {Filename -> "Differencing.h", Contents -> diffHeader}},
+                  {Filename -> "Differencing.h", Contents -> diffHeader}},
                   MapThread[{Filename -> #1, Contents -> #2} &, 
                             {calcFilenames, calcSources}], boundarySources]};
     InfoMessage[Terse, "Creating thorn"];
diff --git a/Tools/CodeGen/Schedule.m b/Tools/CodeGen/Schedule.m
index 66dec5a..f596ad1 100644
--- a/Tools/CodeGen/Schedule.m
+++ b/Tools/CodeGen/Schedule.m
@@ -109,7 +109,7 @@ scheduleCalc[calc_, groups_] :=
         SynchronizedGroups -> If[StringMatchQ[#, "*MoL_CalcRHS*", IgnoreCase -> True] || StringMatchQ[#, "*MoL_RHSBoundaries*", IgnoreCase -> True],
                                  {},
                                  groupsToSync],
-        Language           -> CodeGen`SOURCELANGUAGE, 
+        Language           -> CodeGenC`SOURCELANGUAGE, 
         Comment            -> lookup[calc, Name]
       },
        If[triggered, {TriggerGroups -> lookup[calc, TriggerGroups]},
@@ -147,7 +147,7 @@ scheduleCalc[calc_, groups_] :=
       fnSched = {
         Name               -> lookup[calc, Name],
         SchedulePoint      -> "in " <> groupName,
-        Language           -> CodeGen`SOURCELANGUAGE,
+        Language           -> CodeGenC`SOURCELANGUAGE,
         Comment            -> lookup[calc, Name]
       };
 
@@ -165,7 +165,7 @@ scheduleCalc[calc_, groups_] :=
         SchedulePoint      -> "in " <> bcGroupName,
         SynchronizedGroups -> groupsToSync,
         Options               -> "level",
-        Language           -> CodeGen`SOURCELANGUAGE,
+        Language           -> CodeGenC`SOURCELANGUAGE,
         Comment            -> lookup[calc, Name] <> "_SelectBCs"
       };
 
diff --git a/Tools/CodeGen/Thorn.m b/Tools/CodeGen/Thorn.m
index 8033d2a..afffb47 100644
--- a/Tools/CodeGen/Thorn.m
+++ b/Tools/CodeGen/Thorn.m
@@ -23,7 +23,7 @@
 (* This package provides a set of functions to create the various
    parts of a Cactus thorn and assemble them. *)
 
-BeginPackage["Thorn`", "CodeGen`", "CalculationFunction`",
+BeginPackage["Thorn`", "CodeGen`", "CodeGenC`", "CodeGenCactus`", "CalculationFunction`",
   "CalculationBoundaries`", "MapLookup`", "KrancGroups`", "Helpers`",
   "Errors`", "Kranc`"];
 
@@ -221,6 +221,7 @@ CreateConfiguration[opts:OptionsPattern[]] :=
   {whoWhen["CCL"],
    "REQUIRES GenericFD\n",
    If[OptionValue[UseLoopControl], "REQUIRES LoopControl\n", {}],
+   If[OptionValue[UseOpenCL], "REQUIRES OpenCL OpenCLRunTime\n", {}],
    If[OptionValue[UseVectors], "REQUIRES Vectors\n", {}]
   };
 
@@ -501,11 +502,11 @@ CreateSetterSource[calcs_, debug_, include_, imp_,
   If[!MatchQ[include, _List],
     Throw["CreateSetterSource: Include should be a list but is in fact " <> ToString[include]]];
 
-  {whoWhen[CodeGen`SOURCELANGUAGE],
+  {whoWhen[CodeGenC`SOURCELANGUAGE],
 
-   "#define KRANC_" <> ToUpperCase[CodeGen`SOURCELANGUAGE] <> "\n\n",
+   "#define KRANC_" <> ToUpperCase[CodeGenC`SOURCELANGUAGE] <> "\n\n",
 
-   If[CodeGen`SOURCELANGUAGE == "C",
+   If[CodeGenC`SOURCELANGUAGE == "C",
          {IncludeSystemFile["assert.h"],
           IncludeSystemFile["math.h"],
           IncludeSystemFile["stdio.h"],
@@ -518,6 +519,7 @@ CreateSetterSource[calcs_, debug_, include_, imp_,
                          (*"precomputations.h",*) "GenericFD.h", "Differencing.h"},
                          include,
                          If[OptionValue[UseLoopControl], {"loopcontrol.h"}, {}],
+                         If[OptionValue[UseOpenCL], {"OpenCLRunTime.h"}, {}],
                          If[OptionValue[UseVectors], {"vectors.h"}, {}]]],
    calculationMacros[OptionValue[UseVectors]],
 
@@ -591,8 +593,8 @@ CreateSymmetriesRegistrationSource[thornName_, implementationName_, GFs_, reflec
       Print["Registering Symmetries for: ", GFs];
     ];
 
-  lang = CodeGen`SOURCELANGUAGE;
-  CodeGen`SOURCELANGUAGE = "C";
+  lang = CodeGenC`SOURCELANGUAGE;
+  CodeGenC`SOURCELANGUAGE = "C";
 
   spec = Table[{FullName -> implementationName <> "::" <> ToString@GFs[[j]],
                 Sym      -> If[reflectionSymmetries === False,
@@ -615,7 +617,7 @@ CreateSymmetriesRegistrationSource[thornName_, implementationName_, GFs_, reflec
 ]
   };
 
-  CodeGen`SOURCELANGUAGE = lang;
+  CodeGenC`SOURCELANGUAGE = lang;
 
 tmp
 ];
@@ -642,8 +644,8 @@ CreateMoLRegistrationSource[spec_, debug_] :=
     Print["  Primitive Gridfunctions: ", lookup[spec, PrimitiveGFs] ];
     ];
 	
-    lang = CodeGen`SOURCELANGUAGE;
-    CodeGen`SOURCELANGUAGE= "C";
+    lang = CodeGenC`SOURCELANGUAGE;
+    CodeGenC`SOURCELANGUAGE= "C";
 
     tmp = {whoWhen["C"],
 
@@ -674,7 +676,7 @@ CreateMoLRegistrationSource[spec_, debug_] :=
           lookup[spec, PrimitiveGFs]]],  *)
 	"return;\n"}]};
 
-      CodeGen`SOURCELANGUAGE = lang;
+      CodeGenC`SOURCELANGUAGE = lang;
 
 tmp
 ];
@@ -907,8 +909,8 @@ CreateMoLBoundariesSource[spec_] :=
      "\n}\n"}];
 
 
-   lang = CodeGen`SOURCELANGUAGE;
-   CodeGen`SOURCELANGUAGE = "C";
+   lang = CodeGenC`SOURCELANGUAGE;
+   CodeGenC`SOURCELANGUAGE = "C";
 
   tmp = {whoWhen["C"],
 
@@ -958,7 +960,7 @@ CreateMoLBoundariesSource[spec_] :=
      "*/\n\n"
      };
 
-	 CodeGen`SOURCELANGUAGE = lang;
+	 CodeGenC`SOURCELANGUAGE = lang;
 tmp
 ];
 
@@ -970,8 +972,8 @@ CreateMoLExcisionSource[spec_] :=
  
   Print["Applying excision to GFs: ", gfs];
  
-  currentlang = CodeGen`SOURCELANGUAGE;
-  CodeGen`SOURCELANGUAGE = "Fortran";
+  currentlang = CodeGenC`SOURCELANGUAGE;
+  CodeGenC`SOURCELANGUAGE = "Fortran";
 
   excisionExtrap[gf_] :=  "  call ExcisionExtrapolate(ierr, "
     <> ToString@gf <> ", " <> ToString@gf
@@ -1055,7 +1057,7 @@ CreateMoLExcisionSource[spec_] :=
      "return\n"}]
 };
 
-CodeGen`SOURCELANGUAGE = currentlang;
+CodeGenC`SOURCELANGUAGE = currentlang;
 
 body
 ];
@@ -1242,8 +1244,8 @@ CreateMPCharSource[spec_, debug_] :=
 
   groups = Map[unqualifiedGroupName, lookup[spec, Groups]];
 
-  lang = CodeGen`SOURCELANGUAGE;
-  CodeGen`SOURCELANGUAGE = "C";
+  lang = CodeGenC`SOURCELANGUAGE;
+  CodeGenC`SOURCELANGUAGE = "C";
 
   tmp = {whoWhen["C"],
 
@@ -1305,7 +1307,7 @@ charInfoFunction["C2P", spec, debug],      "\n\n",
 charInfoFunction["WHATEVER", spec, debug]
 };
 
-CodeGen`SOURCELANGUAGE = lang;
+CodeGenC`SOURCELANGUAGE = lang;
 tmp
 ];
 
@@ -1332,8 +1334,8 @@ CreatePrecompMacros[functions_] :=
 CreateStartupFile[thornName_, bannerText_] :=
   Module[{tmp, lang},
   
-  lang = CodeGen`SOURCELANGUAGE;
-  CodeGen`SOURCELANGUAGE = "C";
+  lang = CodeGenC`SOURCELANGUAGE;
+  CodeGenC`SOURCELANGUAGE = "C";
 
   tmp = {whoWhen["C"],
 
@@ -1343,7 +1345,7 @@ CreateStartupFile[thornName_, bannerText_] :=
       "CCTK_RegisterBanner(banner);\n",
       "return 0;\n"}]};
 
-  CodeGen`SOURCELANGUAGE = lang;
+  CodeGenC`SOURCELANGUAGE = lang;
 
   tmp
    ];
diff --git a/Tools/MathematicaMisc/Errors.m b/Tools/MathematicaMisc/Errors.m
index 4de9a65..6032b0a 100644
--- a/Tools/MathematicaMisc/Errors.m
+++ b/Tools/MathematicaMisc/Errors.m
@@ -1,5 +1,5 @@
 
-BeginPackage["Errors`"];
+BeginPackage["Errors`", {"Profile`"}];
 
 PrintError::usage = "";
 ThrowError::usage = "";
@@ -9,12 +9,14 @@ VerifyStringList;
 VerifyList;
 InfoMessage;
 SetDebugLevel;
+ErrorDefinition::usage = "ErrorDefinition[f] creates a default definition of a function f which throws an exception.  This can be used to catch programming errors where f is called with incorrect arguments.";
 
 DebugQuiet = 0;
 Warnings = 1
 Terse = 2;
 Info = 3;
 InfoFull = 4;
+DefFn;
 
 Begin["`Private`"];
 
@@ -35,7 +37,7 @@ PrintStructure[x_]:=
 
 PrintStructure[l_List, prefix_, suffix_] :=
   Module[{},
-    If[StringLength[ToString[l] <> prefix] > 50,
+    If[StringLength[ToString[l,InputForm] <> prefix] > 50,
       Print[prefix, "{"];
       Map[PrintStructure[#, "  " <> prefix, ","] &, l];
       Print[prefix, "}"],
@@ -43,7 +45,7 @@ PrintStructure[l_List, prefix_, suffix_] :=
       Print[prefix, ToString[l,InputForm]]]];
 
 PrintStructure[s_, prefix_, suffix_] :=
-  Print[prefix, s, suffix];
+  Print[prefix, s//InputForm, suffix];
 
 PrintError[err_] :=
   Module[{},
@@ -82,12 +84,24 @@ VerifyList[l_] :=
 InfoMessage[level_, message__] :=
   Module[{args = {message}},
     If[level <= debugLevel,
-      Map[PrintStructure, args]];
+      Map[Print, args]];
   ];
 
 SetDebugLevel[level_] :=
   debugLevel = level;
 
+ErrorDefinition[x_] :=
+  x[args___] :=
+    ThrowError["Invalid arguments to "<>ToString[x], {args}//FullForm];
+
+SetAttributes[DefFn, HoldAll];
+
+DefFn[def:(fn_[args___] := body_)] :=
+  Module[
+    {},
+    ErrorDefinition[fn];
+    fn[args] := (*Profile[fn,*)body(*]*)];
+
 End[];
 
 EndPackage[];
diff --git a/Tools/MathematicaMisc/MapLookup.m b/Tools/MathematicaMisc/MapLookup.m
index 1711cdf..ee0c85b 100644
--- a/Tools/MathematicaMisc/MapLookup.m
+++ b/Tools/MathematicaMisc/MapLookup.m
@@ -4,6 +4,7 @@ BeginPackage["MapLookup`", {"Errors`"}];
 lookup::usage = "";
 mapContains::usage = "";
 mapReplace::usage = "";
+mapReplaceAdd::usage = "";
 mapValueMap::usage = "";
 lookupDefault::usage = "";
 mapValueMapMultiple::usage = "";
@@ -59,6 +60,13 @@ mapReplace[map_, key_, value_] :=
     VerifyMap[map];
     Map[If[First[#] === key, key -> value, #] &, map]];
 
+mapReplaceAdd[map_, key_, value_] :=
+  Module[{},
+    VerifyMap[map];
+  If[mapContains[map, key],
+     mapReplace[map, key, value],
+     mapAdd[map, key, value]]];
+
 mapAdd[map_, key_, value_] :=
   Module[{},
     VerifyMap[map];
diff --git a/Tools/MathematicaMisc/Profile.m b/Tools/MathematicaMisc/Profile.m
new file mode 100644
index 0000000..8df42b8
--- /dev/null
+++ b/Tools/MathematicaMisc/Profile.m
@@ -0,0 +1,113 @@
+(* Copyright (C) 2010 Ian Hinder and Barry Wardell *)
+
+BeginPackage["Profile`"];
+
+Profile;
+ProfileTime;
+ProfileCount;
+ClearProfile;
+Timer;
+GetTimers;
+RemoveTimers;
+ThresholdTimers;
+PrintTimerTree;
+CoalesceTimers;
+
+Begin["`Private`"];
+
+ClearProfile[] :=
+  Module[{},
+    Clear[ProfileTime];
+    Clear[ProfileCount]];
+
+SetAttributes[Profile, HoldAll];
+
+Profile[name_, code_] :=
+  Module[{time, result, name2, result2, subTimers},
+    name2 = Evaluate[name];
+    {result2,subTimers} = Reap[{time, result} = AbsoluteTiming[ReleaseHold[code]]];
+    If[Head[ProfileTime[name2]] === ProfileTime, ProfileTime[name2] = 0.0];
+    If[Head[ProfileCount[name2]] === ProfileCount, ProfileCount[name2] = 0];
+    ProfileTime[name2] += time;
+    ProfileCount[name2] += 1;
+    Sow[Timer[name2,time,If[subTimers === {}, {}, subTimers[[1]]]]];
+    result];
+
+(* Profile[n_, x_] := x; *)
+
+(* The code below is a prototype for handling a tree of timer
+   information from Profile.  It needs to be tidied up and a good
+   interface defined.  Example of use:
+
+     {result, timers} = Reap@MyProfile["Total", <code>];
+
+     treeView@addUntimed@
+       Timer[timers[[1, 1, 1]], timers[[1, 1, 2]],
+         addTimers[timers[[1, 1, 3]], sameq, plus]]
+*)
+
+treeView[t : Timer[n_, v_, c_]] :=
+ Module[{},
+  If[c === {},
+   Row[{"  ", v, " ", n}],
+   OpenerView[{Row[{v, " ", n}],
+     Column[treeView /@ c]}, True]]]
+
+(* Take a list and, for all the members which are "sameq" as each \
+other, replace them with the application of "plus" to them *)
+add[l_, sameq_, plus_] :=
+ Map[plus, Gather[l, sameq]];
+
+addTimers[ts_List, sameq_, plus_] :=
+ add[Map[Timer[#[[1]], #[[2]], addTimers[#[[3]], sameq, plus]] &, ts],
+   sameq, plus];
+
+sameq[Timer[n1_, v1_, s1_], Timer[n2_, v2_, s2_]] :=
+ n1 === n2 && Length[s1] == Length[s2] &&
+  And @@ MapThread[sameq, {s1, s2}];
+
+plus[ts_List] :=
+ Timer[ts[[1, 1]], Plus @@ Map[#[[2]] &, ts],
+  MapThread[plus[{##}] &, Map[#[[3]] &, ts]]];
+
+addUntimed[Timer[n_, v_, s_List]] :=
+ If[Length[s] === 0, Timer[n, v, s],
+  Module[{total},
+   total = Plus @@ Map[#[[2]] &, s];
+   Timer[n, v,
+    Append[Map[addUntimed, s], Timer["untimed", v - total, {}]]]]];
+
+SetAttributes[GetTimers, HoldAll];
+GetTimers[expr_] :=
+  Module[
+    {result,timers},
+    {result,timers} = Reap[Profile["Total", expr]];
+    {result, timers[[1,1]]}];
+
+CoalesceTimers[t_Timer] :=
+  addTimers[{t}, sameq, plus][[1]];
+
+(* GetTimers[expr_] := *)
+(*   addUntimed[ *)
+(*     Timer[timers[[1, 1, 1]], timers[[1, 1, 2]], *)
+(*           addTimers[timers[[1, 1, 3]], sameq, plus]]]; *)
+
+RemoveTimers[timers_, remove_List] :=
+  timers //. Table[Timer[name, _, {sub___}] :> sub, {name, remove}];
+
+ThresholdTimers[t:Timer[name_, value_, sub_], threshold_:0.1] :=
+  thresholdTimers[t, threshold*value];
+
+thresholdTimers[t:Timer[name_, value_, sub_], absThreshold_] :=
+  DeleteCases[t, Timer[n_, v_, s_] /; v < absThreshold, Infinity] /. 
+  Timer[tt_,v_,{Timer["untimed",_,{}]}] :> Timer[tt,v,{}];  
+
+PrintTimerTree[Timer[name_, value_, sub_], indent_:0] :=
+  Module[
+    {},
+    Print[StringJoin@@ConstantArray["| ", indent],"o",PaddedForm[value,{3,2}], "  ", name];
+    Map[PrintTimerTree[#, indent+1]&, sub]];
+
+End[];
+
+EndPackage[];
diff --git a/Tools/MathematicaMisc/RunKranc.m b/Tools/MathematicaMisc/RunKranc.m
index bd6dec6..b554c9c 100644
--- a/Tools/MathematicaMisc/RunKranc.m
+++ b/Tools/MathematicaMisc/RunKranc.m
@@ -8,6 +8,7 @@ $Path = Join[$Path,
    krancDir <> "/Tools/External"}];
 Needs["Errors`"];
 Needs["KrancThorn`"];
+Needs["Profile`"];
 If[Environment["KRANCVERBOSE"] == "yes",
   SetDebugLevel[InfoFull]];
 
@@ -59,7 +60,23 @@ SetOptions["stdout", PageWidth -> Infinity];
 
 exception = Catch[Catch[
   Check[
-    Get[script];None,
+    Block[
+      {$RecursionLimit = Infinity},
+      (*{result,timers} =  GetTimers[ *) Get[script](*]*)];
+
+    (* Put[timers, "timer-output-1.m"]; *)
+
+    (* timers = CoalesceTimers[timers]; *)
+    (* Put[timers, "timer-output-2.m"]; *)
+
+    (* timers = ThresholdTimers[timers,0.1]; *)
+    (* Put[timers, "timer-output-3.m"]; *)
+
+    (* (\* Put[timers2, "timer-output-2.m"]; *\) *)
+
+    (* PrintTimerTree[timers]; *)
+
+    None,
     ThrowError["Messages were generated - aborted"]]], _];
 
 If[exception =!= None,