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

7 files changed, 199 insertions, 904 deletions

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
-