Merge remote-tracking branch 'origin/master' into hydro

25 files changed, 246 insertions, 1678 deletions

diff --git a/Auxiliary/Cactus/IOASCII.m b/Auxiliary/Cactus/IOASCII.m
new file mode 100644
index 0000000..118db8d
--- /dev/null
+++ b/Auxiliary/Cactus/IOASCII.m
@@ -0,0 +1,26 @@
+
+BeginPackage["IOASCII`"];
+
+ReadIOASCII::usage = "ReadIOASCII[file] reads the IOASCII file and parses it into list format";
+MGraph::usage = "MGraph[file] plots the IOASCII file";
+
+Begin["`Private`"];
+
+ReadIOASCII[file_] :=
+ Module[{data1, data2, data3},
+  data1 = Import[file, "Table"];
+  data2 = Select[SplitBy[data1, Length[#] == 0 &], #[[1]] != {} &];
+  data3 = Map[{First[#][[3]], Drop[#, 1]} &, data2]];
+
+MGraph[file_String] :=
+ Module[{data = ReadIOASCII[file]},
+  MGraph[data]];
+
+MGraph[data_List] :=
+ Module[{},
+  Manipulate[
+   ListLinePlot[data[[it, 2]], PlotLabel -> data[[it, 1]]], {{it,1,"it"}, 1, Length[data], 1}]];
+
+End[];
+
+EndPackage[];
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/interface.ccl b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/interface.ccl
deleted file mode 100644
index 8a9294b..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/interface.ccl
+++ /dev/null
@@ -1,8 +0,0 @@
-# file produced by user shusa, 31/3/2004
-
-# $Id$
-
-implements: Ceiling
-
-inherits: Grid 
-
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/param.ccl b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/param.ccl
deleted file mode 100644
index a6dc073..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/param.ccl
+++ /dev/null
@@ -1,32 +0,0 @@
-# file produced by user shusa, 31/3/2004
-# Produced with Mathematica Version 5.0 for Linux (June 9, 2003)
-
-# Mathematica script written by Ian Hinder and Sascha Husa
-
-# $Id$
-
-private:
-BOOLEAN check_active "whether to check ceiling values"
-{
-} "false"
-
-REAL ceiling_value "what value we use for the ceiling"
-{
-*:* :: "with negative values no cutoff is set"
-} -1
-
-KEYWORD type "what type of checking to apply"
-{
-  "absolute" :: "check for absolute value of GF"
-  "differential" :: "check for difference between min & max"
-} "absolute"
-
-STRING vars "List of evolved grid functions that should have dissipation added" STEERABLE=always
-{
-  .* :: "must be a valid list of grid functions"
-} ""
-
-BOOLEAN verbose "produce log output" STEERABLE=always
-{
-} "no"
-
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/schedule.ccl b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/schedule.ccl
deleted file mode 100644
index b9d167c..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/schedule.ccl
+++ /dev/null
@@ -1,21 +0,0 @@
-# file produced by user shusa, 31/3/2004
-# Produced with Mathematica Version 5.0 for Linux (June 9, 2003)
-
-# Mathematica script written by Ian Hinder and Sascha Husa
-
-# $Id$
-
-
-schedule Ceiling_Startup at startup
-{
-LANG: C
-} "ceiling startup message" 
-
-if (check_active)
-{
-  schedule check_ceiling at PostStep
-  {
-    LANG: C
-    
-  } "check ceiling"
-}
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Ceiling.F90 b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Ceiling.F90
deleted file mode 100644
index 9e1d2c1..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Ceiling.F90
+++ /dev/null
@@ -1,68 +0,0 @@
-! file written by s. husa, 5/6/2004
-
-! $Id$
-
-#include "cctk.h"
-
-subroutine apply_check_abs(var, ni, nj, nk, ceiling_value)
-
-implicit none
-
-CCTK_INT,                          intent(in) :: ni, nj, nk
-CCTK_REAL, dimension (ni, nj, nk), intent(in) :: var(ni, nj, nk)
-CCTK_REAL,                         intent(in) :: ceiling_value
-
-CCTK_REAL       :: criterion
-CCTK_REAL, save :: initial_value
-
-CCTK_INT,  save :: counter
-
-counter   = counter + 1
-criterion = maxval(abs(var) + epsilon(1.0d0))
-
-if (counter == 1) then
-  initial_value = criterion
-  write (*,*) "<<<<<< using ceiling initial value", initial_value
-else
-  criterion = criterion / initial_value
-  if ((ceiling_value > 0).AND.(criterion > ceiling_value)) then
-          
-               call CCTK_INFO("Ceiling thorn terminates evolution")
-               call CCTK_TerminateNext(var)
-  endif
-endif
-
-end subroutine apply_check_abs
-
-
-subroutine apply_check_diff(var, ni, nj, nk, ceiling_value)
-
-implicit none
-
-CCTK_INT,                          intent(in) :: ni, nj, nk
-CCTK_REAL, dimension (ni, nj, nk), intent(in) :: var(ni, nj, nk)
-CCTK_REAL,                         intent(in) :: ceiling_value
-
-CCTK_REAL       :: criterion
-CCTK_REAL, save :: initial_value
-
-CCTK_INT,  save :: counter
-
-counter   = counter + 1
-
-criterion = maxval(var) - minval(var)
-
-if (counter == 1) then
-  initial_value = criterion
-  write (*,*) "<<<<<< using ceiling initial value", initial_value
-else
-  criterion = criterion / initial_value
-  if ((ceiling_value > 0).AND.(criterion > ceiling_value)) then
-
-               call CCTK_INFO("Ceiling thorn terminates evolution.")
-               call CCTK_TerminateNext (var)
-  endif
-endif
-
-end subroutine apply_check_diff
-
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Startup.c b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Startup.c
deleted file mode 100644
index b42b2a1..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Startup.c
+++ /dev/null
@@ -1,15 +0,0 @@
-/*  file produced by user shusa, 31/3/2004 */
-/*  Produced with Mathematica Version 5.0 for Linux (June 9, 2003) */
-
-/*  Mathematica script written by Ian Hinder and Sascha Husa */
-
-/*  $Id$ */
-
-#include "cctk.h"
-
-int Ceiling_Startup(void)
-{
-  const char * banner = "Ceiling: abort when solution grows through the roof";
-  CCTK_RegisterBanner(banner);
-  return 0;
-}
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/make.code.defn b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/make.code.defn
deleted file mode 100644
index dd1377a..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/make.code.defn
+++ /dev/null
@@ -1,8 +0,0 @@
-# file produced by user shusa, 31/3/2004
-# Produced with Mathematica Version 5.0 for Linux (June 9, 2003)
-
-# Mathematica script written by Ian Hinder and Sascha Husa
-
-# $Id$
-
-SRCS = Startup.c Ceiling.F90 selectGFs.c
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/selectGFs.c b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/selectGFs.c
deleted file mode 100644
index 86e1196..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/selectGFs.c
+++ /dev/null
@@ -1,94 +0,0 @@
-/* $Header$ */
-
-/* this code is based on Erik Schnetter's dissipation thorn */
-
-#include <assert.h>
-#include <stdlib.h>
-
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-
-void CCTK_FCALL
-CCTK_FNAME(apply_check_abs) (CCTK_REAL const * const var,
-                             int       const * const ni,
-                             int       const * const nj,
-                             int       const * const nk,
-                             CCTK_REAL const * const ceiling_value);
-
-
-void CCTK_FCALL
-CCTK_FNAME(apply_check_diff) (CCTK_REAL const * const var,
-                              int       const * const ni,
-                              int       const * const nj,
-                              int       const * const nk,
-                              CCTK_REAL const * const ceiling_value);
-
-static void
-call_apply_check (int const varindex, char const * const optstring, void * const arg);
-
-void
-check_ceiling (CCTK_ARGUMENTS)
-{
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
-
-  CCTK_TraverseString (vars, call_apply_check, cctkGH, CCTK_GROUP_OR_VAR);
-}
-
-
-void
-call_apply_check (int const varindex, char const * const optstring, void * const arg)
-{
-  cGH const * const cctkGH = (cGH const *) arg;
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
-
-  int vargroup;
-  cGroup vardata;
-
-  CCTK_REAL const * varptr;
-  int ierr /* , terminate */ ;
-
-  assert (varindex >= 0);
-
-  if (verbose) {
-    char * const fullvarname = CCTK_FullName (varindex);
-    assert (fullvarname);
-    CCTK_VInfo (CCTK_THORNSTRING,
-                "Applying ceiling check to \"%s\" ",
-                fullvarname);
-    free (fullvarname);
-  }
-
-  vargroup = CCTK_GroupIndexFromVarI (varindex);
-  assert (vargroup >= 0);
-
-  ierr = CCTK_GroupData (vargroup, &vardata);
-  assert (!ierr);
-
-  assert (vardata.grouptype == CCTK_GF);
-  assert (vardata.vartype == CCTK_VARIABLE_REAL);
-  assert (vardata.dim == cctk_dim);
-
-  varptr = CCTK_VarDataPtrI (cctkGH, 0, varindex);
-  assert (varptr);
-
-  if (CCTK_Equals (type, "absolute"))
-  {
-   CCTK_FNAME(apply_check_abs)
-      (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2], &ceiling_value);
-  }
-  else if (CCTK_Equals (type, "differential"))
-  {
-   CCTK_FNAME(apply_check_diff)
-      (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2], &ceiling_value);
-  }
-  else
-  {
-   CCTK_INFO("keyword ceiling::type only allows values 'absolute' and 'differential'");
-  }
-
-  /* if (terminate > 0) {CCTK_TerminateNext (cctkGH);}  */
-}
-
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/param.ccl b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/param.ccl
index 430c2f2..8e410f7 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/param.ccl
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/param.ccl
@@ -4,38 +4,46 @@
 
 # $Header$ 
 
-private:    
-
-KEYWORD FDscheme "Type of finite differencing to use"
-{                                                     
-  "2nd order centered macro"      :: "centered 2nd order implemented with macros"          
-  "4th order centered macro"      :: "centered 4th order implemented with macros"
-} "2nd order centered macro"
-
 restricted:
 
-CCTK_INT  stencil_width "stencil width used near boundary"
+CCTK_INT  stencil_width "stencil width used near boundary DEPRECATED"
 {
   -1:* :: "outgoing characteristic speed > 0, default of -1 is intentionally invalid"
 } -1
 
-CCTK_INT  stencil_width_x "stencil width used near boundary"
+CCTK_INT  stencil_width_x "stencil width used near boundary DEPRECATED"
 {
   -1:* :: "outgoing characteristic speed > 0, default of -1 is intentionally invalid"
 } -1
 
-CCTK_INT  stencil_width_y "stencil width used near boundary"
+CCTK_INT  stencil_width_y "stencil width used near boundary DEPRECATED"
 {
   -1:* :: "outgoing characteristic speed > 0, default of -1 is intentionally invalid"
 } -1
 
-CCTK_INT  stencil_width_z "stencil width used near boundary"
+CCTK_INT  stencil_width_z "stencil width used near boundary DEPRECATED"
 {
   -1:* :: "outgoing characteristic speed > 0, default of -1 is intentionally invalid"
 } -1
 
-CCTK_INT  boundary_width "width of boundary (fix later to use Cactus boundary calls)"
+CCTK_INT  boundary_width "width of boundary (fix later to use Cactus boundary calls) DEPRECATED"
 {
   -1:* :: "Any integer"
 } 1
 
+restricted:
+CCTK_STRING jacobian_group "Name of group containing Jacobian" STEERABLE=RECOVER
+{
+  "" :: "String of the form <implementation>::<groupname>"
+} ""
+
+restricted:
+CCTK_STRING jacobian_derivative_group "Name of group containing Jacobian derivative" STEERABLE=RECOVER
+{
+  "" :: "String of the form <implementation>::<groupname>"
+} ""
+
+CCTK_INT  jacobian_identity_map "Map number on which the Jacobian should not be applied"
+{
+  -1:* :: "Any integer"
+} -1
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/schedule.ccl b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/schedule.ccl
index a013b00..05ca0e1 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/schedule.ccl
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/schedule.ccl
@@ -3,16 +3,3 @@
 # author: S. Husa
 
 # $Header$ 
-
-
-
-schedule GenericFD_Startup at STARTUP          
-{                            
-  LANG: C                    
-}   "Register Banner"        
-
-# schedule GenericFD_ParamCheck at ParamCheck
-# {
-#   LANG: Fortran
-# }   "check stencil width parameters for consistency"
-
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
index 4e1e774..8c06940 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
@@ -59,18 +59,28 @@ int sgn(CCTK_REAL x)
     return 0;
 }
 
-int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH)
+void GenericFD_GetBoundaryWidths(cGH const * restrict const cctkGH, int nboundaryzones[6])
 {
   int is_internal[6];
   int is_staggered[6];
-  int nboundaryzones[6];
   int shiftout[6];
   int ierr = -1;
 
   if (CCTK_IsFunctionAliased ("MultiPatch_GetBoundarySpecification")) {
     int const map = MultiPatch_GetMap (cctkGH);
+    /* This doesn't make sense in level mode */
     if (map < 0)
-      CCTK_WARN(0, "Could not determine current map");
+    {
+      static int have_warned = 0;
+      if (!have_warned)
+      {
+        CCTK_WARN(1, "GenericFD_GetBoundaryWidths: Could not determine current map (can be caused by calling in LEVEL mode)");
+        have_warned = 1;
+      }
+      for (int i = 0; i < 6; i++)
+        nboundaryzones[i] = 0;
+      return;
+    }
     ierr = MultiPatch_GetBoundarySpecification
       (map, 6, nboundaryzones, is_internal, is_staggered, shiftout);
     if (ierr != 0)
@@ -83,6 +93,12 @@ int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH)
   } else {
     CCTK_WARN(0, "Could not obtain boundary specification");
   }
+}
+
+int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH)
+{
+  int nboundaryzones[6];
+  GenericFD_GetBoundaryWidths(cctkGH, nboundaryzones);
 
   int bw = nboundaryzones[0];
 
@@ -93,20 +109,6 @@ int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH)
   return bw;
 }
 
-/* int GenericFD_BoundaryWidthTable(cGH const * restrict const cctkGH) */
-/* { */
-/*   int nboundaryzones[6]; */
-/*   GenericFD_GetBoundaryWidth(cctkGH, nboundaryzones); */
-
-/*   int table = Util_TableCreate(0); */
-/*   if (table < 0) CCTK_WARN(0, "Could not create table"); */
-
-/*   if (Util_TableSetIntArray(table, 6, nboundaryzones, "BOUNDARY_WIDTH") < 0) */
-/*     CCTK_WARN(0, "Could not set table"); */
-/*   return table; */
-/* } */
-
-
 /* Return the array indices in imin and imax for looping over the
    interior of the grid. imin is the index of the first grid point.
    imax is the index of the last grid point plus 1.  So a loop over
@@ -219,7 +221,7 @@ void GenericFD_GetBoundaryInfo(cGH const * restrict const cctkGH,
 	imin[dir] = npoints;
 	break;
       case 1: /* Upper boundary */
-	imax[dir] = cctk_lssh[CCTK_LSSH_IDX(0,dir)] - npoints;
+	imax[dir] = CCTK_LSSH(0,dir) - npoints;
 	break;
       default:
 	CCTK_WARN(0, "internal error");
@@ -239,7 +241,7 @@ void GenericFD_LoopOverEverything(cGH const * restrict const cctkGH, Kranc_Calcu
   CCTK_REAL  tangentA[] = {0,0,0};
   CCTK_REAL  tangentB[] = {0,0,0};
   int   bmin[] = {0,0,0};
-  int   bmax[] = {cctk_lssh[CCTK_LSSH_IDX(0,0)],cctk_lssh[CCTK_LSSH_IDX(0,1)],cctk_lssh[CCTK_LSSH_IDX(0,2)]};
+  int   bmax[] = {CCTK_LSSH(0,0),CCTK_LSSH(0,1),CCTK_LSSH(0,2)};
 
   calc(cctkGH, dir, face, normal, tangentA, tangentB, bmin, bmax, 0, NULL);
   return;
@@ -301,7 +303,7 @@ void GenericFD_LoopOverBoundary(cGH const * restrict const cctkGH, Kranc_Calcula
             break;
           case +1:
             bmin[d] = imax[d];
-            bmax[d] = cctk_lssh[CCTK_LSSH_IDX(0,d)];
+            bmax[d] = CCTK_LSSH(0,d);
             have_bnd = 1;
             all_physbnd = all_physbnd && is_physbnd[2*d+1];
             break;
@@ -394,7 +396,7 @@ void GenericFD_LoopOverBoundaryWithGhosts(cGH const * restrict const cctkGH, Kra
             break;
           case +1:
             bmin[d] = imax[d];
-            bmax[d] = cctk_lssh[CCTK_LSSH_IDX(0,d)];
+            bmax[d] = CCTK_LSSH(0,d);
             have_bnd = 1;
             have_physbnd = have_physbnd || is_physbnd[2*d+1];
             break;
@@ -474,7 +476,7 @@ void GenericFD_PenaltyPrim2Char(cGH const * restrict const cctkGH, int const dir
   CCTK_REAL  tangentA[] = {0,0,0};
   CCTK_REAL  tangentB[] = {0,0,0};
   int   bmin[] = {0,0,0};
-  int   bmax[] = {cctk_lssh[CCTK_LSSH_IDX(0,0)],cctk_lssh[CCTK_LSSH_IDX(0,1)],cctk_lssh[CCTK_LSSH_IDX(0,2)]};
+  int   bmax[] = {cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]};
   CCTK_REAL  **all_vars;
   int        i = 0;
 
@@ -499,3 +501,100 @@ void GenericFD_PenaltyPrim2Char(cGH const * restrict const cctkGH, int const dir
   
   return;
 }
+
+void GenericFD_AssertGroupStorage(cGH const * restrict const cctkGH, const char *calc,
+                                  int ngroups, const char *group_names[ngroups])
+{
+  for (int i = 0; i < ngroups; i++)
+  {
+    int result = CCTK_QueryGroupStorage(cctkGH, group_names[i]);
+    if (result == 0)
+    {
+      CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                 "Error in %s: Group \"%s\" does not have storage", calc, group_names[i]);
+    }
+    else if (result < 0)
+    {
+      CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                 "Error in %s: Invalid group name \"%s\"", calc, group_names[i]);
+    }
+  }
+}
+
+/* Return a list of pointers to the members of a named group */
+void GenericFD_GroupDataPointers(cGH const * restrict const cctkGH, const char *group_name,
+                                 int nvars, CCTK_REAL const *restrict *ptrs)
+{
+  int group_index, status;
+  cGroup  group_info;
+
+  group_index = CCTK_GroupIndex(group_name);
+  if (group_index < 0)
+    CCTK_VWarn(CCTK_WARN_ABORT,   __LINE__, __FILE__, CCTK_THORNSTRING,
+               "Error return %d trying to get group index for group \'%s\'",
+               group_index,
+               group_name);
+
+  status = CCTK_GroupData(group_index, &group_info);
+  if (status < 0)
+    CCTK_VWarn(CCTK_WARN_ABORT,   __LINE__, __FILE__, CCTK_THORNSTRING,
+               "Error return %d trying to get info for group \'%s\'",
+               status,
+               group_name);
+
+  if (group_info.numvars != nvars)
+  {
+    CCTK_VWarn(CCTK_WARN_ABORT,   __LINE__, __FILE__, CCTK_THORNSTRING,
+               "Group \'%s\' has %d variables but %d were expected",
+               group_name, group_info.numvars, nvars);
+  }
+
+  int v1 = CCTK_FirstVarIndex(group_name);
+
+  for (int v = 0; v < nvars; v++)
+  {
+    ptrs[v] = (CCTK_REAL const *) CCTK_VarDataPtrI(cctkGH, 0 /* timelevel */, v1+v);
+  }
+}
+
+void GenericFD_EnsureStencilFits(cGH const * restrict const cctkGH, const char *calc, int ni, int nj, int nk)
+{
+  DECLARE_CCTK_ARGUMENTS
+
+  int bws[6];
+  GenericFD_GetBoundaryWidths(cctkGH, bws);
+
+  int ns[] = {ni, nj, nk};
+  const char *dirs[] = {"x", "y", "z"};
+  const char *faces[] = {"lower", "upper"};
+  int abort = 0;
+
+  for (int dir = 0; dir < 3; dir++)
+  {
+    for (int face = 0; face < 2; face++)
+    {
+      int bw = bws[2*dir+face];
+      if (bw < ns[dir])
+      {
+        CCTK_VInfo(CCTK_THORNSTRING,
+                   "The stencil for %s requires %d points, but the %s %s boundary has only %d points.",
+                   calc, ns[dir], faces[face], dirs[dir], bw);
+        abort = 1;
+      }
+    }
+    int gz = cctk_nghostzones[dir];
+    if (gz < ns[dir])
+    {
+      CCTK_VInfo(CCTK_THORNSTRING,
+                 "The stencil for %s requires %d points, but there are only %d ghost zones in the %s direction.",
+                 calc, ns[dir], gz, dirs[dir]);
+      abort = 1;
+    }
+  }
+
+  if (abort)
+  {
+    CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "Insufficient ghost or boundary points for %s", calc);
+  }
+}
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
index 041347d..401b4e0 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
@@ -27,655 +27,33 @@
     along with Kranc; if not, write to the Free Software
     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */
-                         
-#ifndef NOPRECOMPUTE
-#define PRECOMPUTE
-#endif
-
-#if defined(FD_C0) || defined(FD_C2) || defined(FD_C4) || defined(FD_C2C4)
-#define FD_SET_BY_USER
-#endif
-#ifndef FD_SET_BY_USER
-#define FD_C2
-#endif	
 
+#include "cctk.h"
 
-#if defined(FD_C0)
-#define FD_METHOD_DESC  "FD method: replace derivatives by zero"
+#ifdef __cplusplus
+extern "C" {
 #endif
 
-#if defined(FD_C2)
-#define FD_METHOD_DESC  "FD method: second order centered finite differences"
-#endif
-
-#if defined(FD_C4)
-#define FD_METHOD_DESC  "FD method: fourth order centered finite differences"
-#endif
-
-#if defined(FD_C2C4)
-#define FD_METHOD_DESC  "FD method: weighted 2nd/4th order centered finite differences"
-#endif
-
-
-
-/* utility functions  */
-#if defined(KRANC_C)
-#define string(d,f) d ## f
-#else
-#define string(d,f) d/**/f
+#ifdef __cplusplus
+#  ifdef CCTK_CXX_RESTRICT
+#    define restrict CCTK_CXX_RESTRICT
+#  endif
 #endif
 
-#if defined(KRANC_C)
-#define IMAX(int1, int2) ((int1) > (int2) ? (int1) : (int2))
-#endif
-
-
 #include "MathematicaCompat.h"
 
-/* finite differencing macros */
-
-/*                                   */
-/* add method argument to shorthands */
-/*                                   */
-
-/* third derivatives */
-#define D111x(gf)  string(D111,gf)
-#define D211x(gf)  string(D211,gf)
-#define D311x(gf)  string(D311,gf)
-#define D221x(gf)  string(D221,gf)
-#define D321x(gf)  string(D321,gf)
-#define D331x(gf)  string(D331,gf)
-#define D222x(gf)  string(D222,gf)
-#define D322x(gf)  string(D322,gf)
-#define D332x(gf)  string(D332,gf)
-#define D333x(gf)  string(D333,gf)
-
-/* second derivatives */
-#define D11x(gf)  string(D11,gf)
-#define D22x(gf)  string(D22,gf)
-#define D33x(gf)  string(D33,gf)
-#define D21x(gf)  string(D21,gf)
-#define D32x(gf)  string(D32,gf)
-#define D31x(gf)  string(D31,gf)
-
-/* first derivatives                */
-#define D1x(gf)   string(D1,gf)
-#define D2x(gf)   string(D2,gf)
-#define D3x(gf)   string(D3,gf)
-
-#ifdef PRECOMPUTE
-
-/* third derivatives */
-#define D111(gf)  string(D111,gf)
-#define D211(gf)  string(D211,gf)
-#define D311(gf)  string(D311,gf)
-#define D221(gf)  string(D221,gf)
-#define D321(gf)  string(D321,gf)
-#define D331(gf)  string(D331,gf)
-#define D222(gf)  string(D222,gf)
-#define D322(gf)  string(D322,gf)
-#define D332(gf)  string(D332,gf)
-#define D333(gf)  string(D333,gf)
-
-/* second derivatives */
-#define D11(gf,i,j,k)  string(D11,gf)
-#define D22(gf,i,j,k)  string(D22,gf)
-#define D33(gf,i,j,k)  string(D33,gf)
-#define D21(gf,i,j,k)  string(D21,gf)
-#define D32(gf,i,j,k)  string(D32,gf)
-#define D31(gf,i,j,k)  string(D31,gf) 
-
-/* first derivatives                */
-#define D1(gf,i,j,k)   string(D1,gf)
-#define D2(gf,i,j,k)   string(D2,gf) 
-#define D3(gf,i,j,k)   string(D3,gf) 
-
-#else
-
-/* third derivatives */
-#define D111(gf)  D111gf(gf,i,j,k)
-#define D211(gf)  D211gf(gf,i,j,k)
-#define D311(gf)  D311gf(gf,i,j,k)
-#define D221(gf)  D221gf(gf,i,j,k)
-#define D321(gf)  D321gf(gf,i,j,k)
-#define D331(gf)  D331gf(gf,i,j,k)
-#define D222(gf)  D222gf(gf,i,j,k)
-#define D322(gf)  D322gf(gf,i,j,k)
-#define D332(gf)  D332gf(gf,i,j,k)
-#define D333(gf)  D333gf(gf,i,j,k)
-
-/* second derivatives */
-#define D11(gf,i,j,k)  D11gf(gf,i,j,k)
-#define D22(gf,i,j,k)  D22gf(gf,i,j,k)
-#define D33(gf,i,j,k)  D33gf(gf,i,j,k)
-#define D21(gf,i,j,k)  D21gf(gf,i,j,k)
-#define D32(gf,i,j,k)  D32gf(gf,i,j,k)
-#define D31(gf,i,j,k)  D31gf(gf,i,j,k)
-
-/* first derivatives                */
-#define D1(gf,i,j,k)   D1gf(gf, i,j,k)
-#define D2(gf,i,j,k)   D2gf(gf, i,j,k)
-#define D3(gf,i,j,k)   D3gf(gf, i,j,k)
-
-#endif
-
-
-#ifdef FD_C0
-/* third derivatives */
-#define D111gf(gf,i,j,k)  D111_c0(gf,i,j,k)
-#define D211gf(gf,i,j,k)  D211_c0(gf,i,j,k)
-#define D311gf(gf,i,j,k)  D311_c0(gf,i,j,k)
-#define D221gf(gf,i,j,k)  D221_c0(gf,i,j,k)
-#define D321gf(gf,i,j,k)  D321_c0(gf,i,j,k)
-#define D331gf(gf,i,j,k)  D331_c0(gf,i,j,k)
-#define D222gf(gf,i,j,k)  D222_c0(gf,i,j,k)
-#define D322gf(gf,i,j,k)  D322_c0(gf,i,j,k)
-#define D332gf(gf,i,j,k)  D332_c0(gf,i,j,k)
-#define D333gf(gf,i,j,k)  D333_c0(gf,i,j,k)
-
-/* second derivatives */
-#define D11gf(gf,i,j,k)  D11_c0(gf,i,j,k)
-#define D22gf(gf,i,j,k)  D22_c0(gf,i,j,k)
-#define D33gf(gf,i,j,k)  D33_c0(gf,i,j,k)
-#define D21gf(gf,i,j,k)  D21_c0(gf,i,j,k)
-#define D32gf(gf,i,j,k)  D32_c0(gf,i,j,k)
-#define D31gf(gf,i,j,k)  D31_c0(gf,i,j,k)
-
-/* first derivatives                */
-#define D1gf(gf,i,j,k)   D1_c0(gf, i,j,k)
-#define D2gf(gf,i,j,k)   D2_c0(gf, i,j,k)
-#define D3gf(gf,i,j,k)   D3_c0(gf, i,j,k)
-#endif
-
-
-
-#ifdef FD_C2
-/* third derivatives */
-#define D111gf(gf,i,j,k)  D111_c2(gf,i,j,k)
-#define D211gf(gf,i,j,k)  D211_c2(gf,i,j,k)
-#define D311gf(gf,i,j,k)  D311_c2(gf,i,j,k)
-#define D221gf(gf,i,j,k)  D221_c2(gf,i,j,k)
-#define D321gf(gf,i,j,k)  D321_c2(gf,i,j,k)
-#define D331gf(gf,i,j,k)  D331_c2(gf,i,j,k)
-#define D222gf(gf,i,j,k)  D222_c2(gf,i,j,k)
-#define D322gf(gf,i,j,k)  D322_c2(gf,i,j,k)
-#define D332gf(gf,i,j,k)  D332_c2(gf,i,j,k)
-#define D333gf(gf,i,j,k)  D333_c2(gf,i,j,k)
-
-/* second derivatives */
-#define D11gf(gf,i,j,k)  D11_c2(gf,i,j,k)
-#define D22gf(gf,i,j,k)  D22_c2(gf,i,j,k)
-#define D33gf(gf,i,j,k)  D33_c2(gf,i,j,k)
-#define D21gf(gf,i,j,k)  D21_c2(gf,i,j,k)
-#define D32gf(gf,i,j,k)  D32_c2(gf,i,j,k)
-#define D31gf(gf,i,j,k)  D31_c2(gf,i,j,k)
-
-/* first derivatives                */
-#define D1gf(gf,i,j,k)   D1_c2(gf, i,j,k)
-#define D2gf(gf,i,j,k)   D2_c2(gf, i,j,k)
-#define D3gf(gf,i,j,k)   D3_c2(gf, i,j,k)
-#endif
-
-
-
-#ifdef FD_C4
-/* third derivatives */
-#define D111gf(gf,i,j,k)  D111_c4(gf,i,j,k)
-#define D211gf(gf,i,j,k)  D211_c4(gf,i,j,k)
-#define D311gf(gf,i,j,k)  D311_c4(gf,i,j,k)
-#define D221gf(gf,i,j,k)  D221_c4(gf,i,j,k)
-#define D321gf(gf,i,j,k)  D321_c4(gf,i,j,k)
-#define D331gf(gf,i,j,k)  D331_c4(gf,i,j,k)
-#define D222gf(gf,i,j,k)  D222_c4(gf,i,j,k)
-#define D322gf(gf,i,j,k)  D322_c4(gf,i,j,k)
-#define D332gf(gf,i,j,k)  D332_c4(gf,i,j,k)
-#define D333gf(gf,i,j,k)  D333_c4(gf,i,j,k)
-
-/* second derivatives */
-#define D11gf(gf,i,j,k)  D11_c4(gf,i,j,k) 
-#define D22gf(gf,i,j,k)  D22_c4(gf,i,j,k)
-#define D33gf(gf,i,j,k)  D33_c4(gf,i,j,k)
-#define D21gf(gf,i,j,k)  D21_c4(gf,i,j,k) 
-#define D32gf(gf,i,j,k)  D32_c4(gf,i,j,k) 
-#define D31gf(gf,i,j,k)  D31_c4(gf,i,j,k)
-
-/* first derivatives                */
-#define D1gf(gf,i,j,k)   D1_c4(gf, i,j,k) 
-#define D2gf(gf,i,j,k)   D2_c4(gf, i,j,k)
-#define D3gf(gf,i,j,k)   D3_c4(gf, i,j,k)
-#endif
-
-
-#ifdef FD_C2C4
-/* third derivatives */
-#define D111gf(gf,i,j,k)  D111_c2c4(gf,i,j,k)
-#define D211gf(gf,i,j,k)  D211_c2c4(gf,i,j,k)
-#define D311gf(gf,i,j,k)  D311_c2c4(gf,i,j,k)
-#define D221gf(gf,i,j,k)  D221_c2c4(gf,i,j,k)
-#define D321gf(gf,i,j,k)  D321_c2c4(gf,i,j,k)
-#define D331gf(gf,i,j,k)  D331_c2c4(gf,i,j,k)
-#define D222gf(gf,i,j,k)  D222_c2c4(gf,i,j,k)
-#define D322gf(gf,i,j,k)  D322_c2c4(gf,i,j,k)
-#define D332gf(gf,i,j,k)  D332_c2c4(gf,i,j,k)
-#define D333gf(gf,i,j,k)  D333_c2c4(gf,i,j,k)
-
-/* second derivatives */
-#define D11gf(gf,i,j,k)  D11_c2c4(gf,i,j,k)
-#define D22gf(gf,i,j,k)  D22_c2c4(gf,i,j,k)
-#define D33gf(gf,i,j,k)  D33_c2c4(gf,i,j,k)
-#define D21gf(gf,i,j,k)  D21_c2c4(gf,i,j,k)
-#define D32gf(gf,i,j,k)  D32_c2c4(gf,i,j,k)
-#define D31gf(gf,i,j,k)  D31_c2c4(gf,i,j,k)
-
-/* first derivatives                */
-#define D1gf(gf,i,j,k)   D1_c2c4(gf, i,j,k)
-#define D2gf(gf,i,j,k)   D2_c2c4(gf, i,j,k)
-#define D3gf(gf,i,j,k)   D3_c2c4(gf, i,j,k)
-#endif
-
-
-
-
-/*****************************************************/
-/*                                                   */
-/*             METHODS                               */
-/*                                                   */
-/*****************************************************/
-
-/*  c0                              */
-
-/*  set all derivatives = 0         */
-/*  for debugging and benchmarking  */
-
-/* third derivatives */
-
-#define D111_c0(gf,i,j,k)   0.
-#define D211_c0(gf,i,j,k)   0.
-#define D311_c0(gf,i,j,k)   0.
-#define D221_c0(gf,i,j,k)   0.
-#define D321_c0(gf,i,j,k)   0.
-#define D331_c0(gf,i,j,k)   0.
-#define D222_c0(gf,i,j,k)   0.
-#define D322_c0(gf,i,j,k)   0.
-#define D332_c0(gf,i,j,k)   0.
-#define D333_c0(gf,i,j,k)   0.
-
-/*  second derivatives              */
-
-#define D11_c0(gf,i,j,k)   0.
-#define D22_c0(gf,i,j,k)   0.
-#define D33_c0(gf,i,j,k)   0.
-#define D21_c0(gf,i,j,k)   0.
-#define D32_c0(gf,i,j,k)   0.
-#define D31_c0(gf,i,j,k)   0.
-
-/* first derivatives                */
-
-#define D1_c0(gf,i,j,k)    0.
-#define D2_c0(gf,i,j,k)    0.
-#define D3_c0(gf,i,j,k)    0.
-
-
-
-#ifndef KRANC_C
-
-/*  c2                  */
-/*                      */
-/*  2nd order centered  */
-/*                      */
-
-/* third derivatives, centered, 2nd order */
-
-#define D111_c2(gf,i,j,k)   ((- gf(i+2,j,k) + 2*gf(i+1,j,k) - 2*gf(i-1,j,k) + gf(i-2,j,k)) * dxi*dxi*dxi * (1.0/2.0))
-#define D211_c2(gf,i,j,k)   ((gf(i+1,j+1,k) - 2*gf(i,j+1,k) + gf(i-1,j+1,k) - gf(i+1,j-1,k) + 2*gf(i,j-1,k) - gf(i-1,j-1,k)) * dxi*dxi*dyi * (1.0/2.0))
-#define D311_c2(gf,i,j,k)   ((gf(i+1,j,k+1) - 2*gf(i,j,k+1) + gf(i-1,j,k+1) - gf(i+1,j,k-1) + 2*gf(i,j,k-1) - gf(i-1,j,k-1)) * dxi*dxi*dzi * (1.0/2.0))
-#define D221_c2(gf,i,j,k)   ((gf(i+1,j+1,k) - 2*gf(i+1,j,k) + gf(i+1,j-1,k) - gf(i-1,j+1,k) + 2*gf(i-1,j,k) - gf(i-1,j-1,k)) * dxi*dyi*dyi * (1.0/2.0))
-#define D321_c2(gf,i,j,k)   ((gf(i+1,j+1,k+1) - gf(i-1,j+1,k+1) - gf(i+1,j-1,k+1) + gf(i-1,j-1,k+1) - gf(i+1,j+1,k-1) + gf(i-1,j+1,k-1) + gf(i+1,j-1,k-1) - gf(i-1,j-1,k-1)) * dxi*dyi*dzi * (1.0/8.0))
-#define D331_c2(gf,i,j,k)   ((gf(i+1,j,k+1) - 2*gf(i+1,j,k) + gf(i+1,j,k-1) - gf(i-1,j,k+1) + 2*gf(i-1,j,k) - gf(i-1,j,k-1)) * dxi*dzi*dzi * (1.0/2.0))
-#define D222_c2(gf,i,j,k)   ((- gf(i,j+2,k) + 2*gf(i,j+1,k) - 2*gf(i,j-1,k) + gf(i,j-2,k)) * dyi*dyi*dyi * (1.0/2.0))
-#define D322_c2(gf,i,j,k)   ((gf(i,j+1,k+1) - 2*gf(i,j,k+1) + gf(i,j-1,k+1) - gf(i,j+1,k-1) + 2*gf(i,j,k-1) - gf(i,j-1,k-1)) * dyi*dyi*dzi * (1.0/2.0))
-#define D332_c2(gf,i,j,k)   ((gf(i,j+1,k+1) - 2*gf(i,j+1,k) + gf(i,j+1,k-1) - gf(i,j-1,k+1) + 2*gf(i,j-1,k) - gf(i,j-1,k-1)) * dyi*dzi*dzi * (1.0/2.0))
-#define D333_c2(gf,i,j,k)   ((- gf(i,j,k+2) + 2*gf(i,j,k+1) - 2*gf(i,j,k-1) + gf(i,j,k-2)) * dzi*dzi*dzi * (1.0/2.0))
-
-/* second derivatives, centered, 2nd order */
-
-#define D11_c2(gf,i,j,k)                         \
-	 ((   gf(i+1,j,k) \
-	 - 2.*gf(i,  j,k) \
-	 +    gf(i-1,j,k)) * dxi * dxi )
-
-#define D22_c2(gf,i,j,k)                         \
-	 ((   gf(i,j+1,k) \
-	 - 2.*gf(i,j,  k) \
-	 +    gf(i,j-1,k)) * dyi * dyi )
-
-#define D33_c2(gf,i,j,k)                         \
-	 ((   gf(i,j,k+1) \
-	 - 2.*gf(i,j,k  ) \
-	 +    gf(i,j,k-1)) * dzi * dzi )
-
-#define D21_c2(gf,i,j,k)                        \
-	 ((gf(i+1,j+1,k) \
-	 + gf(i-1,j-1,k) \
-	 - gf(i+1,j-1,k) \
-	 - gf(i-1,j+1,k)) * hdxi * hdyi )
-
-#define D32_c2(gf,i,j,k)                        \
-	 ((gf(i,j+1,k+1) \
-	 + gf(i,j-1,k-1) \
-	 - gf(i,j+1,k-1) \
-	 - gf(i,j-1,k+1)) * hdyi * hdzi )
-
-#define D31_c2(gf,i,j,k)                        \
-	 ((gf(i+1,j,k+1) \
-	 + gf(i-1,j,k-1) \
-	 - gf(i+1,j,k-1) \
-	 - gf(i-1,j,k+1)) * hdxi * hdzi )
-
-/* first derivatives, centered, 2nd order */
-
-#define D1_c2(gf,i,j,k)                       \
-	 ((gf(i+1,j,k) \
-	 - gf(i-1,j,k)) * hdxi)
-
-#define D2_c2(gf,i,j,k)                       \
-	 ((gf(i,j+1,k) \
-	 - gf(i,j-1,k)) * hdyi)
-
-#define D3_c2(gf,i,j,k)                       \
-	 ((gf(i,j,k+1) \
-	 - gf(i,j,k-1)) * hdzi)
-	
-#else
-
-
-#define D11_c2(gf,i,j,k)                         \
-	 ((   gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k)] \
-	 - 2.*gf[CCTK_GFINDEX3D(cctkGH,i,  j,k)] \
-	 +    gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k)]) * dxi * dxi )
-
-#define D22_c2(gf,i,j,k)                         \
-	 ((   gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k)] \
-	 - 2.*gf[CCTK_GFINDEX3D(cctkGH,i,j,  k)] \
-	 +    gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k)]) * dyi * dyi )
-
-#define D33_c2(gf,i,j,k)                         \
-	 ((   gf[CCTK_GFINDEX3D(cctkGH,i,j,k+1)] \
-	 - 2.*gf[CCTK_GFINDEX3D(cctkGH,i,j,k  )] \
-	 +    gf[CCTK_GFINDEX3D(cctkGH,i,j,k-1)]) * dzi * dzi )
-
-#define D21_c2(gf,i,j,k)                        \
-	 ((gf[CCTK_GFINDEX3D(cctkGH,i+1,j+1,k)] \
-	 + gf[CCTK_GFINDEX3D(cctkGH,i-1,j-1,k)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i+1,j-1,k)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i-1,j+1,k)]) * hdxi * hdyi )
-
-#define D32_c2(gf,i,j,k)                        \
-	 ((gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k+1)] \
-	 + gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k-1)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k-1)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k+1)]) * hdyi * hdzi )
-
-#define D31_c2(gf,i,j,k)                        \
-	 ((gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k+1)] \
-	 + gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k-1)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k-1)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k+1)]) * hdxi * hdzi )
-
-/* first derivatives, centered, 2nd order */
-
-#define D1_c2(gf,i,j,k)                       \
-	 ((gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k)]) * hdxi)
-
-#define D2_c2(gf,i,j,k)                       \
-	 ((gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k)]) * hdyi)
-
-#define D3_c2(gf,i,j,k)                       \
-	 ((gf[CCTK_GFINDEX3D(cctkGH,i,j,k+1)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i,j,k-1)]) * hdzi)
-
-#endif	
-
-	
-/*  c4                  */
-/*                      */
-/*  4th order centered  */
-/*                      */
-
-/* second derivatives, centered, 4th order */
-
-#ifndef KRANC_C
-
-#define D11_c4(gf,i,j,k)                    \
-  ((-       gf(i+2,j,k)			      \
-	   + 16. * gf(i+1,j,k) \
-	   - 30. * gf(i,  j,k) \
-	   + 16. * gf(i-1,j,k) \
-	   -       gf(i-2,j,k)) * dxi * dxi / 12.)
-
-
-#define D22_c4(gf,i,j,k)                    \
-	 ((-       gf(i,j+2,k) \
-	   + 16. * gf(i,j+1,k) \
-	   - 30. * gf(i,j,  k) \
-	   + 16. * gf(i,j-1,k) \
-	   -       gf(i,j-2,k)) * dyi * dyi / 12.)
-
-
-#define D33_c4(gf,i,j,k)                              \
-	 ((-       gf(i,j,k+2) \
-	   + 16. * gf(i,j,k+1) \
-	   - 30. * gf(i,j,k  ) \
-	   + 16. * gf(i,j,k-1) \
-	   -       gf(i,j,k-2)) * dzi * dzi / 12.)
-
-#define D21_c4(gf,i,j,k)                                    \
-	     ((-       gf(i+2,j+2,k) \
-               +       gf(i+2,j-2,k) \
-	       +       gf(i-2,j+2,k) \
-               -       gf(i-2,j-2,k) \
-               + 16. * gf(i+1,j+1,k) \
-    	       - 16. * gf(i+1,j-1,k) \
-	       - 16. * gf(i-1,j+1,k) \
-	       + 16. * gf(i-1,j-1,k)) * dxi * dyi / 48.)
-
-#define D31_c4(gf,i,j,k)                                    \
-             ((-       gf(i+2,j,k+2) \
-               +       gf(i+2,j,k-2) \
-               +       gf(i-2,j,k+2) \
-               -       gf(i-2,j,k-2) \
-               + 16. * gf(i+1,j,k+1) \
-               - 16. * gf(i+1,j,k-1) \
-               - 16. * gf(i-1,j,k+1) \
-               + 16. * gf(i-1,j,k-1)) * dxi * dzi / 48.)
-
-	
-#define D32_c4(gf,i,j,k)                                    \
-             ((-       gf(i,j+2,k+2) \
-               +       gf(i,j+2,k-2) \
-               +       gf(i,j-2,k+2) \
-               -       gf(i,j-2,k-2) \
-               + 16. * gf(i,j+1,k+1) \
-               - 16. * gf(i,j+1,k-1) \
-               - 16. * gf(i,j-1,k+1) \
-               + 16. * gf(i,j-1,k-1)) * dzi * dyi / 48.)
-
-
-/* first derivatives, centered, 4th order */
-
-#define D1_c4(gf,i,j,k)                            \
-       ((-      gf(i+2,j,k) \
-         + 8. * gf(i+1,j,k) \
-         - 8. * gf(i-1,j,k) \
-	 +      gf(i-2,j,k)) * (dxi / 12.))
-
-#define D2_c4(gf,i,j,k)                            \
-       ((-      gf(i,j+2,k) \
-         + 8. * gf(i,j+1,k) \
-         - 8. * gf(i,j-1,k) \
-	 +      gf(i,j-2,k)) * (dyi / 12.))
-
-#define D3_c4(gf,i,j,k)                            \
-       ((-      gf(i,j,k+2) \
-         + 8. * gf(i,j,k+1) \
-         - 8. * gf(i,j,k-1) \
-	 +      gf(i,j,k-2)) * (dzi / 12.))
-
-
-#else
-
-#define D11_c4(gf,i,j,k)                    \
-	 ((-       gf[CCTK_GFINDEX3D(cctkGH,i+2,j,k)] \
-	   + 16. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k)] \
-	   - 30. * gf[CCTK_GFINDEX3D(cctkGH,i,  j,k)] \
-	   + 16. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k)] \
-	   -       gf[CCTK_GFINDEX3D(cctkGH,i-2,j,k)]) * dxi * dxi / 12.)
-
-
-#define D22_c4(gf,i,j,k)                    \
-	 ((-       gf[CCTK_GFINDEX3D(cctkGH,i,j+2,k)] \
-	   + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k)] \
-	   - 30. * gf[CCTK_GFINDEX3D(cctkGH,i,j,  k)] \
-	   + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k)] \
-	   -       gf[CCTK_GFINDEX3D(cctkGH,i,j-2,k)]) * dyi * dyi / 12.)
-
-
-#define D33_c4(gf,i,j,k)                              \
-	 ((-       gf[CCTK_GFINDEX3D(cctkGH,i,j,k+2)] \
-	   + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j,k+1)] \
-	   - 30. * gf[CCTK_GFINDEX3D(cctkGH,i,j,k  )] \
-	   + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j,k-1)] \
-	   -       gf[CCTK_GFINDEX3D(cctkGH,i,j,k-2)]) * dzi * dzi / 12.)
-
-#define D21_c4(gf,i,j,k)                                    \
-	     ((-       gf[CCTK_GFINDEX3D(cctkGH,i+2,j+2,k)] \
-               +       gf[CCTK_GFINDEX3D(cctkGH,i+2,j-2,k)] \
-	       +       gf[CCTK_GFINDEX3D(cctkGH,i-2,j+2,k)] \
-               -       gf[CCTK_GFINDEX3D(cctkGH,i-2,j-2,k)] \
-               + 16. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j+1,k)] \
-    	       - 16. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j-1,k)] \
-	       - 16. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j+1,k)] \
-	       + 16. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j-1,k)]) * dxi * dyi / 48.)
-
-#define D31_c4(gf,i,j,k)                                    \
-             ((-       gf[CCTK_GFINDEX3D(cctkGH,i+2,j,k+2)] \
-               +       gf[CCTK_GFINDEX3D(cctkGH,i+2,j,k-2)] \
-               +       gf[CCTK_GFINDEX3D(cctkGH,i-2,j,k+2)] \
-               -       gf[CCTK_GFINDEX3D(cctkGH,i-2,j,k-2)] \
-               + 16. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k+1)] \
-               - 16. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k-1)] \
-               - 16. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k+1)] \
-               + 16. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k-1)]) * dxi * dzi / 48.)
-
-	
-#define D32_c4(gf,i,j,k)                                    \
-             ((-       gf[CCTK_GFINDEX3D(cctkGH,i,j+2,k+2)] \
-               +       gf[CCTK_GFINDEX3D(cctkGH,i,j+2,k-2)] \
-               +       gf[CCTK_GFINDEX3D(cctkGH,i,j-2,k+2)] \
-               -       gf[CCTK_GFINDEX3D(cctkGH,i,j-2,k-2)] \
-               + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k+1)] \
-               - 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k-1)] \
-               - 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k+1)] \
-               + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k-1)]) * dzi * dyi / 48.)
-
-
-/* first derivatives, centered, 4th order */
-
-#define D1_c4(gf,i,j,k)                            \
-       ((-      gf[CCTK_GFINDEX3D(cctkGH,i+2,j,k)] \
-         + 8. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k)] \
-         - 8. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k)] \
-	 +      gf[CCTK_GFINDEX3D(cctkGH,i-2,j,k)]) * (dxi / 12.))
-
-#define D2_c4(gf,i,j,k)                            \
-       ((-      gf[CCTK_GFINDEX3D(cctkGH,i,j+2,k)] \
-         + 8. * gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k)] \
-         - 8. * gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k)] \
-	 +      gf[CCTK_GFINDEX3D(cctkGH,i,j-2,k)]) * (dyi / 12.))
-
-#define D3_c4(gf,i,j,k)                            \
-       ((-      gf[CCTK_GFINDEX3D(cctkGH,i,j,k+2)] \
-         + 8. * gf[CCTK_GFINDEX3D(cctkGH,i,j,k+1)] \
-         - 8. * gf[CCTK_GFINDEX3D(cctkGH,i,j,k-1)] \
-	 +      gf[CCTK_GFINDEX3D(cctkGH,i,j,k-2)]) * (dzi / 12.))
-
-#endif
-
-/*****************************************************/
-/*                                                    */
-/*         DERIVED METHODS                            */
-/*                                                    */
-/******************************************************/
-
-
-/* blend c2 and c4 */
-/* second derivatives */
-#define D11_c2c4(gf,i,j,k) (fdweight_c2*D11_c2(gf,i,j,k) + fdweight_c4*D11_c4(gf,i,j,k))
-#define D22_c2c4(gf,i,j,k) (fdweight_c2*D22_c2(gf,i,j,k) + fdweight_c4*D22_c4(gf,i,j,k))
-#define D33_c2c4(gf,i,j,k) (fdweight_c2*D33_c2(gf,i,j,k) + fdweight_c4*D33_c4(gf,i,j,k))
-#define D21_c2c4(gf,i,j,k) (fdweight_c2*D21_c2(gf,i,j,k) + fdweight_c4*D21_c4(gf,i,j,k))
-#define D32_c2c4(gf,i,j,k) (fdweight_c2*D32_c2(gf,i,j,k) + fdweight_c4*D32_c4(gf,i,j,k))
-#define D31_c2c4(gf,i,j,k) (fdweight_c2*D31_c2(gf,i,j,k) + fdweight_c4*D31_c4(gf,i,j,k))
-
-/* first derivatives */
-#define D1_c2c4(gf,i,j,k)  (fdweight_c2*D1_c2(gf, i,j,k) + fdweight_c4*D1_c4(gf,i,j,k))
-#define D2_c2c4(gf,i,j,k)  (fdweight_c2*D2_c2(gf, i,j,k) + fdweight_c4*D2_c4(gf,i,j,k))
-#define D3_c2c4(gf,i,j,k)  (fdweight_c2*D3_c2(gf, i,j,k) + fdweight_c4*D3_c4(gf,i,j,k))
-
-
-/*****************************************************/
-/*                                                    */
-/*  Poor man's one-sided derivatives                  */
-/*                                                    */
-/******************************************************/
-
-
-#define Dplus1(gf,i,j,k)   string(Dplus1,gf)
-#define Dplus2(gf,i,j,k)   string(Dplus2,gf)
-#define Dplus3(gf,i,j,k)   string(Dplus3,gf)
-
-
-#define Dplus1gf(gf,i,j,k)   Dplus1x(gf, i,j,k)
-#define Dplus2gf(gf,i,j,k)   Dplus2x(gf, i,j,k)
-#define Dplus3gf(gf,i,j,k)   Dplus3x(gf, i,j,k)
-
-#ifdef KRANC_C
-
-#define Dplus1x(gf,i,j,k)                     \
-	 ((gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dxi)
-
-#define Dplus2x(gf,i,j,k)                     \
-	 ((gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dyi)
-
-#define Dplus3x(gf,i,j,k)                     \
-	 ((gf[CCTK_GFINDEX3D(cctkGH,i,j,k+1)] \
-	 - gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dzi)
-
-#else
-
-#define Dplus1x(gf,i,j,k)  ( ( gf(i+1, j,     k)     - gf(i,j,k) ) * dxi )
-#define Dplus2x(gf,i,j,k)  ( ( gf(i,   j + 1, k)     - gf(i,j,k) ) * dxi )
-#define Dplus3x(gf,i,j,k)  ( ( gf(i,   j,     k + 1) - gf(i,j,k) ) * dxi )
-
-#endif
-
 #ifdef KRANC_C
 int sgn(CCTK_REAL x);
 
-#define Dupwind1(gf,dir,i,j,k) ((dir * gf[CCTK_GFINDEX3D(cctkGH,i+dir,j,k)] \
-	 - dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dxi)
-#define Dupwind2(gf,dir,i,j,k) ((dir * gf[CCTK_GFINDEX3D(cctkGH,i,j+dir,k)] \
-	 - dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dxi)
-#define Dupwind3(gf,dir,i,j,k) ((dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k+dir)] \
-	 - dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dxi)
-
 int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH);
 
-/* int GenericFD_BoundaryWidthTable(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,
@@ -688,24 +66,8 @@ void GenericFD_GetBoundaryInfo(cGH const * restrict cctkGH,
 			       int * restrict is_physbnd,
                                int * restrict is_ipbnd);
 
-#if 0
-/* Finite differencing near boundaries */
-
-/* The array var is to be accessed at the location
-   [i+ioff,j+joff,k+koff].  idir,jdir,kdir specify whether there is a
-   lower (dir<0), upper (dir>0), or no boundary nearby.  If a boundary
-   is in the way, the value 0 is returned instead of the array
-   content.  */
-#define CCTK_GFACCESS3D(cctkGH, var, i,j,k, ioff,joff,koff, idir,jdir,kdir) \
-  (((idir)<0 && (ioff)<0) ||                                            \
-   ((jdir)<0 && (joff)<0) ||                                            \
-   ((kdir)<0 && (koff)<0) ||                                            \
-   ((idir)>0 && (ioff)>0) ||                                            \
-   ((jdir)>0 && (joff)>0) ||                                            \
-   ((kdir)>0 && (koff)>0) ||                                            \
-   ? 0                                                                  \
-   : (var)[CCTK_GFINDEX3D((cctkGH), (i)+(ioff),(j)+(joff),(k)+(koff))])
-#endif
+void GenericFD_AssertGroupStorage(cGH const * restrict const cctkGH, const char *calc,
+                                  int ngroups, const char *group_names[]);
 
 /* Summation by parts */
 
@@ -790,40 +152,13 @@ void GenericFD_LoopOverBoundary(cGH const * restrict cctkGH, Kranc_Calculation c
 void GenericFD_LoopOverBoundaryWithGhosts(cGH const * restrict cctkGH, Kranc_Calculation calc);
 void GenericFD_LoopOverInterior(cGH const * restrict cctkGH, Kranc_Calculation calc);
 
+void GenericFD_GroupDataPointers(cGH const * restrict const cctkGH, const char *group_name,
+                                 int nvars, CCTK_REAL const *restrict *ptrs);
+void GenericFD_EnsureStencilFits(cGH const * restrict const cctkGH, const char *calc, int ni, int nj, int nk);
 
 
-/* Vectorisation of memory accesses  */
-
-#include <stdlib.h>
-#include <cctk.h>
-
-#if defined(__SSE2__) && defined(CCTK_REAL_PRECISION_8)
-
-#include <emmintrin.h>
-
-/* A vector type corresponding to CCTK_REAL */
-typedef __m128d CCTK_REAL_VEC;
-
-/* Really only SSE is required, but there doesn't seem to be a
-   preprocessing flag to check for this */
-#elif defined(__SSE2__) && defined(CCTK_REAL_PRECISION_4)
-
-#include <emmintrin.h>
-
-/* A vector type corresponding to CCTK_REAL */
-typedef __m128 CCTK_REAL_VEC;
-
-#else
-
-/* There is no vector type corresponding to CCTK_REAL */
-typedef CCTK_REAL CCTK_REAL_VEC;
-
+#ifdef __cplusplus
+} /* extern "C" */
 #endif
 
-/* The number of vector elements in a CCTK_REAL_VEC */
-static
-size_t const CCTK_REAL_VEC_SIZE = sizeof(CCTK_REAL_VEC) / sizeof(CCTK_REAL);
-
-
-
 #endif
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h
index 86b70eb..e48daa2 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h
@@ -1,19 +1,19 @@
  
-#define Power(x, y)   (pow((x), (y)))
+#define Power(x, y)   (pow(x,y))
 #define Sqrt(x)       (sqrt(x))
 
 
 #ifdef KRANC_C
-#define Abs(x)        (fabs(x))
-#define Min(x, y)     (fmin((x), (y)))
-#define Min3(x, y, z)     (fmin(fmin((x), (y)), (z)))
-#define Max(x, y)     (fmax((x), (y)))
-#define IfThen(x,y,z) ((x) ? (y) : (z))
+#  define Abs(x)        (fabs(x))
+#  define Min(x, y)     (fmin(x,y))
+#  define Min3(x, y, z) (fmin(fmin((x), (y)), (z)))
+#  define Max(x, y)     (fmax(x,y))
+#  define IfThen(x,y,z) ((x) ? (y) : (z))
 #else
-#define Abs(x)        (abs(x))
-#define Min(x, y)     (min((x), (y)))
-#define Max(x, y)     (max((x), (y)))
-/* IfThen cannot be expressed in Fortran */
+#  define Abs(x)        (abs(x))
+#  define Min(x, y)     (min(x,y))
+#  define Max(x, y)     (max(x,y))
+#  define IfThen(x,y,z) ((x)*(y) + (1-(x))*(z))
 #endif
 
 #ifdef KRANC_C
@@ -42,11 +42,52 @@
 #define Tanh(x)       (tanh(x))
 
 #ifdef KRANC_C
-#define E               M_E
-#define Pi              M_PI
+#  define Sign(x)     (copysign(1.0,(x)))
+#  define ToReal(x)   ((CCTK_REAL)(x))
 #else
-#define E               2.71828182845904523536029d0
-#define Pi              3.14159265358979323846264d0
+#  define Sign(x)     (sgn(x))
+#  define ToReal(x)   (real((x),kind(khalf)))
 #endif
 
-#define StepFunction(x) ( (x) > 0 ? 1 : 0 )
+#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
+
+#else
+#  define E           2.71828182845904523536029d0
+#  define Pi          3.14159265358979323846264d0
+#endif
+
+#define StepFunction(x) ((x)>0)
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/ParamCheck.F90 b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/ParamCheck.F90
deleted file mode 100644
index c2baf82..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/ParamCheck.F90
+++ /dev/null
@@ -1,96 +0,0 @@
-/*@@                                         
-   @file      GenericFD/src/ParamCheck.F90
-   @date      October 20 2004
-   @author    S. Husa                           
-   @desc                                        
-   Check consistency of parameters associated with stencil widths
-
-   $Header$                                  
-
-   @enddesc                                     
- @@*/                                           
-                                                
-/*  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
-*/
-
-#define KRANC_FORTRAN
-                                           
-#include "cctk.h"         
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-#include "GenericFD.h"
-
-/*@@                                         
-   @routine   GenericFD_ParamCheck 
-   @date      October 20 2004
-   @author    S. Husa                           
-   @desc                     
-   Check consistency of parameters associated with stencil widths                   
-   
-   @enddesc                                     
-   @calls                                       
-   @calledby                                    
-   @history                                     
-   @endhistory                                  
- @@*/                                           
-                                                
-                                              
-
-SUBROUTINE GenericFD_ParamCheck(CCTK_ARGUMENTS)
-
-implicit none
-
-DECLARE_CCTK_ARGUMENTS
-DECLARE_CCTK_PARAMETERS
-
-if (stencil_width < 0) then
-  call CCTK_WARN(0, "stencil_width < 0  - set GenericFD::stencil_width > 0 in par file!")
-endif
-
-if ((stencil_width_x < 0).AND.(stencil_width < 0)) then
-  call CCTK_WARN(0, "stencil_width_x < 0, set GenericFD::stencil_width_x (or stencil_width) > 0!")
-endif
-
-if ((stencil_width_y < 0).AND.(stencil_width < 0)) then
-  call CCTK_WARN(0, "stencil_width_y < 0, set GenericFD::stencil_width_x (or stencil_width) > 0!")
-endif
-
-if ((stencil_width_z < 0).AND.(stencil_width < 0)) then
-  call CCTK_WARN(0, "stencil_width_z < 0, set GenericFD::stencil_width_x (or stencil_width) > 0!")
-endif
-
-
-if (stencil_width > maxval(cctk_nghostzones)) then
- call CCTK_WARN(0, "stencil_width is larger than max(cctk_nghostzones)!")
-endif
-
-if (stencil_width_x > cctk_nghostzones(1)) then
- call CCTK_WARN(0, "stencil_width is smaller than cctk_nghostzones(1)!")
-endif
-
-if (stencil_width_y > cctk_nghostzones(2)) then
- call CCTK_WARN(0, "stencil_width is smaller than cctk_nghostzones(2)!")
-endif
-
-if (stencil_width_z > cctk_nghostzones(3)) then
- call CCTK_WARN(0, "stencil_width is smaller than cctk_nghostzones(3)!")
-endif
-
-
-END SUBROUTINE GenericFD_ParamCheck
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/Startup.c b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/Startup.c
deleted file mode 100644
index 78e9915..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/Startup.c
+++ /dev/null
@@ -1,60 +0,0 @@
-/*@@                                         
-   @file      GenericFD/src/Startup.c
-   @date      June 16 2002
-   @author    S. Husa                           
-   @desc                                        
-   Register Banner - straight copy of WaveToy   
-
-   $Header$                                  
-
-   @enddesc                                     
- @@*/                                           
-                                                
-/*  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
-*/
-                                           
-#include "cctk.h"                           
-#include "GenericFD.h"
-
-int GenericFD_Startup(void);
-
-/*@@                                         
-   @routine   GenericFD_Startup 
-   @date      June 16 2002
-   @author    S. Husa                           
-   @desc                                        
-   
-   @enddesc                                     
-   @calls                                       
-   @calledby                                    
-   @history                                     
-   @endhistory                                  
- @@*/                                           
-                                                
-                                              
-
-int GenericFD_Startup(void)
-{                                          
-    const char *banner = "GenericFD: generic finite differencing";
-    CCTK_RegisterBanner(banner);              
-
-
-    CCTK_INFO(FD_METHOD_DESC);
-    return 0; 
- }
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/make.code.defn b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/make.code.defn
index fde3275..c9ef2f1 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/make.code.defn
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/make.code.defn
@@ -2,8 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = Startup.c GenericFD.c 
-#ParamCheck.F90
+SRCS = GenericFD.c
 
 # Subdirectories containing source files
 SUBDIRS =
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/testmacros.c b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/testmacros.c
deleted file mode 100644
index 3208ea0..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/testmacros.c
+++ /dev/null
@@ -1,23 +0,0 @@
-/* process me with something like
- *  cpp -DFD_C2 testmacros.c
- */
-
-
-#include "GenericFD.h"
-
-
-dummy =  D11loc(testfunc,i,j,k);
-
-dummy =  D11gf(testfunc,i,j,k);
-
-dummy =  D11_c2c4(testfunc,i,j,k);
-
-dummy = D11(testfunc,i,j,k);
-
-dummy = D31(testfunc,i,j,k);
-
-dummy = D1(testfunc,i,j,k);
-
-
-DECLAREVARS
-
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/README b/Auxiliary/Cactus/KrancNumericalTools/Perturb/README
deleted file mode 100644
index 899015f..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/README
+++ /dev/null
@@ -1,5 +0,0 @@
-
-Purpose of the thorn:
-
-This thorn adds a perturbation to grid variables at initial data and/or
-boundaries, in the style of Denis' noise thorn.
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/interface.ccl b/Auxiliary/Cactus/KrancNumericalTools/Perturb/interface.ccl
deleted file mode 100644
index a745e5c..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/interface.ccl
+++ /dev/null
@@ -1,7 +0,0 @@
-# Interface definition for thorn Perturb
-# $Header$
-
-IMPLEMENTS: perturb
-
-INHERITS: grid
-
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/param.ccl b/Auxiliary/Cactus/KrancNumericalTools/Perturb/param.ccl
deleted file mode 100644
index 8e87dc6..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/param.ccl
+++ /dev/null
@@ -1,45 +0,0 @@
-# Parameter definitions for thorn Perturb
-# $Header$
-
-#------------------------------------------------------------------------------
-# Private:
-#------------------------------------------------------------------------------
-private:
-
-BOOLEAN apply_id_perturb "Add perturbation to initial data"
-{
-} "no"
-
-BOOLEAN apply_bc_perturb "Add perturbation to boundary data"
-{
-} "no"
-
-STRING id_vars "Initial data variables to modify with perturbation"
-{
-  .* :: "A regex which matches everything"
-} ""
-
-STRING bc_vars "Variables to modify with perturbation at boundary"
-{
-  .* :: "A regex which matches everything"
-} ""
-
-BOOLEAN perturb_boundaries[6] "At which boundaries to apply perturbation"
-{
-} "yes"
-
-INT perturb_stencil[3] "Number of boundary points"
-{
-  0:* :: "0:*"
-} 1
-
-REAL amplitude "Amplitude of perturbation data"
-{
-  0: :: "Positive number"
-} 0.000001
-
-REAL period "period of perturbation data"
-{
-  0: :: "Positive number"
-} 1.0
-
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/schedule.ccl b/Auxiliary/Cactus/KrancNumericalTools/Perturb/schedule.ccl
deleted file mode 100644
index 268f81d..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/schedule.ccl
+++ /dev/null
@@ -1,21 +0,0 @@
-# Schedule definitions for thorn Perturb
-# $Header$
-
-if (apply_id_perturb) {
-  SCHEDULE id_perturb AT CCTK_POSTINITIAL
-  {
-    LANG: C
-  } "Add perturb to initial data"
-}
-
-if (apply_bc_perturb) {
-  SCHEDULE bc_perturb AT CCTK_POSTSTEP
-  {
-     LANG: C
-  } "Add perturb to boundary condition"
-
-  SCHEDULE bc_perturb AT CCTK_POSTRESTRICT
-  {
-     LANG: C
-  } "Add perturb to boundary condition"
-}
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/bc_perturb.c b/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/bc_perturb.c
deleted file mode 100644
index bf56039..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/bc_perturb.c
+++ /dev/null
@@ -1,322 +0,0 @@
-/* $Header$ */
-
-#include <assert.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-
-
-#include "perturb.h"
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-#include "cctk_FortranString.h"
-
-#include "Symmetry.h"
-
-
-
-/* #define DEBUG_BOUNDARY 1 */
-
-static int ApplyBndperturb (const cGH *GH,
-			  int stencil_dir,
-			  const int *stencil,
-			  int dir,
-			  int first_var,
-			  int num_vars);
-
-
-
-int
-BndperturbVI (const cGH *GH, const int *stencil, int vi)
-{
-  int retval;
-  retval = ApplyBndperturb (GH, -1, stencil, 0, vi, 1);
-  return retval;
-}
-
-void
-CCTK_FCALL CCTK_FNAME (BndperturbVI) (int *ierr, const cGH **GH,
-				    const int *stencil, const int *vi)
-{
-  *ierr = BndperturbVI (*GH, stencil, *vi);
-}
-
-int
-BndperturbVN (const cGH *GH, const int *stencil, const char *vn)
-{
-  int vi, retval;
-  vi = CCTK_VarIndex(vn);
-  retval = BndperturbVI (GH, stencil, vi);
-  return retval;
-}
-
-void
-CCTK_FCALL CCTK_FNAME (BndperturbVN) (int *ierr, const cGH **GH,
-				    const int *stencil, ONE_FORTSTRING_ARG)
-{
-  ONE_FORTSTRING_CREATE (vn);
-  *ierr = BndperturbVN (*GH, stencil, vn);
-  free (vn);
-}
-
-
-
-int
-BndperturbGI (const cGH *GH, const int *stencil, int gi)
-{
-  int first_vi, retval;
-
-  first_vi = CCTK_FirstVarIndexI (gi);
-  if (first_vi >= 0)
-    {
-      retval = ApplyBndperturb (GH, -1, stencil, 0, first_vi,
-			      CCTK_NumVarsInGroupI (gi));
-    }
-  else
-    {
-      CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
-		  "Invalid group index %d in BndperturbGI", gi);
-      retval = -1;
-    }
-
-  return (retval);
-}
-
-
-
-void
-CCTK_FCALL CCTK_FNAME (BndperturbGI) (int *ierr, const cGH **GH,
-				    const int *stencil, const int *gi)
-{
-  *ierr = BndperturbGI (*GH, stencil, *gi);
-}
-
-
-int
-BndperturbGN (const cGH *GH, const int *stencil, const char *gn)
-{
-  int gi, retval;
-
-  gi = CCTK_GroupIndex (gn);
-  if (gi >= 0)
-    {
-      retval = BndperturbGI (GH, stencil, gi);
-    }
-  else
-    {
-      CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
-		  "Invalid group name '%s' in BndperturbGN", gn);
-      retval = -1;
-    }
-
-  return (retval);
-}
-
-void CCTK_FCALL CCTK_FNAME (BndperturbGN)
-     (int *ierr,
-      const cGH **GH,
-      const int *stencil,
-      ONE_FORTSTRING_ARG)
-{
-  ONE_FORTSTRING_CREATE (gn)
-    *ierr = BndperturbGN (*GH, stencil, gn);
-  free (gn);
-}
-
-
-
-
-
-static int ApplyBndperturb (const cGH *GH,
-			  int stencil_dir,
-			  const int *stencil,
-			  int dir,
-			  int first_var,
-			  int num_vars)
-{
-  DECLARE_CCTK_PARAMETERS;
-  int i, j, k;
-  int var, vtypesize, gindex, gdim, timelvl;
-  int doBC[2*MAXDIM], dstag[MAXDIM], lsh[MAXDIM], lssh[MAXDIM];
-  SymmetryGHex *sGHex;
-  int type;
-
-  /* This argument is unused an undocumented; better make sure people
-     don't try to use it for something.  */
-  assert (stencil_dir == -1);
-
-  /* get the group index of the variables */
-  gindex = CCTK_GroupIndexFromVarI (first_var);
-
-  /* get the number of dimensions and the size of the variables' type */
-  gdim      = CCTK_GroupDimI (gindex);
-  vtypesize = CCTK_VarTypeSize (CCTK_VarTypeI (first_var));
-
-  /* make sure we can deal with this number of dimensions */
-  if (gdim > MAXDIM)
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		  "ApplyBndperturb: Variable dimension of %d not supported", gdim);
-      return (-1);
-    }
-
-  /* check the direction parameter */
-  if (abs (dir) > gdim)
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		  "ApplyBndperturb: direction %d greater than dimension %d",
-		  dir, gdim);
-      return (-2);
-    }
-
-  /* initialize arrays for variables with less dimensions than MAXDIM
-     so that we can use the INDEX_3D macro later on */
-  for (i = gdim; i < MAXDIM; i++)
-    {
-      lsh[i]  = 1;
-      lssh[i] = 1;
-    }
-
-  /* get the directional staggering of the group */
-  CCTK_GroupStaggerDirArrayGI (dstag, gdim, gindex);
-
-  /* get the current timelevel */
-  timelvl = 0;
-
-  /* see if we have a symmetry array */
-  sGHex = (SymmetryGHex *) CCTK_GHExtension (GH, "Symmetry");
-
-
-  /* now loop over all variables */
-  for (var = first_var; var < first_var + num_vars; var++)
-    {
-      /* Apply condition if:
-	 + boundary is not a symmetry boundary (no symmetry or unset(=unsed))
-	 + boundary is a physical boundary
-	 + have enough grid points
-      */
-      memset (doBC, 1, sizeof (doBC));
-      if (sGHex)
-	{
-	  for (i = 0; i < 2 * gdim; i++)
-	    {
-	      doBC[i] = sGHex->GFSym[var][i] == GFSYM_NOSYM ||
-		sGHex->GFSym[var][i] == GFSYM_UNSET;
-	    }
-	}
-      for (i = 0; i < gdim; i++)
-	{
-	  lsh[i]       = GH->cctk_lsh[i];
-	  lssh[i]      = GH->cctk_lssh[CCTK_LSSH_IDX (dstag[i], i)];
-	  doBC[i*2]   &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2];
-	  doBC[i*2+1] &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2+1];
-	  if (dir != 0)
-	    {
-	      doBC[i*2]   &= i+1 == -dir;
-	      doBC[i*2+1] &= i+1 ==  dir;
-	    }
-	}
-
-      /* now apply the boundaries face by face */
-      if (gdim > 0)
-	{
-#ifdef DEBUG_BOUNDARY
-	  if (doBC[0])
-	    {
-	      printf("Boundary: Applying lower x perturb to boundary\n");
-	    }
-	  if (doBC[1])
-	    {
-	      printf("Boundary: Applying upper x perturb to boundary\n");
-	    }
-#endif /* DEBUG_BOUNDARY */
-	  /* lower x */
-	  BOUNDARY_PERTURB (doBC[0], stencil[0], lssh[1], lssh[2],
-			  i, j, k);
-	  /* upper x */
-	  BOUNDARY_PERTURB (doBC[1], stencil[0], lssh[1], lssh[2],
-			  lssh[0]-i-1, j, k);
-
-	}
-      if (gdim > 1)
-
-	{
-#ifdef DEBUG_BOUNDARY
-	  if (doBC[2])
-	    {
-	      printf("Boundary: Applying lower y perturb to boundary\n");
-	    }
-	  if (doBC[3])
-	    {
-	      printf("Boundary: Applying upper y perturb to boundary\n");
-	    }
-#endif /* DEBUG_BOUNDARY */
-	  /* lower y */
-	  BOUNDARY_PERTURB (doBC[2], lssh[0], stencil[1], lssh[2],
-			  i, j, k);
-	  /* upper y */
-	  BOUNDARY_PERTURB (doBC[3], lssh[0], stencil[1], lssh[2],
-			  i, lssh[1]-j-1, k);
-	}
-      if (gdim > 2)
-	{
-#ifdef DEBUG_BOUNDARY
-	  if (doBC[4])
-	    {
-	      printf("Boundary: Applying lower z perturb to boundary\n");
-	    }
-	  if (doBC[5])
-	    {
-	      printf("Boundary: Applying upper z perturb to boundary\n");
-	    }
-#endif /* DEBUG_BOUNDARY */
-	  /* lower z */
-	  BOUNDARY_PERTURB (doBC[4], lssh[0], lssh[1], stencil[2],
-			  i, j, k);
-	  /* upper z */
-	  BOUNDARY_PERTURB (doBC[5], lssh[0], lssh[1], stencil[2],
-			  i, j, lssh[2]-k-1);
-	}
-    }
-
-  return(0);
-}
-
-
-static void
-add_bc_perturb_to_var (int idx, const char* optstring, void* cctkGH)
-{
-  DECLARE_CCTK_PARAMETERS;
-  cGH* GH = cctkGH;
-  int sw[3];
-
-  /* Change type from CCTK_INT to int */
-  sw[0] = perturb_stencil[0];
-  sw[1] = perturb_stencil[1];
-  sw[2] = perturb_stencil[2];
-
-  if (perturb_boundaries[0]) ApplyBndperturb (GH, -1, sw, -1, idx, 1);
-  if (perturb_boundaries[1]) ApplyBndperturb (GH, -1, sw, +1, idx, 1);
-  if (perturb_boundaries[2]) ApplyBndperturb (GH, -1, sw, -2, idx, 1);
-  if (perturb_boundaries[3]) ApplyBndperturb (GH, -1, sw, +2, idx, 1);
-  if (perturb_boundaries[4]) ApplyBndperturb (GH, -1, sw, -3, idx, 1);
-  if (perturb_boundaries[5]) ApplyBndperturb (GH, -1, sw, +3, idx, 1);
-}
-
-
-void
-bc_perturb(CCTK_ARGUMENTS)
-{
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_PARAMETERS
-
-/*   Boundary_MakeSureThatTheSelectionIsEmpty(); */
-  if (CCTK_TraverseString(bc_vars, add_bc_perturb_to_var, cctkGH,
-      		    CCTK_GROUP_OR_VAR) < 0)
-    {
-      CCTK_WARN (1, "Failed to parse 'perturb::bc_vars' parameter");
-    }
-}
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/id_perturb.c b/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/id_perturb.c
deleted file mode 100644
index 07bd106..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/id_perturb.c
+++ /dev/null
@@ -1,49 +0,0 @@
-/*
-  $Header$
-*/
-
-#include <stdlib.h>
-
-#include "perturb.h"
-
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-
-void
-add_perturb_to_var (int idx, const char* optstring, void* cctkGH)
-{
-  DECLARE_CCTK_PARAMETERS;
-  int i, j, k, ijk;
-  CCTK_REAL* data;
-  cGH* GH = cctkGH;
-  int type;
-
-  data = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, idx);
-
-  for (k=1; k< GH->cctk_lsh[2]-1; ++k)
-    {
-      for (j=1; j< GH->cctk_lsh[1]-1; ++j)
-	{
-	  for (i=1; i< GH->cctk_lsh[0]-1; ++i)
-	    {
-	      ijk = CCTK_GFINDEX3D(GH, i, j, k);
-
-	      data[ijk] += amplitude * pow(-1, i/period) * pow(-1, j/period) * pow(-1, k/period);
-      	    }
-	}
-    }
-}
-
-
-void
-id_perturb(CCTK_ARGUMENTS)
-{
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_PARAMETERS
-
-  if (CCTK_TraverseString(id_vars, add_perturb_to_var, cctkGH,
-			  CCTK_GROUP_OR_VAR) < 0)
-    {
-      CCTK_WARN (1, "Failed to parse 'IDRandom::noisy_id_vars' parameter");
-    }
-}
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/make.code.defn b/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/make.code.defn
deleted file mode 100644
index 795dd51..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/make.code.defn
+++ /dev/null
@@ -1,8 +0,0 @@
-# Main make.code.defn file for thorn perturb
-# $Header$
-
-# Source files in this directory
-SRCS = id_perturb.c bc_perturb.c
-
-# Subdirectories containing source files
-SUBDIRS = 
diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/perturb.h b/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/perturb.h
deleted file mode 100644
index e7da060..0000000
--- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/perturb.h
+++ /dev/null
@@ -1,45 +0,0 @@
-/* $Header$ */
-
-#ifndef PERTURB_H
-#define PERTURB_H
-
-#include <assert.h>
-#include <stdlib.h>
-#include <string.h>
-#include <time.h>
-#include "cctk.h"
-
-
-/* constants */
-
-#define MAXDIM  3
-
-
-/* macros */
-
-#define RAND_VAL ((random()*(1.0/RAND_MAX)-0.5)*amplitude)
-
-#define BOUNDARY_PERTURB(doBC,                                            \
-                       iend, jend, kend,                                \
-                       ii, jj, kk)                                      \
-{                                                                       \
-  if (doBC)                                                             \
-  {                                                                     \
-    CCTK_REAL* v= CCTK_VarDataPtrI(GH, timelvl, var);                   \
-    for (k = 0; k < kend; k++)                                          \
-    {                                                                   \
-      for (j = 0; j < jend; j++)                                        \
-      {                                                                 \
-        for (i = 0; i < iend; i++)                                      \
-        {                                                               \
-          const int _index = CCTK_GFINDEX3D(GH, (ii), (jj), (kk));      \
-          v[_index] += RAND_VAL;                                        \
-        }                                                               \
-      }                                                                 \
-    }                                                                   \
-  }                                                                     \
-}
-
-
-
-#endif /* !define(PERTURB_H) */