Regenerate code

19 files changed, 552 insertions, 679 deletions

diff --git a/ML_BSSN/param.ccl b/ML_BSSN/param.ccl
index 4df76a8..8b629ba 100644
--- a/ML_BSSN/param.ccl
+++ b/ML_BSSN/param.ccl
@@ -249,12 +249,6 @@ CCTK_INT ML_BSSN_RHSStaticBoundary_calc_every "ML_BSSN_RHSStaticBoundary_calc_ev
 } 1
 
 restricted:
-CCTK_INT ML_BSSN_RHSRadiativeBoundary_calc_every "ML_BSSN_RHSRadiativeBoundary_calc_every"
-{
-  *:* :: ""
-} 1
-
-restricted:
 CCTK_INT ML_BSSN_enforce_calc_every "ML_BSSN_enforce_calc_every"
 {
   *:* :: ""
@@ -291,13 +285,13 @@ CCTK_INT ML_BSSN_convertToADMBaseFakeDtLapseShift_calc_every "ML_BSSN_convertToA
 } 1
 
 restricted:
-CCTK_INT ML_BSSN_constraints_calc_every "ML_BSSN_constraints_calc_every"
+CCTK_INT ML_BSSN_constraints1_calc_every "ML_BSSN_constraints1_calc_every"
 {
   *:* :: ""
 } 1
 
 restricted:
-CCTK_INT ML_BSSN_constraints_boundary_calc_every "ML_BSSN_constraints_boundary_calc_every"
+CCTK_INT ML_BSSN_constraints2_calc_every "ML_BSSN_constraints2_calc_every"
 {
   *:* :: ""
 } 1
@@ -339,12 +333,6 @@ CCTK_INT ML_BSSN_RHSStaticBoundary_calc_offset "ML_BSSN_RHSStaticBoundary_calc_o
 } 0
 
 restricted:
-CCTK_INT ML_BSSN_RHSRadiativeBoundary_calc_offset "ML_BSSN_RHSRadiativeBoundary_calc_offset"
-{
-  *:* :: ""
-} 0
-
-restricted:
 CCTK_INT ML_BSSN_enforce_calc_offset "ML_BSSN_enforce_calc_offset"
 {
   *:* :: ""
@@ -381,13 +369,13 @@ CCTK_INT ML_BSSN_convertToADMBaseFakeDtLapseShift_calc_offset "ML_BSSN_convertTo
 } 0
 
 restricted:
-CCTK_INT ML_BSSN_constraints_calc_offset "ML_BSSN_constraints_calc_offset"
+CCTK_INT ML_BSSN_constraints1_calc_offset "ML_BSSN_constraints1_calc_offset"
 {
   *:* :: ""
 } 0
 
 restricted:
-CCTK_INT ML_BSSN_constraints_boundary_calc_offset "ML_BSSN_constraints_boundary_calc_offset"
+CCTK_INT ML_BSSN_constraints2_calc_offset "ML_BSSN_constraints2_calc_offset"
 {
   *:* :: ""
 } 0
@@ -702,7 +690,7 @@ KEYWORD ML_curv_bound "Boundary condition to implement"
   "scalar" :: "Dirichlet boundary condition"
   "newrad" :: "Improved radiative boundary condition"
   "skip" :: "skip boundary condition code"
-} "skip"
+} "none"
 
 private:
 KEYWORD ML_dtlapse_bound "Boundary condition to implement"
@@ -714,7 +702,7 @@ KEYWORD ML_dtlapse_bound "Boundary condition to implement"
   "scalar" :: "Dirichlet boundary condition"
   "newrad" :: "Improved radiative boundary condition"
   "skip" :: "skip boundary condition code"
-} "skip"
+} "none"
 
 private:
 KEYWORD ML_dtshift_bound "Boundary condition to implement"
@@ -726,7 +714,7 @@ KEYWORD ML_dtshift_bound "Boundary condition to implement"
   "scalar" :: "Dirichlet boundary condition"
   "newrad" :: "Improved radiative boundary condition"
   "skip" :: "skip boundary condition code"
-} "skip"
+} "none"
 
 private:
 KEYWORD ML_Gamma_bound "Boundary condition to implement"
@@ -738,7 +726,7 @@ KEYWORD ML_Gamma_bound "Boundary condition to implement"
   "scalar" :: "Dirichlet boundary condition"
   "newrad" :: "Improved radiative boundary condition"
   "skip" :: "skip boundary condition code"
-} "skip"
+} "none"
 
 private:
 KEYWORD ML_lapse_bound "Boundary condition to implement"
@@ -750,7 +738,7 @@ KEYWORD ML_lapse_bound "Boundary condition to implement"
   "scalar" :: "Dirichlet boundary condition"
   "newrad" :: "Improved radiative boundary condition"
   "skip" :: "skip boundary condition code"
-} "skip"
+} "none"
 
 private:
 KEYWORD ML_log_confac_bound "Boundary condition to implement"
@@ -762,7 +750,7 @@ KEYWORD ML_log_confac_bound "Boundary condition to implement"
   "scalar" :: "Dirichlet boundary condition"
   "newrad" :: "Improved radiative boundary condition"
   "skip" :: "skip boundary condition code"
-} "skip"
+} "none"
 
 private:
 KEYWORD ML_metric_bound "Boundary condition to implement"
@@ -774,7 +762,7 @@ KEYWORD ML_metric_bound "Boundary condition to implement"
   "scalar" :: "Dirichlet boundary condition"
   "newrad" :: "Improved radiative boundary condition"
   "skip" :: "skip boundary condition code"
-} "skip"
+} "none"
 
 private:
 KEYWORD ML_shift_bound "Boundary condition to implement"
@@ -786,7 +774,7 @@ KEYWORD ML_shift_bound "Boundary condition to implement"
   "scalar" :: "Dirichlet boundary condition"
   "newrad" :: "Improved radiative boundary condition"
   "skip" :: "skip boundary condition code"
-} "skip"
+} "none"
 
 private:
 KEYWORD ML_trace_curv_bound "Boundary condition to implement"
@@ -798,7 +786,7 @@ KEYWORD ML_trace_curv_bound "Boundary condition to implement"
   "scalar" :: "Dirichlet boundary condition"
   "newrad" :: "Improved radiative boundary condition"
   "skip" :: "skip boundary condition code"
-} "skip"
+} "none"
 
 private:
 CCTK_REAL At11_bound_speed "characteristic speed at boundary"
diff --git a/ML_BSSN/schedule.ccl b/ML_BSSN/schedule.ccl
index 516fd9a..a6227b2 100644
--- a/ML_BSSN/schedule.ccl
+++ b/ML_BSSN/schedule.ccl
@@ -312,15 +312,6 @@ if (CCTK_EQUALS(my_rhs_boundary_condition, "static"))
   } "ML_BSSN_RHSStaticBoundary"
 }
 
-
-if (CCTK_EQUALS(my_rhs_boundary_condition, "radiative"))
-{
-  schedule ML_BSSN_RHSRadiativeBoundary IN MoL_CalcRHS
-  {
-    LANG: C
-  } "ML_BSSN_RHSRadiativeBoundary"
-}
-
 schedule ML_BSSN_enforce IN MoL_PostStep BEFORE ML_BSSN_SelectBoundConds
 {
   LANG: C
@@ -369,20 +360,82 @@ if (CCTK_EQUALS(dt_lapse_shift_method, "noLapseShiftAdvection"))
   } "ML_BSSN_convertToADMBaseFakeDtLapseShift"
 }
 
-schedule ML_BSSN_constraints IN ML_BSSN_constraintsCalcGroup
+schedule group ML_BSSN_constraints1_group in MoL_PseudoEvolution
+{
+  # no language specified
+} "ML_BSSN_constraints1"
+
+schedule ML_BSSN_constraints1 in ML_BSSN_constraints1_group
+{
+  LANG: C
+} "ML_BSSN_constraints1"
+
+schedule ML_BSSN_constraints1_SelectBCs in ML_BSSN_constraints1_bc_group
 {
   LANG: C
+  OPTIONS: level
   SYNC: ML_cons_detg
   SYNC: ML_cons_Gamma
-  SYNC: ML_cons_traceA
   SYNC: ML_Ham
-  SYNC: ML_mom
-} "ML_BSSN_constraints"
+} "ML_BSSN_constraints1_SelectBCs"
+
+schedule group ApplyBCs as ML_BSSN_constraints1_ApplyBCs in ML_BSSN_constraints1_bc_group after ML_BSSN_constraints1_SelectBCs
+{
+  # no language specified
+} "Apply BCs for groups set in ML_BSSN_constraints1"
+
+schedule group ML_BSSN_constraints1_bc_group in ML_BSSN_constraints1_group after ML_BSSN_constraints1
+{
+  # no language specified
+} "ML_BSSN_constraints1"
+
+schedule group ML_BSSN_constraints1_bc_group at CCTK_POSTRESTRICT
+{
+  # no language specified
+} "ML_BSSN_constraints1"
 
-schedule ML_BSSN_constraints_boundary IN ML_BSSN_constraintsCalcGroup AFTER ML_BSSN_constraints
+schedule group ML_BSSN_constraints1_bc_group at CCTK_POSTRESTRICTINITIAL
+{
+  # no language specified
+} "ML_BSSN_constraints1"
+
+schedule group ML_BSSN_constraints2_group in MoL_PseudoEvolution
+{
+  # no language specified
+} "ML_BSSN_constraints2"
+
+schedule ML_BSSN_constraints2 in ML_BSSN_constraints2_group
 {
   LANG: C
-} "ML_BSSN_constraints_boundary"
+} "ML_BSSN_constraints2"
+
+schedule ML_BSSN_constraints2_SelectBCs in ML_BSSN_constraints2_bc_group
+{
+  LANG: C
+  OPTIONS: level
+  SYNC: ML_cons_traceA
+  SYNC: ML_mom
+} "ML_BSSN_constraints2_SelectBCs"
+
+schedule group ApplyBCs as ML_BSSN_constraints2_ApplyBCs in ML_BSSN_constraints2_bc_group after ML_BSSN_constraints2_SelectBCs
+{
+  # no language specified
+} "Apply BCs for groups set in ML_BSSN_constraints2"
+
+schedule group ML_BSSN_constraints2_bc_group in ML_BSSN_constraints2_group after ML_BSSN_constraints2
+{
+  # no language specified
+} "ML_BSSN_constraints2"
+
+schedule group ML_BSSN_constraints2_bc_group at CCTK_POSTRESTRICT
+{
+  # no language specified
+} "ML_BSSN_constraints2"
+
+schedule group ML_BSSN_constraints2_bc_group at CCTK_POSTRESTRICTINITIAL
+{
+  # no language specified
+} "ML_BSSN_constraints2"
 
 schedule ML_BSSN_SelectBoundConds in MoL_PostStep
 {
diff --git a/ML_BSSN/src/ML_BSSN_Minkowski.c b/ML_BSSN/src/ML_BSSN_Minkowski.c
index 6700d1d..7564156 100644
--- a/ML_BSSN/src/ML_BSSN_Minkowski.c
+++ b/ML_BSSN/src/ML_BSSN_Minkowski.c
@@ -38,6 +38,9 @@ void ML_BSSN_Minkowski_Body(cGH const * restrict const cctkGH, int const dir, in
     return;
   }
   
+  const char *groups[] = {"ML_BSSN::ML_curv","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtshift","ML_BSSN::ML_Gamma","ML_BSSN::ML_lapse","ML_BSSN::ML_log_confac","ML_BSSN::ML_metric","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_Minkowski", 9, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
diff --git a/ML_BSSN/src/ML_BSSN_RHS1.c b/ML_BSSN/src/ML_BSSN_RHS1.c
index 22daf60..41adf25 100644
--- a/ML_BSSN/src/ML_BSSN_RHS1.c
+++ b/ML_BSSN/src/ML_BSSN_RHS1.c
@@ -71,6 +71,9 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     return;
   }
   
+  const char *groups[] = {"grid::coordinates","Grid::coordinates","ML_BSSN::ML_curv","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtlapserhs","ML_BSSN::ML_dtshift","ML_BSSN::ML_dtshiftrhs","ML_BSSN::ML_Gamma","ML_BSSN::ML_Gammarhs","ML_BSSN::ML_lapse","ML_BSSN::ML_lapserhs","ML_BSSN::ML_log_confac","ML_BSSN::ML_log_confacrhs","ML_BSSN::ML_metric","ML_BSSN::ML_metricrhs","ML_BSSN::ML_shift","ML_BSSN::ML_shiftrhs","ML_BSSN::ML_trace_curv","ML_BSSN::ML_trace_curvrhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_RHS1", 19, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
diff --git a/ML_BSSN/src/ML_BSSN_RHS2.c b/ML_BSSN/src/ML_BSSN_RHS2.c
index 672e384..60716d2 100644
--- a/ML_BSSN/src/ML_BSSN_RHS2.c
+++ b/ML_BSSN/src/ML_BSSN_RHS2.c
@@ -50,6 +50,9 @@ void ML_BSSN_RHS2_Body(cGH const * restrict const cctkGH, int const dir, int con
     return;
   }
   
+  const char *groups[] = {"ML_BSSN::ML_curv","ML_BSSN::ML_curvrhs","ML_BSSN::ML_Gamma","ML_BSSN::ML_lapse","ML_BSSN::ML_log_confac","ML_BSSN::ML_metric","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_RHS2", 8, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
diff --git a/ML_BSSN/src/ML_BSSN_RHSRadiativeBoundary.c b/ML_BSSN/src/ML_BSSN_RHSRadiativeBoundary.c
deleted file mode 100644
index 5572dae..0000000
--- a/ML_BSSN/src/ML_BSSN_RHSRadiativeBoundary.c
+++ /dev/null
@@ -1,331 +0,0 @@
-/*  File produced by Kranc */
-
-#define KRANC_C
-
-#include <assert.h>
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-#include "GenericFD.h"
-#include "Differencing.h"
-#include "loopcontrol.h"
-
-/* Define macros used in calculations */
-#define INITVALUE  (42)
-#define INV(x) ((1.0) / (x))
-#define SQR(x) ((x) * (x))
-#define CUB(x) ((x) * (x) * (x))
-#define QAD(x) ((x) * (x) * (x) * (x))
-
-void ML_BSSN_RHSRadiativeBoundary_SelectBCs(CCTK_ARGUMENTS)
-{
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
-  
-  CCTK_INT ierr = 0;
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_curvrhs","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_curvrhs.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_dtlapserhs","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_dtlapserhs.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_dtshiftrhs","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_dtshiftrhs.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_Gammarhs","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_Gammarhs.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_lapserhs","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_lapserhs.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_log_confacrhs","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_log_confacrhs.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_metricrhs","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_metricrhs.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_shiftrhs","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_shiftrhs.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_trace_curvrhs","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_trace_curvrhs.");
-  return;
-}
-
-void ML_BSSN_RHSRadiativeBoundary_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
-{
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
-  
-  
-  /* Declare finite differencing variables */
-  
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_RHSRadiativeBoundary_Body");
-  }
-  
-  if (cctk_iteration % ML_BSSN_RHSRadiativeBoundary_calc_every != ML_BSSN_RHSRadiativeBoundary_calc_offset)
-  {
-    return;
-  }
-  
-  /* Include user-supplied include files */
-  
-  /* Initialise finite differencing variables */
-  CCTK_REAL const dx = CCTK_DELTA_SPACE(0);
-  CCTK_REAL const dy = CCTK_DELTA_SPACE(1);
-  CCTK_REAL const dz = CCTK_DELTA_SPACE(2);
-  int const di = 1;
-  int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
-  int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
-  CCTK_REAL const dxi = 1.0 / dx;
-  CCTK_REAL const dyi = 1.0 / dy;
-  CCTK_REAL const dzi = 1.0 / dz;
-  CCTK_REAL const khalf = 0.5;
-  CCTK_REAL const kthird = 1/3.0;
-  CCTK_REAL const ktwothird = 2.0/3.0;
-  CCTK_REAL const kfourthird = 4.0/3.0;
-  CCTK_REAL const keightthird = 8.0/3.0;
-  CCTK_REAL const hdxi = 0.5 * dxi;
-  CCTK_REAL const hdyi = 0.5 * dyi;
-  CCTK_REAL const hdzi = 0.5 * dzi;
-  
-  /* Initialize predefined quantities */
-  CCTK_REAL const p1o12dx = INV(dx)/12.;
-  CCTK_REAL const p1o12dy = INV(dy)/12.;
-  CCTK_REAL const p1o12dz = INV(dz)/12.;
-  CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.;
-  CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.;
-  CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.;
-  CCTK_REAL const p1o24dx = INV(dx)/24.;
-  CCTK_REAL const p1o24dy = INV(dy)/24.;
-  CCTK_REAL const p1o24dz = INV(dz)/24.;
-  CCTK_REAL const p1o64dx = INV(dx)/64.;
-  CCTK_REAL const p1o64dy = INV(dy)/64.;
-  CCTK_REAL const p1o64dz = INV(dz)/64.;
-  CCTK_REAL const p1odx = INV(dx);
-  CCTK_REAL const p1ody = INV(dy);
-  CCTK_REAL const p1odz = INV(dz);
-  CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.;
-  CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.;
-  CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.;
-  
-  /* Loop over the grid points */
-  #pragma omp parallel
-  LC_LOOP3 (ML_BSSN_RHSRadiativeBoundary,
-            i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],
-            cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
-  {
-    // int index = INITVALUE;
-    int const index = CCTK_GFINDEX3D(cctkGH,i,j,k);
-    /* Declare derivatives */
-    
-    /* Assign local copies of grid functions */
-    CCTK_REAL  AL = A[index];
-    CCTK_REAL  alphaL = alpha[index];
-    CCTK_REAL  At11L = At11[index];
-    CCTK_REAL  At12L = At12[index];
-    CCTK_REAL  At13L = At13[index];
-    CCTK_REAL  At22L = At22[index];
-    CCTK_REAL  At23L = At23[index];
-    CCTK_REAL  At33L = At33[index];
-    CCTK_REAL  B1L = B1[index];
-    CCTK_REAL  B2L = B2[index];
-    CCTK_REAL  B3L = B3[index];
-    CCTK_REAL  beta1L = beta1[index];
-    CCTK_REAL  beta2L = beta2[index];
-    CCTK_REAL  beta3L = beta3[index];
-    CCTK_REAL  gt11L = gt11[index];
-    CCTK_REAL  gt12L = gt12[index];
-    CCTK_REAL  gt13L = gt13[index];
-    CCTK_REAL  gt22L = gt22[index];
-    CCTK_REAL  gt23L = gt23[index];
-    CCTK_REAL  gt33L = gt33[index];
-    CCTK_REAL  phiL = phi[index];
-    CCTK_REAL  trKL = trK[index];
-    CCTK_REAL  Xt1L = Xt1[index];
-    CCTK_REAL  Xt2L = Xt2[index];
-    CCTK_REAL  Xt3L = Xt3[index];
-    
-    /* Include user supplied include files */
-    
-    /* Precompute derivatives */
-    
-    /* Calculate temporaries and grid functions */
-    int dir1 = Sign(normal[0]);
-    
-    int dir2 = Sign(normal[1]);
-    
-    int dir3 = Sign(normal[2]);
-    
-    CCTK_REAL detgt = 1;
-    
-    CCTK_REAL gtu11 = INV(detgt)*(gt22L*gt33L - SQR(gt23L));
-    
-    CCTK_REAL gtu12 = (gt13L*gt23L - gt12L*gt33L)*INV(detgt);
-    
-    CCTK_REAL gtu13 = (-(gt13L*gt22L) + gt12L*gt23L)*INV(detgt);
-    
-    CCTK_REAL gtu22 = INV(detgt)*(gt11L*gt33L - SQR(gt13L));
-    
-    CCTK_REAL gtu23 = (gt12L*gt13L - gt11L*gt23L)*INV(detgt);
-    
-    CCTK_REAL gtu33 = INV(detgt)*(gt11L*gt22L - SQR(gt12L));
-    
-    CCTK_REAL em4phi = IfThen(conformalMethod,SQR(phiL),exp(-4*phiL));
-    
-    CCTK_REAL gu11 = em4phi*gtu11;
-    
-    CCTK_REAL gu12 = em4phi*gtu12;
-    
-    CCTK_REAL gu13 = em4phi*gtu13;
-    
-    CCTK_REAL gu22 = em4phi*gtu22;
-    
-    CCTK_REAL gu23 = em4phi*gtu23;
-    
-    CCTK_REAL gu33 = em4phi*gtu33;
-    
-    CCTK_REAL nn1 = normal[0];
-    
-    CCTK_REAL nn2 = normal[1];
-    
-    CCTK_REAL nn3 = normal[2];
-    
-    CCTK_REAL nu1 = gu11*nn1 + gu12*nn2 + gu13*nn3;
-    
-    CCTK_REAL nu2 = gu12*nn1 + gu22*nn2 + gu23*nn3;
-    
-    CCTK_REAL nu3 = gu13*nn1 + gu23*nn2 + gu33*nn3;
-    
-    CCTK_REAL nlen2 = nn1*nu1 + nn2*nu2 + nn3*nu3;
-    
-    CCTK_REAL nlen = pow(nlen2,0.5);
-    
-    CCTK_REAL su1 = nu1*INV(nlen);
-    
-    CCTK_REAL su2 = nu2*INV(nlen);
-    
-    CCTK_REAL su3 = nu3*INV(nlen);
-    
-    CCTK_REAL vg = pow(harmonicF,0.5);
-    
-    CCTK_REAL phirhsL = -((PDonesided1(phi, i, j, k)*su1 + 
-      PDonesided2(phi, i, j, k)*su2 + PDonesided3(phi, i, j, k)*su3)*vg);
-    
-    CCTK_REAL gt11rhsL = -(PDonesided1(gt11, i, j, k)*su1) - 
-      PDonesided2(gt11, i, j, k)*su2 - PDonesided3(gt11, i, j, k)*su3;
-    
-    CCTK_REAL gt12rhsL = -(PDonesided1(gt12, i, j, k)*su1) - 
-      PDonesided2(gt12, i, j, k)*su2 - PDonesided3(gt12, i, j, k)*su3;
-    
-    CCTK_REAL gt13rhsL = -(PDonesided1(gt13, i, j, k)*su1) - 
-      PDonesided2(gt13, i, j, k)*su2 - PDonesided3(gt13, i, j, k)*su3;
-    
-    CCTK_REAL gt22rhsL = -(PDonesided1(gt22, i, j, k)*su1) - 
-      PDonesided2(gt22, i, j, k)*su2 - PDonesided3(gt22, i, j, k)*su3;
-    
-    CCTK_REAL gt23rhsL = -(PDonesided1(gt23, i, j, k)*su1) - 
-      PDonesided2(gt23, i, j, k)*su2 - PDonesided3(gt23, i, j, k)*su3;
-    
-    CCTK_REAL gt33rhsL = -(PDonesided1(gt33, i, j, k)*su1) - 
-      PDonesided2(gt33, i, j, k)*su2 - PDonesided3(gt33, i, j, k)*su3;
-    
-    CCTK_REAL trKrhsL = -((PDonesided1(trK, i, j, k)*su1 + 
-      PDonesided2(trK, i, j, k)*su2 + PDonesided3(trK, i, j, k)*su3)*vg);
-    
-    CCTK_REAL At11rhsL = -(PDonesided1(At11, i, j, k)*su1) - 
-      PDonesided2(At11, i, j, k)*su2 - PDonesided3(At11, i, j, k)*su3;
-    
-    CCTK_REAL At12rhsL = -(PDonesided1(At12, i, j, k)*su1) - 
-      PDonesided2(At12, i, j, k)*su2 - PDonesided3(At12, i, j, k)*su3;
-    
-    CCTK_REAL At13rhsL = -(PDonesided1(At13, i, j, k)*su1) - 
-      PDonesided2(At13, i, j, k)*su2 - PDonesided3(At13, i, j, k)*su3;
-    
-    CCTK_REAL At22rhsL = -(PDonesided1(At22, i, j, k)*su1) - 
-      PDonesided2(At22, i, j, k)*su2 - PDonesided3(At22, i, j, k)*su3;
-    
-    CCTK_REAL At23rhsL = -(PDonesided1(At23, i, j, k)*su1) - 
-      PDonesided2(At23, i, j, k)*su2 - PDonesided3(At23, i, j, k)*su3;
-    
-    CCTK_REAL At33rhsL = -(PDonesided1(At33, i, j, k)*su1) - 
-      PDonesided2(At33, i, j, k)*su2 - PDonesided3(At33, i, j, k)*su3;
-    
-    CCTK_REAL Xt1rhsL = -(PDonesided1(Xt1, i, j, k)*su1) - 
-      PDonesided2(Xt1, i, j, k)*su2 - PDonesided3(Xt1, i, j, k)*su3;
-    
-    CCTK_REAL Xt2rhsL = -(PDonesided1(Xt2, i, j, k)*su1) - 
-      PDonesided2(Xt2, i, j, k)*su2 - PDonesided3(Xt2, i, j, k)*su3;
-    
-    CCTK_REAL Xt3rhsL = -(PDonesided1(Xt3, i, j, k)*su1) - 
-      PDonesided2(Xt3, i, j, k)*su2 - PDonesided3(Xt3, i, j, k)*su3;
-    
-    CCTK_REAL alpharhsL = -((PDonesided1(alpha, i, j, k)*su1 + 
-      PDonesided2(alpha, i, j, k)*su2 + PDonesided3(alpha, i, j, 
-      k)*su3)*vg);
-    
-    CCTK_REAL ArhsL = -((PDonesided1(A, i, j, k)*su1 + PDonesided2(A, 
-      i, j, k)*su2 + PDonesided3(A, i, j, k)*su3)*vg);
-    
-    CCTK_REAL beta1rhsL = -(PDonesided1(beta1, i, j, k)*su1) - 
-      PDonesided2(beta1, i, j, k)*su2 - PDonesided3(beta1, i, j, k)*su3;
-    
-    CCTK_REAL beta2rhsL = -(PDonesided1(beta2, i, j, k)*su1) - 
-      PDonesided2(beta2, i, j, k)*su2 - PDonesided3(beta2, i, j, k)*su3;
-    
-    CCTK_REAL beta3rhsL = -(PDonesided1(beta3, i, j, k)*su1) - 
-      PDonesided2(beta3, i, j, k)*su2 - PDonesided3(beta3, i, j, k)*su3;
-    
-    CCTK_REAL B1rhsL = -(PDonesided1(B1, i, j, k)*su1) - 
-      PDonesided2(B1, i, j, k)*su2 - PDonesided3(B1, i, j, k)*su3;
-    
-    CCTK_REAL B2rhsL = -(PDonesided1(B2, i, j, k)*su1) - 
-      PDonesided2(B2, i, j, k)*su2 - PDonesided3(B2, i, j, k)*su3;
-    
-    CCTK_REAL B3rhsL = -(PDonesided1(B3, i, j, k)*su1) - 
-      PDonesided2(B3, i, j, k)*su2 - PDonesided3(B3, i, j, k)*su3;
-    
-    
-    /* Copy local copies back to grid functions */
-    alpharhs[index] = alpharhsL;
-    Arhs[index] = ArhsL;
-    At11rhs[index] = At11rhsL;
-    At12rhs[index] = At12rhsL;
-    At13rhs[index] = At13rhsL;
-    At22rhs[index] = At22rhsL;
-    At23rhs[index] = At23rhsL;
-    At33rhs[index] = At33rhsL;
-    B1rhs[index] = B1rhsL;
-    B2rhs[index] = B2rhsL;
-    B3rhs[index] = B3rhsL;
-    beta1rhs[index] = beta1rhsL;
-    beta2rhs[index] = beta2rhsL;
-    beta3rhs[index] = beta3rhsL;
-    gt11rhs[index] = gt11rhsL;
-    gt12rhs[index] = gt12rhsL;
-    gt13rhs[index] = gt13rhsL;
-    gt22rhs[index] = gt22rhsL;
-    gt23rhs[index] = gt23rhsL;
-    gt33rhs[index] = gt33rhsL;
-    phirhs[index] = phirhsL;
-    trKrhs[index] = trKrhsL;
-    Xt1rhs[index] = Xt1rhsL;
-    Xt2rhs[index] = Xt2rhsL;
-    Xt3rhs[index] = Xt3rhsL;
-  }
-  LC_ENDLOOP3 (ML_BSSN_RHSRadiativeBoundary);
-}
-
-void ML_BSSN_RHSRadiativeBoundary(CCTK_ARGUMENTS)
-{
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
-  
-  GenericFD_LoopOverBoundary(cctkGH, &ML_BSSN_RHSRadiativeBoundary_Body);
-}
diff --git a/ML_BSSN/src/ML_BSSN_RHSStaticBoundary.c b/ML_BSSN/src/ML_BSSN_RHSStaticBoundary.c
index a322308..7f177c5 100644
--- a/ML_BSSN/src/ML_BSSN_RHSStaticBoundary.c
+++ b/ML_BSSN/src/ML_BSSN_RHSStaticBoundary.c
@@ -74,6 +74,9 @@ void ML_BSSN_RHSStaticBoundary_Body(cGH const * restrict const cctkGH, int const
     return;
   }
   
+  const char *groups[] = {"ML_BSSN::ML_curvrhs","ML_BSSN::ML_dtlapserhs","ML_BSSN::ML_dtshiftrhs","ML_BSSN::ML_Gammarhs","ML_BSSN::ML_lapserhs","ML_BSSN::ML_log_confacrhs","ML_BSSN::ML_metricrhs","ML_BSSN::ML_shiftrhs","ML_BSSN::ML_trace_curvrhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_RHSStaticBoundary", 9, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
diff --git a/ML_BSSN/src/ML_BSSN_boundary.c b/ML_BSSN/src/ML_BSSN_boundary.c
index eac45fc..0abd9b1 100644
--- a/ML_BSSN/src/ML_BSSN_boundary.c
+++ b/ML_BSSN/src/ML_BSSN_boundary.c
@@ -74,6 +74,9 @@ void ML_BSSN_boundary_Body(cGH const * restrict const cctkGH, int const dir, int
     return;
   }
   
+  const char *groups[] = {"ML_BSSN::ML_curv","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtshift","ML_BSSN::ML_Gamma","ML_BSSN::ML_lapse","ML_BSSN::ML_log_confac","ML_BSSN::ML_metric","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_boundary", 9, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
diff --git a/ML_BSSN/src/ML_BSSN_constraints.c b/ML_BSSN/src/ML_BSSN_constraints1.c
index eb5642a..9f344f3 100644
--- a/ML_BSSN/src/ML_BSSN_constraints.c
+++ b/ML_BSSN/src/ML_BSSN_constraints1.c
@@ -20,7 +20,7 @@
 #define CUB(x) ((x) * (x) * (x))
 #define QAD(x) ((x) * (x) * (x) * (x))
 
-void ML_BSSN_constraints_SelectBCs(CCTK_ARGUMENTS)
+void ML_BSSN_constraints1_SelectBCs(CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
@@ -32,19 +32,13 @@ void ML_BSSN_constraints_SelectBCs(CCTK_ARGUMENTS)
   ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_cons_Gamma","flat");
   if (ierr < 0)
     CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_cons_Gamma.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_cons_traceA","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_cons_traceA.");
   ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_Ham","flat");
   if (ierr < 0)
     CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_Ham.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_mom","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_mom.");
   return;
 }
 
-void ML_BSSN_constraints_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+void ML_BSSN_constraints1_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
@@ -54,14 +48,17 @@ void ML_BSSN_constraints_Body(cGH const * restrict const cctkGH, int const dir,
   
   if (verbose > 1)
   {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_constraints_Body");
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_constraints1_Body");
   }
   
-  if (cctk_iteration % ML_BSSN_constraints_calc_every != ML_BSSN_constraints_calc_offset)
+  if (cctk_iteration % ML_BSSN_constraints1_calc_every != ML_BSSN_constraints1_calc_offset)
   {
     return;
   }
   
+  const char *groups[] = {"ML_BSSN::ML_cons_detg","ML_BSSN::ML_cons_Gamma","ML_BSSN::ML_curv","ML_BSSN::ML_Gamma","ML_BSSN::ML_Ham","ML_BSSN::ML_lapse","ML_BSSN::ML_log_confac","ML_BSSN::ML_metric","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_constraints1", 10, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
@@ -105,31 +102,13 @@ void ML_BSSN_constraints_Body(cGH const * restrict const cctkGH, int const dir,
   
   /* Loop over the grid points */
   #pragma omp parallel
-  LC_LOOP3 (ML_BSSN_constraints,
+  LC_LOOP3 (ML_BSSN_constraints1,
             i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],
             cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
     // int index = INITVALUE;
     int const index = CCTK_GFINDEX3D(cctkGH,i,j,k);
     /* Declare derivatives */
-    // CCTK_REAL PDstandardNth1At11 = INITVALUE;
-    // CCTK_REAL PDstandardNth2At11 = INITVALUE;
-    // CCTK_REAL PDstandardNth3At11 = INITVALUE;
-    // CCTK_REAL PDstandardNth1At12 = INITVALUE;
-    // CCTK_REAL PDstandardNth2At12 = INITVALUE;
-    // CCTK_REAL PDstandardNth3At12 = INITVALUE;
-    // CCTK_REAL PDstandardNth1At13 = INITVALUE;
-    // CCTK_REAL PDstandardNth2At13 = INITVALUE;
-    // CCTK_REAL PDstandardNth3At13 = INITVALUE;
-    // CCTK_REAL PDstandardNth1At22 = INITVALUE;
-    // CCTK_REAL PDstandardNth2At22 = INITVALUE;
-    // CCTK_REAL PDstandardNth3At22 = INITVALUE;
-    // CCTK_REAL PDstandardNth1At23 = INITVALUE;
-    // CCTK_REAL PDstandardNth2At23 = INITVALUE;
-    // CCTK_REAL PDstandardNth3At23 = INITVALUE;
-    // CCTK_REAL PDstandardNth1At33 = INITVALUE;
-    // CCTK_REAL PDstandardNth2At33 = INITVALUE;
-    // CCTK_REAL PDstandardNth3At33 = INITVALUE;
     // CCTK_REAL PDstandardNth1gt11 = INITVALUE;
     // CCTK_REAL PDstandardNth2gt11 = INITVALUE;
     // CCTK_REAL PDstandardNth3gt11 = INITVALUE;
@@ -193,9 +172,6 @@ void ML_BSSN_constraints_Body(cGH const * restrict const cctkGH, int const dir,
     // CCTK_REAL PDstandardNth12phi = INITVALUE;
     // CCTK_REAL PDstandardNth13phi = INITVALUE;
     // CCTK_REAL PDstandardNth23phi = INITVALUE;
-    // CCTK_REAL PDstandardNth1trK = INITVALUE;
-    // CCTK_REAL PDstandardNth2trK = INITVALUE;
-    // CCTK_REAL PDstandardNth3trK = INITVALUE;
     // CCTK_REAL PDstandardNth1Xt1 = INITVALUE;
     // CCTK_REAL PDstandardNth2Xt1 = INITVALUE;
     // CCTK_REAL PDstandardNth3Xt1 = INITVALUE;
@@ -242,24 +218,6 @@ void ML_BSSN_constraints_Body(cGH const * restrict const cctkGH, int const dir,
     /* Include user supplied include files */
     
     /* Precompute derivatives */
-    CCTK_REAL const PDstandardNth1At11 = PDstandardNth1(At11, i, j, k);
-    CCTK_REAL const PDstandardNth2At11 = PDstandardNth2(At11, i, j, k);
-    CCTK_REAL const PDstandardNth3At11 = PDstandardNth3(At11, i, j, k);
-    CCTK_REAL const PDstandardNth1At12 = PDstandardNth1(At12, i, j, k);
-    CCTK_REAL const PDstandardNth2At12 = PDstandardNth2(At12, i, j, k);
-    CCTK_REAL const PDstandardNth3At12 = PDstandardNth3(At12, i, j, k);
-    CCTK_REAL const PDstandardNth1At13 = PDstandardNth1(At13, i, j, k);
-    CCTK_REAL const PDstandardNth2At13 = PDstandardNth2(At13, i, j, k);
-    CCTK_REAL const PDstandardNth3At13 = PDstandardNth3(At13, i, j, k);
-    CCTK_REAL const PDstandardNth1At22 = PDstandardNth1(At22, i, j, k);
-    CCTK_REAL const PDstandardNth2At22 = PDstandardNth2(At22, i, j, k);
-    CCTK_REAL const PDstandardNth3At22 = PDstandardNth3(At22, i, j, k);
-    CCTK_REAL const PDstandardNth1At23 = PDstandardNth1(At23, i, j, k);
-    CCTK_REAL const PDstandardNth2At23 = PDstandardNth2(At23, i, j, k);
-    CCTK_REAL const PDstandardNth3At23 = PDstandardNth3(At23, i, j, k);
-    CCTK_REAL const PDstandardNth1At33 = PDstandardNth1(At33, i, j, k);
-    CCTK_REAL const PDstandardNth2At33 = PDstandardNth2(At33, i, j, k);
-    CCTK_REAL const PDstandardNth3At33 = PDstandardNth3(At33, i, j, k);
     CCTK_REAL const PDstandardNth1gt11 = PDstandardNth1(gt11, i, j, k);
     CCTK_REAL const PDstandardNth2gt11 = PDstandardNth2(gt11, i, j, k);
     CCTK_REAL const PDstandardNth3gt11 = PDstandardNth3(gt11, i, j, k);
@@ -323,9 +281,6 @@ void ML_BSSN_constraints_Body(cGH const * restrict const cctkGH, int const dir,
     CCTK_REAL const PDstandardNth12phi = PDstandardNth12(phi, i, j, k);
     CCTK_REAL const PDstandardNth13phi = PDstandardNth13(phi, i, j, k);
     CCTK_REAL const PDstandardNth23phi = PDstandardNth23(phi, i, j, k);
-    CCTK_REAL const PDstandardNth1trK = PDstandardNth1(trK, i, j, k);
-    CCTK_REAL const PDstandardNth2trK = PDstandardNth2(trK, i, j, k);
-    CCTK_REAL const PDstandardNth3trK = PDstandardNth3(trK, i, j, k);
     CCTK_REAL const PDstandardNth1Xt1 = PDstandardNth1(Xt1, i, j, k);
     CCTK_REAL const PDstandardNth2Xt1 = PDstandardNth2(Xt1, i, j, k);
     CCTK_REAL const PDstandardNth3Xt1 = PDstandardNth3(Xt1, i, j, k);
@@ -927,75 +882,11 @@ void ML_BSSN_constraints_Body(cGH const * restrict const cctkGH, int const dir,
       beta2L*beta3L*eTyzL) + eTxxL*SQR(beta1L) + eTyyL*SQR(beta2L) + 
       eTzzL*SQR(beta3L));
     
-    CCTK_REAL S1 = (-eTtxL + beta1L*eTxxL + beta2L*eTxyL + 
-      beta3L*eTxzL)*INV(alphaL);
-    
-    CCTK_REAL S2 = (-eTtyL + beta1L*eTxyL + beta2L*eTyyL + 
-      beta3L*eTyzL)*INV(alphaL);
-    
-    CCTK_REAL S3 = (-eTtzL + beta1L*eTxzL + beta2L*eTyzL + 
-      beta3L*eTzzL)*INV(alphaL);
-    
     CCTK_REAL HL = -2.*(Atm12*Atm21 + Atm13*Atm31 + Atm23*Atm32) - 
       50.26548245743669181540229413247204614715*rho + trR - 1.*(SQR(Atm11) + 
       SQR(Atm22) + SQR(Atm33)) + 
       0.6666666666666666666666666666666666666667*SQR(trKL);
     
-    CCTK_REAL M1L = -2.*((At12L*Gt211 + At13L*Gt311)*gtu11 + 
-      At11L*Gt123*gtu23) + At12L*((6.*cdphi1 - 1.*Gt111)*gtu12 - 
-      3.*Gt213*gtu13 + 6.*(cdphi2*gtu22 + cdphi3*gtu23) - 1.*Gt233*gtu33) - 
-      1.*((At22L*Gt212 + At12L*(Gt112 + Gt222) + At23L*Gt312)*gtu22 + 
-      (At13L*Gt112 + At12L*Gt113 + At23L*Gt212)*gtu23 + (At13L*Gt113 + 
-      At23L*Gt213 + At33L*Gt313)*gtu33) + At11L*((6.*cdphi1 - 2.*Gt111)*gtu11 
-      + 6.*(cdphi2*gtu12 + cdphi3*gtu13) - 1.*(Gt122*gtu22 + Gt133*gtu33)) + 
-      At13L*((6.*cdphi1 - 1.*Gt111)*gtu13 + 6.*(cdphi2*gtu23 + cdphi3*gtu33) 
-      - 1.*(Gt322*gtu22 + Gt333*gtu33)) + gtu11*PDstandardNth1At11 - 
-      0.6666666666666666666666666666666666666667*PDstandardNth1trK + 
-      gtu12*(-1.*(At22L*Gt211 + At23L*Gt311) - 3.*(At11L*Gt112 + At12L*Gt212 
-      + At13L*Gt312) + PDstandardNth1At12 + PDstandardNth2At11) + 
-      gtu22*PDstandardNth2At12 + gtu13*(-1.*(At23L*Gt211 + At33L*Gt311) - 
-      3.*(At11L*Gt113 + At13L*Gt313) + PDstandardNth1At13 + 
-      PDstandardNth3At11) + gtu23*(-1.*(At22L*Gt213 + At33L*Gt312 + 
-      At23L*Gt313) - 2.*(At12L*Gt223 + At13L*Gt323) + PDstandardNth2At13 + 
-      PDstandardNth3At12) + gtu33*PDstandardNth3At13 - 
-      25.13274122871834590770114706623602307358*S1;
-    
-    CCTK_REAL M2L = At12L*((6.*cdphi1 - 1.*Gt111)*gtu11 + 6.*(cdphi2*gtu12 
-      + cdphi3*gtu13) - 2.*Gt122*gtu22 - 3.*Gt123*gtu23 - 1.*Gt133*gtu33) + 
-      At22L*((6.*cdphi2 - 2.*Gt222)*gtu22 + 6.*cdphi3*gtu23 - 1.*Gt233*gtu33) 
-      + At23L*(-2.*Gt322*gtu22 - 1.*Gt333*gtu33 + 6.*(cdphi2*gtu23 + 
-      cdphi3*gtu33)) - 1.*((At11L*Gt112 + At22L*Gt211 + At12L*Gt212 + 
-      At23L*Gt311 + At13L*Gt312)*gtu11 + Gt122*(At11L*gtu12 + At13L*gtu23) + 
-      (At23L*Gt223 + At33L*Gt323)*gtu33 + At13L*(Gt112*gtu13 + Gt123*gtu33)) 
-      + gtu11*PDstandardNth1At12 + gtu12*(At22L*(6.*cdphi1 - 3.*Gt212) + 
-      At12L*(-3.*Gt112 - 1.*Gt222) - 3.*At23L*Gt312 - 1.*At13L*Gt322 + 
-      PDstandardNth1At22 + PDstandardNth2At12) + gtu22*PDstandardNth2At22 - 
-      0.6666666666666666666666666666666666666667*PDstandardNth2trK + 
-      gtu13*(-2.*(At12L*Gt113 + At22L*Gt213) + At23L*(6.*cdphi1 - 1.*Gt212 - 
-      2.*Gt313) - 1.*(At11L*Gt123 + At12L*Gt223 + At33L*Gt312 + At13L*Gt323) 
-      + PDstandardNth1At23 + PDstandardNth3At12) + gtu23*(-1.*(At23L*Gt222 + 
-      At33L*Gt322) - 3.*(At22L*Gt223 + At23L*Gt323) + PDstandardNth2At23 + 
-      PDstandardNth3At22) + gtu33*PDstandardNth3At23 - 
-      25.13274122871834590770114706623602307358*S2;
-    
-    CCTK_REAL M3L = -1.*((At11L*Gt113 + At23L*Gt211 + At12L*Gt213 + 
-      At33L*Gt311)*gtu11 + (At22L*Gt223 + At33L*Gt322 + At23L*Gt323)*gtu22 + 
-      At12L*(Gt113*gtu12 + Gt123*gtu22) + Gt133*(At11L*gtu13 + At12L*gtu23)) 
-      + At13L*((6.*cdphi1 - 1.*(Gt111 + Gt313))*gtu11 + 6.*(cdphi2*gtu12 + 
-      cdphi3*gtu13) - 1.*Gt122*gtu22 - 3.*Gt123*gtu23 - 2.*Gt133*gtu33) + 
-      At23L*((6.*cdphi2 - 1.*Gt222)*gtu22 + 6.*cdphi3*gtu23 - 2.*Gt233*gtu33) 
-      + gtu11*PDstandardNth1At13 + gtu12*(-2.*(At13L*Gt112 + At33L*Gt312) + 
-      At23L*(6.*cdphi1 - 2.*Gt212 - 1.*Gt313) - 1.*(At11L*Gt123 + At22L*Gt213 
-      + At12L*Gt223 + At13L*Gt323) + PDstandardNth1At23 + PDstandardNth2At13) 
-      + gtu22*PDstandardNth2At23 + gtu13*(-3.*(At13L*Gt113 + At23L*Gt213) + 
-      At33L*(6.*cdphi1 - 3.*Gt313) - 1.*(At12L*Gt233 + At13L*Gt333) + 
-      PDstandardNth1At33 + PDstandardNth3At13) + gtu23*(-1.*At22L*Gt233 + 
-      At33L*(6.*cdphi2 - 3.*Gt323) + At23L*(-3.*Gt223 - 1.*Gt333) + 
-      PDstandardNth2At33 + PDstandardNth3At23) + gtu33*(At33L*(6.*cdphi3 - 
-      2.*Gt333) + PDstandardNth3At33) - 
-      0.6666666666666666666666666666666666666667*PDstandardNth3trK - 
-      25.13274122871834590770114706623602307358*S3;
-    
     CCTK_REAL cSL = Log(detgt);
     
     CCTK_REAL cXt1L = Gt111*gtu11 + Gt122*gtu22 + 2*(Gt112*gtu12 + 
@@ -1007,28 +898,21 @@ void ML_BSSN_constraints_Body(cGH const * restrict const cctkGH, int const dir,
     CCTK_REAL cXt3L = Gt311*gtu11 + Gt322*gtu22 + 2*(Gt312*gtu12 + 
       Gt313*gtu13 + Gt323*gtu23) + Gt333*gtu33 - Xt3L;
     
-    CCTK_REAL cAL = At11L*gtu11 + At22L*gtu22 + 2*(At12L*gtu12 + 
-      At13L*gtu13 + At23L*gtu23) + At33L*gtu33;
-    
     
     /* Copy local copies back to grid functions */
-    cA[index] = cAL;
     cS[index] = cSL;
     cXt1[index] = cXt1L;
     cXt2[index] = cXt2L;
     cXt3[index] = cXt3L;
     H[index] = HL;
-    M1[index] = M1L;
-    M2[index] = M2L;
-    M3[index] = M3L;
   }
-  LC_ENDLOOP3 (ML_BSSN_constraints);
+  LC_ENDLOOP3 (ML_BSSN_constraints1);
 }
 
-void ML_BSSN_constraints(CCTK_ARGUMENTS)
+void ML_BSSN_constraints1(CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
-  GenericFD_LoopOverInterior(cctkGH, &ML_BSSN_constraints_Body);
+  GenericFD_LoopOverInterior(cctkGH, &ML_BSSN_constraints1_Body);
 }
diff --git a/ML_BSSN/src/ML_BSSN_constraints2.c b/ML_BSSN/src/ML_BSSN_constraints2.c
new file mode 100644
index 0000000..dd6174c
--- /dev/null
+++ b/ML_BSSN/src/ML_BSSN_constraints2.c
@@ -0,0 +1,405 @@
+/*  File produced by Kranc */
+
+#define KRANC_C
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "GenericFD.h"
+#include "Differencing.h"
+#include "loopcontrol.h"
+
+/* Define macros used in calculations */
+#define INITVALUE  (42)
+#define INV(x) ((1.0) / (x))
+#define SQR(x) ((x) * (x))
+#define CUB(x) ((x) * (x) * (x))
+#define QAD(x) ((x) * (x) * (x) * (x))
+
+void ML_BSSN_constraints2_SelectBCs(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_INT ierr = 0;
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_cons_traceA","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_cons_traceA.");
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_mom","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_mom.");
+  return;
+}
+
+void ML_BSSN_constraints2_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  
+  /* Declare finite differencing variables */
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_constraints2_Body");
+  }
+  
+  if (cctk_iteration % ML_BSSN_constraints2_calc_every != ML_BSSN_constraints2_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"ML_BSSN::ML_cons_traceA","ML_BSSN::ML_curv","ML_BSSN::ML_lapse","ML_BSSN::ML_log_confac","ML_BSSN::ML_metric","ML_BSSN::ML_mom","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_constraints2", 8, groups);
+  
+  /* Include user-supplied include files */
+  
+  /* Initialise finite differencing variables */
+  CCTK_REAL const dx = CCTK_DELTA_SPACE(0);
+  CCTK_REAL const dy = CCTK_DELTA_SPACE(1);
+  CCTK_REAL const dz = CCTK_DELTA_SPACE(2);
+  int const di = 1;
+  int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  CCTK_REAL const dxi = 1.0 / dx;
+  CCTK_REAL const dyi = 1.0 / dy;
+  CCTK_REAL const dzi = 1.0 / dz;
+  CCTK_REAL const khalf = 0.5;
+  CCTK_REAL const kthird = 1/3.0;
+  CCTK_REAL const ktwothird = 2.0/3.0;
+  CCTK_REAL const kfourthird = 4.0/3.0;
+  CCTK_REAL const keightthird = 8.0/3.0;
+  CCTK_REAL const hdxi = 0.5 * dxi;
+  CCTK_REAL const hdyi = 0.5 * dyi;
+  CCTK_REAL const hdzi = 0.5 * dzi;
+  
+  /* Initialize predefined quantities */
+  CCTK_REAL const p1o12dx = INV(dx)/12.;
+  CCTK_REAL const p1o12dy = INV(dy)/12.;
+  CCTK_REAL const p1o12dz = INV(dz)/12.;
+  CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.;
+  CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.;
+  CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.;
+  CCTK_REAL const p1o24dx = INV(dx)/24.;
+  CCTK_REAL const p1o24dy = INV(dy)/24.;
+  CCTK_REAL const p1o24dz = INV(dz)/24.;
+  CCTK_REAL const p1o64dx = INV(dx)/64.;
+  CCTK_REAL const p1o64dy = INV(dy)/64.;
+  CCTK_REAL const p1o64dz = INV(dz)/64.;
+  CCTK_REAL const p1odx = INV(dx);
+  CCTK_REAL const p1ody = INV(dy);
+  CCTK_REAL const p1odz = INV(dz);
+  CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.;
+  CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.;
+  CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.;
+  
+  /* Loop over the grid points */
+  #pragma omp parallel
+  LC_LOOP3 (ML_BSSN_constraints2,
+            i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],
+            cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+  {
+    // int index = INITVALUE;
+    int const index = CCTK_GFINDEX3D(cctkGH,i,j,k);
+    /* Declare derivatives */
+    // CCTK_REAL PDstandardNth1At11 = INITVALUE;
+    // CCTK_REAL PDstandardNth2At11 = INITVALUE;
+    // CCTK_REAL PDstandardNth3At11 = INITVALUE;
+    // CCTK_REAL PDstandardNth1At12 = INITVALUE;
+    // CCTK_REAL PDstandardNth2At12 = INITVALUE;
+    // CCTK_REAL PDstandardNth3At12 = INITVALUE;
+    // CCTK_REAL PDstandardNth1At13 = INITVALUE;
+    // CCTK_REAL PDstandardNth2At13 = INITVALUE;
+    // CCTK_REAL PDstandardNth3At13 = INITVALUE;
+    // CCTK_REAL PDstandardNth1At22 = INITVALUE;
+    // CCTK_REAL PDstandardNth2At22 = INITVALUE;
+    // CCTK_REAL PDstandardNth3At22 = INITVALUE;
+    // CCTK_REAL PDstandardNth1At23 = INITVALUE;
+    // CCTK_REAL PDstandardNth2At23 = INITVALUE;
+    // CCTK_REAL PDstandardNth3At23 = INITVALUE;
+    // CCTK_REAL PDstandardNth1At33 = INITVALUE;
+    // CCTK_REAL PDstandardNth2At33 = INITVALUE;
+    // CCTK_REAL PDstandardNth3At33 = INITVALUE;
+    // CCTK_REAL PDstandardNth1gt11 = INITVALUE;
+    // CCTK_REAL PDstandardNth2gt11 = INITVALUE;
+    // CCTK_REAL PDstandardNth3gt11 = INITVALUE;
+    // CCTK_REAL PDstandardNth1gt12 = INITVALUE;
+    // CCTK_REAL PDstandardNth2gt12 = INITVALUE;
+    // CCTK_REAL PDstandardNth3gt12 = INITVALUE;
+    // CCTK_REAL PDstandardNth1gt13 = INITVALUE;
+    // CCTK_REAL PDstandardNth2gt13 = INITVALUE;
+    // CCTK_REAL PDstandardNth3gt13 = INITVALUE;
+    // CCTK_REAL PDstandardNth1gt22 = INITVALUE;
+    // CCTK_REAL PDstandardNth2gt22 = INITVALUE;
+    // CCTK_REAL PDstandardNth3gt22 = INITVALUE;
+    // CCTK_REAL PDstandardNth1gt23 = INITVALUE;
+    // CCTK_REAL PDstandardNth2gt23 = INITVALUE;
+    // CCTK_REAL PDstandardNth3gt23 = INITVALUE;
+    // CCTK_REAL PDstandardNth1gt33 = INITVALUE;
+    // CCTK_REAL PDstandardNth2gt33 = INITVALUE;
+    // CCTK_REAL PDstandardNth3gt33 = INITVALUE;
+    // CCTK_REAL PDstandardNth1phi = INITVALUE;
+    // CCTK_REAL PDstandardNth2phi = INITVALUE;
+    // CCTK_REAL PDstandardNth3phi = INITVALUE;
+    // CCTK_REAL PDstandardNth1trK = INITVALUE;
+    // CCTK_REAL PDstandardNth2trK = INITVALUE;
+    // CCTK_REAL PDstandardNth3trK = INITVALUE;
+    
+    /* Assign local copies of grid functions */
+    CCTK_REAL  alphaL = alpha[index];
+    CCTK_REAL  At11L = At11[index];
+    CCTK_REAL  At12L = At12[index];
+    CCTK_REAL  At13L = At13[index];
+    CCTK_REAL  At22L = At22[index];
+    CCTK_REAL  At23L = At23[index];
+    CCTK_REAL  At33L = At33[index];
+    CCTK_REAL  beta1L = beta1[index];
+    CCTK_REAL  beta2L = beta2[index];
+    CCTK_REAL  beta3L = beta3[index];
+    CCTK_REAL  eTtxL = (*stress_energy_state) ? (eTtx[index]) : 0.0;
+    CCTK_REAL  eTtyL = (*stress_energy_state) ? (eTty[index]) : 0.0;
+    CCTK_REAL  eTtzL = (*stress_energy_state) ? (eTtz[index]) : 0.0;
+    CCTK_REAL  eTxxL = (*stress_energy_state) ? (eTxx[index]) : 0.0;
+    CCTK_REAL  eTxyL = (*stress_energy_state) ? (eTxy[index]) : 0.0;
+    CCTK_REAL  eTxzL = (*stress_energy_state) ? (eTxz[index]) : 0.0;
+    CCTK_REAL  eTyyL = (*stress_energy_state) ? (eTyy[index]) : 0.0;
+    CCTK_REAL  eTyzL = (*stress_energy_state) ? (eTyz[index]) : 0.0;
+    CCTK_REAL  eTzzL = (*stress_energy_state) ? (eTzz[index]) : 0.0;
+    CCTK_REAL  gt11L = gt11[index];
+    CCTK_REAL  gt12L = gt12[index];
+    CCTK_REAL  gt13L = gt13[index];
+    CCTK_REAL  gt22L = gt22[index];
+    CCTK_REAL  gt23L = gt23[index];
+    CCTK_REAL  gt33L = gt33[index];
+    CCTK_REAL  phiL = phi[index];
+    CCTK_REAL  trKL = trK[index];
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDstandardNth1At11 = PDstandardNth1(At11, i, j, k);
+    CCTK_REAL const PDstandardNth2At11 = PDstandardNth2(At11, i, j, k);
+    CCTK_REAL const PDstandardNth3At11 = PDstandardNth3(At11, i, j, k);
+    CCTK_REAL const PDstandardNth1At12 = PDstandardNth1(At12, i, j, k);
+    CCTK_REAL const PDstandardNth2At12 = PDstandardNth2(At12, i, j, k);
+    CCTK_REAL const PDstandardNth3At12 = PDstandardNth3(At12, i, j, k);
+    CCTK_REAL const PDstandardNth1At13 = PDstandardNth1(At13, i, j, k);
+    CCTK_REAL const PDstandardNth2At13 = PDstandardNth2(At13, i, j, k);
+    CCTK_REAL const PDstandardNth3At13 = PDstandardNth3(At13, i, j, k);
+    CCTK_REAL const PDstandardNth1At22 = PDstandardNth1(At22, i, j, k);
+    CCTK_REAL const PDstandardNth2At22 = PDstandardNth2(At22, i, j, k);
+    CCTK_REAL const PDstandardNth3At22 = PDstandardNth3(At22, i, j, k);
+    CCTK_REAL const PDstandardNth1At23 = PDstandardNth1(At23, i, j, k);
+    CCTK_REAL const PDstandardNth2At23 = PDstandardNth2(At23, i, j, k);
+    CCTK_REAL const PDstandardNth3At23 = PDstandardNth3(At23, i, j, k);
+    CCTK_REAL const PDstandardNth1At33 = PDstandardNth1(At33, i, j, k);
+    CCTK_REAL const PDstandardNth2At33 = PDstandardNth2(At33, i, j, k);
+    CCTK_REAL const PDstandardNth3At33 = PDstandardNth3(At33, i, j, k);
+    CCTK_REAL const PDstandardNth1gt11 = PDstandardNth1(gt11, i, j, k);
+    CCTK_REAL const PDstandardNth2gt11 = PDstandardNth2(gt11, i, j, k);
+    CCTK_REAL const PDstandardNth3gt11 = PDstandardNth3(gt11, i, j, k);
+    CCTK_REAL const PDstandardNth1gt12 = PDstandardNth1(gt12, i, j, k);
+    CCTK_REAL const PDstandardNth2gt12 = PDstandardNth2(gt12, i, j, k);
+    CCTK_REAL const PDstandardNth3gt12 = PDstandardNth3(gt12, i, j, k);
+    CCTK_REAL const PDstandardNth1gt13 = PDstandardNth1(gt13, i, j, k);
+    CCTK_REAL const PDstandardNth2gt13 = PDstandardNth2(gt13, i, j, k);
+    CCTK_REAL const PDstandardNth3gt13 = PDstandardNth3(gt13, i, j, k);
+    CCTK_REAL const PDstandardNth1gt22 = PDstandardNth1(gt22, i, j, k);
+    CCTK_REAL const PDstandardNth2gt22 = PDstandardNth2(gt22, i, j, k);
+    CCTK_REAL const PDstandardNth3gt22 = PDstandardNth3(gt22, i, j, k);
+    CCTK_REAL const PDstandardNth1gt23 = PDstandardNth1(gt23, i, j, k);
+    CCTK_REAL const PDstandardNth2gt23 = PDstandardNth2(gt23, i, j, k);
+    CCTK_REAL const PDstandardNth3gt23 = PDstandardNth3(gt23, i, j, k);
+    CCTK_REAL const PDstandardNth1gt33 = PDstandardNth1(gt33, i, j, k);
+    CCTK_REAL const PDstandardNth2gt33 = PDstandardNth2(gt33, i, j, k);
+    CCTK_REAL const PDstandardNth3gt33 = PDstandardNth3(gt33, i, j, k);
+    CCTK_REAL const PDstandardNth1phi = PDstandardNth1(phi, i, j, k);
+    CCTK_REAL const PDstandardNth2phi = PDstandardNth2(phi, i, j, k);
+    CCTK_REAL const PDstandardNth3phi = PDstandardNth3(phi, i, j, k);
+    CCTK_REAL const PDstandardNth1trK = PDstandardNth1(trK, i, j, k);
+    CCTK_REAL const PDstandardNth2trK = PDstandardNth2(trK, i, j, k);
+    CCTK_REAL const PDstandardNth3trK = PDstandardNth3(trK, i, j, k);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL detgt = 1;
+    
+    CCTK_REAL gtu11 = INV(detgt)*(gt22L*gt33L - SQR(gt23L));
+    
+    CCTK_REAL gtu12 = (gt13L*gt23L - gt12L*gt33L)*INV(detgt);
+    
+    CCTK_REAL gtu13 = (-(gt13L*gt22L) + gt12L*gt23L)*INV(detgt);
+    
+    CCTK_REAL gtu22 = INV(detgt)*(gt11L*gt33L - SQR(gt13L));
+    
+    CCTK_REAL gtu23 = (gt12L*gt13L - gt11L*gt23L)*INV(detgt);
+    
+    CCTK_REAL gtu33 = INV(detgt)*(gt11L*gt22L - SQR(gt12L));
+    
+    CCTK_REAL Gt111 = khalf*(gtu11*PDstandardNth1gt11 + 
+      2*(gtu12*PDstandardNth1gt12 + gtu13*PDstandardNth1gt13) - 
+      gtu12*PDstandardNth2gt11 - gtu13*PDstandardNth3gt11);
+    
+    CCTK_REAL Gt211 = khalf*(gtu12*PDstandardNth1gt11 + 
+      2*(gtu22*PDstandardNth1gt12 + gtu23*PDstandardNth1gt13) - 
+      gtu22*PDstandardNth2gt11 - gtu23*PDstandardNth3gt11);
+    
+    CCTK_REAL Gt311 = khalf*(gtu13*PDstandardNth1gt11 + 
+      2*(gtu23*PDstandardNth1gt12 + gtu33*PDstandardNth1gt13) - 
+      gtu23*PDstandardNth2gt11 - gtu33*PDstandardNth3gt11);
+    
+    CCTK_REAL Gt112 = khalf*(gtu12*PDstandardNth1gt22 + 
+      gtu11*PDstandardNth2gt11 + gtu13*(PDstandardNth1gt23 + 
+      PDstandardNth2gt13 - PDstandardNth3gt12));
+    
+    CCTK_REAL Gt212 = khalf*(gtu22*PDstandardNth1gt22 + 
+      gtu12*PDstandardNth2gt11 + gtu23*(PDstandardNth1gt23 + 
+      PDstandardNth2gt13 - PDstandardNth3gt12));
+    
+    CCTK_REAL Gt312 = khalf*(gtu23*PDstandardNth1gt22 + 
+      gtu13*PDstandardNth2gt11 + gtu33*(PDstandardNth1gt23 + 
+      PDstandardNth2gt13 - PDstandardNth3gt12));
+    
+    CCTK_REAL Gt113 = khalf*(gtu13*PDstandardNth1gt33 + 
+      gtu11*PDstandardNth3gt11 + gtu12*(PDstandardNth1gt23 - 
+      PDstandardNth2gt13 + PDstandardNth3gt12));
+    
+    CCTK_REAL Gt213 = khalf*(gtu23*PDstandardNth1gt33 + 
+      gtu12*PDstandardNth3gt11 + gtu22*(PDstandardNth1gt23 - 
+      PDstandardNth2gt13 + PDstandardNth3gt12));
+    
+    CCTK_REAL Gt313 = khalf*(gtu33*PDstandardNth1gt33 + 
+      gtu13*PDstandardNth3gt11 + gtu23*(PDstandardNth1gt23 - 
+      PDstandardNth2gt13 + PDstandardNth3gt12));
+    
+    CCTK_REAL Gt122 = khalf*(gtu11*(-PDstandardNth1gt22 + 
+      2*PDstandardNth2gt12) + gtu12*PDstandardNth2gt22 + 
+      gtu13*(2*PDstandardNth2gt23 - PDstandardNth3gt22));
+    
+    CCTK_REAL Gt222 = khalf*(gtu12*(-PDstandardNth1gt22 + 
+      2*PDstandardNth2gt12) + gtu22*PDstandardNth2gt22 + 
+      gtu23*(2*PDstandardNth2gt23 - PDstandardNth3gt22));
+    
+    CCTK_REAL Gt322 = khalf*(gtu13*(-PDstandardNth1gt22 + 
+      2*PDstandardNth2gt12) + gtu23*PDstandardNth2gt22 + 
+      gtu33*(2*PDstandardNth2gt23 - PDstandardNth3gt22));
+    
+    CCTK_REAL Gt123 = khalf*(gtu13*PDstandardNth2gt33 + 
+      gtu11*(-PDstandardNth1gt23 + PDstandardNth2gt13 + PDstandardNth3gt12) + 
+      gtu12*PDstandardNth3gt22);
+    
+    CCTK_REAL Gt223 = khalf*(gtu23*PDstandardNth2gt33 + 
+      gtu12*(-PDstandardNth1gt23 + PDstandardNth2gt13 + PDstandardNth3gt12) + 
+      gtu22*PDstandardNth3gt22);
+    
+    CCTK_REAL Gt323 = khalf*(gtu33*PDstandardNth2gt33 + 
+      gtu13*(-PDstandardNth1gt23 + PDstandardNth2gt13 + PDstandardNth3gt12) + 
+      gtu23*PDstandardNth3gt22);
+    
+    CCTK_REAL Gt133 = khalf*(-(gtu11*PDstandardNth1gt33) - 
+      gtu12*PDstandardNth2gt33 + 2*gtu11*PDstandardNth3gt13 + 
+      2*gtu12*PDstandardNth3gt23 + gtu13*PDstandardNth3gt33);
+    
+    CCTK_REAL Gt233 = khalf*(-(gtu12*PDstandardNth1gt33) - 
+      gtu22*PDstandardNth2gt33 + 2*gtu12*PDstandardNth3gt13 + 
+      2*gtu22*PDstandardNth3gt23 + gtu23*PDstandardNth3gt33);
+    
+    CCTK_REAL Gt333 = khalf*(-(gtu13*PDstandardNth1gt33) - 
+      gtu23*PDstandardNth2gt33 + 2*gtu13*PDstandardNth3gt13 + 
+      2*gtu23*PDstandardNth3gt23 + gtu33*PDstandardNth3gt33);
+    
+    CCTK_REAL fac1 = IfThen(conformalMethod,-(khalf*INV(phiL)),1);
+    
+    CCTK_REAL cdphi1 = fac1*PDstandardNth1phi;
+    
+    CCTK_REAL cdphi2 = fac1*PDstandardNth2phi;
+    
+    CCTK_REAL cdphi3 = fac1*PDstandardNth3phi;
+    
+    CCTK_REAL S1 = (-eTtxL + beta1L*eTxxL + beta2L*eTxyL + 
+      beta3L*eTxzL)*INV(alphaL);
+    
+    CCTK_REAL S2 = (-eTtyL + beta1L*eTxyL + beta2L*eTyyL + 
+      beta3L*eTyzL)*INV(alphaL);
+    
+    CCTK_REAL S3 = (-eTtzL + beta1L*eTxzL + beta2L*eTyzL + 
+      beta3L*eTzzL)*INV(alphaL);
+    
+    CCTK_REAL M1L = -2.*((At12L*Gt211 + At13L*Gt311)*gtu11 + 
+      At11L*Gt123*gtu23) + At12L*((6.*cdphi1 - 1.*Gt111)*gtu12 - 
+      3.*Gt213*gtu13 + 6.*(cdphi2*gtu22 + cdphi3*gtu23) - 1.*Gt233*gtu33) - 
+      1.*((At22L*Gt212 + At12L*(Gt112 + Gt222) + At23L*Gt312)*gtu22 + 
+      (At13L*Gt112 + At12L*Gt113 + At23L*Gt212)*gtu23 + (At13L*Gt113 + 
+      At23L*Gt213 + At33L*Gt313)*gtu33) + At11L*((6.*cdphi1 - 2.*Gt111)*gtu11 
+      + 6.*(cdphi2*gtu12 + cdphi3*gtu13) - 1.*(Gt122*gtu22 + Gt133*gtu33)) + 
+      At13L*((6.*cdphi1 - 1.*Gt111)*gtu13 + 6.*(cdphi2*gtu23 + cdphi3*gtu33) 
+      - 1.*(Gt322*gtu22 + Gt333*gtu33)) + gtu11*PDstandardNth1At11 - 
+      0.6666666666666666666666666666666666666667*PDstandardNth1trK + 
+      gtu12*(-1.*(At22L*Gt211 + At23L*Gt311) - 3.*(At11L*Gt112 + At12L*Gt212 
+      + At13L*Gt312) + PDstandardNth1At12 + PDstandardNth2At11) + 
+      gtu22*PDstandardNth2At12 + gtu13*(-1.*(At23L*Gt211 + At33L*Gt311) - 
+      3.*(At11L*Gt113 + At13L*Gt313) + PDstandardNth1At13 + 
+      PDstandardNth3At11) + gtu23*(-1.*(At22L*Gt213 + At33L*Gt312 + 
+      At23L*Gt313) - 2.*(At12L*Gt223 + At13L*Gt323) + PDstandardNth2At13 + 
+      PDstandardNth3At12) + gtu33*PDstandardNth3At13 - 
+      25.13274122871834590770114706623602307358*S1;
+    
+    CCTK_REAL M2L = At12L*((6.*cdphi1 - 1.*Gt111)*gtu11 + 6.*(cdphi2*gtu12 
+      + cdphi3*gtu13) - 2.*Gt122*gtu22 - 3.*Gt123*gtu23 - 1.*Gt133*gtu33) + 
+      At22L*((6.*cdphi2 - 2.*Gt222)*gtu22 + 6.*cdphi3*gtu23 - 1.*Gt233*gtu33) 
+      + At23L*(-2.*Gt322*gtu22 - 1.*Gt333*gtu33 + 6.*(cdphi2*gtu23 + 
+      cdphi3*gtu33)) - 1.*((At11L*Gt112 + At22L*Gt211 + At12L*Gt212 + 
+      At23L*Gt311 + At13L*Gt312)*gtu11 + Gt122*(At11L*gtu12 + At13L*gtu23) + 
+      (At23L*Gt223 + At33L*Gt323)*gtu33 + At13L*(Gt112*gtu13 + Gt123*gtu33)) 
+      + gtu11*PDstandardNth1At12 + gtu12*(At22L*(6.*cdphi1 - 3.*Gt212) + 
+      At12L*(-3.*Gt112 - 1.*Gt222) - 3.*At23L*Gt312 - 1.*At13L*Gt322 + 
+      PDstandardNth1At22 + PDstandardNth2At12) + gtu22*PDstandardNth2At22 - 
+      0.6666666666666666666666666666666666666667*PDstandardNth2trK + 
+      gtu13*(-2.*(At12L*Gt113 + At22L*Gt213) + At23L*(6.*cdphi1 - 1.*Gt212 - 
+      2.*Gt313) - 1.*(At11L*Gt123 + At12L*Gt223 + At33L*Gt312 + At13L*Gt323) 
+      + PDstandardNth1At23 + PDstandardNth3At12) + gtu23*(-1.*(At23L*Gt222 + 
+      At33L*Gt322) - 3.*(At22L*Gt223 + At23L*Gt323) + PDstandardNth2At23 + 
+      PDstandardNth3At22) + gtu33*PDstandardNth3At23 - 
+      25.13274122871834590770114706623602307358*S2;
+    
+    CCTK_REAL M3L = -1.*((At11L*Gt113 + At23L*Gt211 + At12L*Gt213 + 
+      At33L*Gt311)*gtu11 + (At22L*Gt223 + At33L*Gt322 + At23L*Gt323)*gtu22 + 
+      At12L*(Gt113*gtu12 + Gt123*gtu22) + Gt133*(At11L*gtu13 + At12L*gtu23)) 
+      + At13L*((6.*cdphi1 - 1.*(Gt111 + Gt313))*gtu11 + 6.*(cdphi2*gtu12 + 
+      cdphi3*gtu13) - 1.*Gt122*gtu22 - 3.*Gt123*gtu23 - 2.*Gt133*gtu33) + 
+      At23L*((6.*cdphi2 - 1.*Gt222)*gtu22 + 6.*cdphi3*gtu23 - 2.*Gt233*gtu33) 
+      + gtu11*PDstandardNth1At13 + gtu12*(-2.*(At13L*Gt112 + At33L*Gt312) + 
+      At23L*(6.*cdphi1 - 2.*Gt212 - 1.*Gt313) - 1.*(At11L*Gt123 + At22L*Gt213 
+      + At12L*Gt223 + At13L*Gt323) + PDstandardNth1At23 + PDstandardNth2At13) 
+      + gtu22*PDstandardNth2At23 + gtu13*(-3.*(At13L*Gt113 + At23L*Gt213) + 
+      At33L*(6.*cdphi1 - 3.*Gt313) - 1.*(At12L*Gt233 + At13L*Gt333) + 
+      PDstandardNth1At33 + PDstandardNth3At13) + gtu23*(-1.*At22L*Gt233 + 
+      At33L*(6.*cdphi2 - 3.*Gt323) + At23L*(-3.*Gt223 - 1.*Gt333) + 
+      PDstandardNth2At33 + PDstandardNth3At23) + gtu33*(At33L*(6.*cdphi3 - 
+      2.*Gt333) + PDstandardNth3At33) - 
+      0.6666666666666666666666666666666666666667*PDstandardNth3trK - 
+      25.13274122871834590770114706623602307358*S3;
+    
+    CCTK_REAL cAL = At11L*gtu11 + At22L*gtu22 + 2*(At12L*gtu12 + 
+      At13L*gtu13 + At23L*gtu23) + At33L*gtu33;
+    
+    
+    /* Copy local copies back to grid functions */
+    cA[index] = cAL;
+    M1[index] = M1L;
+    M2[index] = M2L;
+    M3[index] = M3L;
+  }
+  LC_ENDLOOP3 (ML_BSSN_constraints2);
+}
+
+void ML_BSSN_constraints2(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  GenericFD_LoopOverInterior(cctkGH, &ML_BSSN_constraints2_Body);
+}
diff --git a/ML_BSSN/src/ML_BSSN_constraints_boundary.c b/ML_BSSN/src/ML_BSSN_constraints_boundary.c
deleted file mode 100644
index 798a19c..0000000
--- a/ML_BSSN/src/ML_BSSN_constraints_boundary.c
+++ /dev/null
@@ -1,162 +0,0 @@
-/*  File produced by Kranc */
-
-#define KRANC_C
-
-#include <assert.h>
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-#include "GenericFD.h"
-#include "Differencing.h"
-#include "loopcontrol.h"
-
-/* Define macros used in calculations */
-#define INITVALUE  (42)
-#define INV(x) ((1.0) / (x))
-#define SQR(x) ((x) * (x))
-#define CUB(x) ((x) * (x) * (x))
-#define QAD(x) ((x) * (x) * (x) * (x))
-
-void ML_BSSN_constraints_boundary_SelectBCs(CCTK_ARGUMENTS)
-{
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
-  
-  CCTK_INT ierr = 0;
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_cons_detg","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_cons_detg.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_cons_Gamma","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_cons_Gamma.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_cons_traceA","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_cons_traceA.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_Ham","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_Ham.");
-  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_mom","flat");
-  if (ierr < 0)
-    CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_mom.");
-  return;
-}
-
-void ML_BSSN_constraints_boundary_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
-{
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
-  
-  
-  /* Declare finite differencing variables */
-  
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_constraints_boundary_Body");
-  }
-  
-  if (cctk_iteration % ML_BSSN_constraints_boundary_calc_every != ML_BSSN_constraints_boundary_calc_offset)
-  {
-    return;
-  }
-  
-  /* Include user-supplied include files */
-  
-  /* Initialise finite differencing variables */
-  CCTK_REAL const dx = CCTK_DELTA_SPACE(0);
-  CCTK_REAL const dy = CCTK_DELTA_SPACE(1);
-  CCTK_REAL const dz = CCTK_DELTA_SPACE(2);
-  int const di = 1;
-  int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
-  int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
-  CCTK_REAL const dxi = 1.0 / dx;
-  CCTK_REAL const dyi = 1.0 / dy;
-  CCTK_REAL const dzi = 1.0 / dz;
-  CCTK_REAL const khalf = 0.5;
-  CCTK_REAL const kthird = 1/3.0;
-  CCTK_REAL const ktwothird = 2.0/3.0;
-  CCTK_REAL const kfourthird = 4.0/3.0;
-  CCTK_REAL const keightthird = 8.0/3.0;
-  CCTK_REAL const hdxi = 0.5 * dxi;
-  CCTK_REAL const hdyi = 0.5 * dyi;
-  CCTK_REAL const hdzi = 0.5 * dzi;
-  
-  /* Initialize predefined quantities */
-  CCTK_REAL const p1o12dx = INV(dx)/12.;
-  CCTK_REAL const p1o12dy = INV(dy)/12.;
-  CCTK_REAL const p1o12dz = INV(dz)/12.;
-  CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.;
-  CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.;
-  CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.;
-  CCTK_REAL const p1o24dx = INV(dx)/24.;
-  CCTK_REAL const p1o24dy = INV(dy)/24.;
-  CCTK_REAL const p1o24dz = INV(dz)/24.;
-  CCTK_REAL const p1o64dx = INV(dx)/64.;
-  CCTK_REAL const p1o64dy = INV(dy)/64.;
-  CCTK_REAL const p1o64dz = INV(dz)/64.;
-  CCTK_REAL const p1odx = INV(dx);
-  CCTK_REAL const p1ody = INV(dy);
-  CCTK_REAL const p1odz = INV(dz);
-  CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.;
-  CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.;
-  CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.;
-  
-  /* Loop over the grid points */
-  #pragma omp parallel
-  LC_LOOP3 (ML_BSSN_constraints_boundary,
-            i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],
-            cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
-  {
-    // int index = INITVALUE;
-    int const index = CCTK_GFINDEX3D(cctkGH,i,j,k);
-    /* Declare derivatives */
-    
-    /* Assign local copies of grid functions */
-    
-    /* Include user supplied include files */
-    
-    /* Precompute derivatives */
-    
-    /* Calculate temporaries and grid functions */
-    CCTK_REAL HL = 0;
-    
-    CCTK_REAL M1L = 0;
-    
-    CCTK_REAL M2L = 0;
-    
-    CCTK_REAL M3L = 0;
-    
-    CCTK_REAL cSL = 0;
-    
-    CCTK_REAL cXt1L = 0;
-    
-    CCTK_REAL cXt2L = 0;
-    
-    CCTK_REAL cXt3L = 0;
-    
-    CCTK_REAL cAL = 0;
-    
-    
-    /* Copy local copies back to grid functions */
-    cA[index] = cAL;
-    cS[index] = cSL;
-    cXt1[index] = cXt1L;
-    cXt2[index] = cXt2L;
-    cXt3[index] = cXt3L;
-    H[index] = HL;
-    M1[index] = M1L;
-    M2[index] = M2L;
-    M3[index] = M3L;
-  }
-  LC_ENDLOOP3 (ML_BSSN_constraints_boundary);
-}
-
-void ML_BSSN_constraints_boundary(CCTK_ARGUMENTS)
-{
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
-  
-  GenericFD_LoopOverBoundaryWithGhosts(cctkGH, &ML_BSSN_constraints_boundary_Body);
-}
diff --git a/ML_BSSN/src/ML_BSSN_convertFromADMBase.c b/ML_BSSN/src/ML_BSSN_convertFromADMBase.c
index b0f0b93..91129d0 100644
--- a/ML_BSSN/src/ML_BSSN_convertFromADMBase.c
+++ b/ML_BSSN/src/ML_BSSN_convertFromADMBase.c
@@ -38,6 +38,9 @@ void ML_BSSN_convertFromADMBase_Body(cGH const * restrict const cctkGH, int cons
     return;
   }
   
+  const char *groups[] = {"ADMBase::curv","ADMBase::lapse","ADMBase::metric","ADMBase::shift","ML_BSSN::ML_curv","ML_BSSN::ML_lapse","ML_BSSN::ML_log_confac","ML_BSSN::ML_metric","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_convertFromADMBase", 10, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
diff --git a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
index 6d6832c..a8d3f99 100644
--- a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
+++ b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
@@ -56,6 +56,9 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH const * restrict const cctkGH, int
     return;
   }
   
+  const char *groups[] = {"ADMBase::dtlapse","ADMBase::dtshift","grid::coordinates","Grid::coordinates","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtshift","ML_BSSN::ML_Gamma","ML_BSSN::ML_lapse","ML_BSSN::ML_metric","ML_BSSN::ML_shift"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_convertFromADMBaseGamma", 10, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
@@ -344,11 +347,11 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH const * restrict const cctkGH, int
     CCTK_REAL B1L = 12*IfThen(ShiftBCoeff*ShiftGammaCoeff != 
       0,dtbetaxL*INV(ShiftGammaCoeff)*INV(theta),0) + 
       IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetaxL - 
-      beta1L*PDupwindNthAnti1beta1*ShiftAdvectionCoeff - 
-      PDupwindNthSymm1beta1*ShiftAdvectionCoeff*Abs(beta1L))*INV(ShiftGammaCoeff)*INV(theta),0) 
+      ShiftAdvectionCoeff*(beta1L*PDupwindNthAnti1beta1 + 
+      PDupwindNthSymm1beta1*Abs(beta1L)))*INV(ShiftGammaCoeff)*INV(theta),0) 
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetaxL - 
-      beta2L*PDupwindNthAnti2beta1*ShiftAdvectionCoeff - 
-      PDupwindNthSymm2beta1*ShiftAdvectionCoeff*Abs(beta2L))*INV(ShiftGammaCoeff)*INV(theta),0) 
+      ShiftAdvectionCoeff*(beta2L*PDupwindNthAnti2beta1 + 
+      PDupwindNthSymm2beta1*Abs(beta2L)))*INV(ShiftGammaCoeff)*INV(theta),0) 
       + 2*(IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetaxL - 
       beta1L*PDupwindNthAnti1beta1*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff)*INV(theta),0) 
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetaxL - 
@@ -362,17 +365,17 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH const * restrict const cctkGH, int
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetaxL - 
       PDupwindNthSymm3beta1*ShiftAdvectionCoeff*Abs(beta3L))*INV(ShiftGammaCoeff)*INV(theta),0)) 
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetaxL - 
-      beta3L*PDupwindNthAnti3beta1*ShiftAdvectionCoeff - 
-      PDupwindNthSymm3beta1*ShiftAdvectionCoeff*Abs(beta3L))*INV(ShiftGammaCoeff)*INV(theta),0);
+      ShiftAdvectionCoeff*(beta3L*PDupwindNthAnti3beta1 + 
+      PDupwindNthSymm3beta1*Abs(beta3L)))*INV(ShiftGammaCoeff)*INV(theta),0);
     
     CCTK_REAL B2L = 12*IfThen(ShiftBCoeff*ShiftGammaCoeff != 
       0,dtbetayL*INV(ShiftGammaCoeff)*INV(theta),0) + 
       IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetayL - 
-      beta1L*PDupwindNthAnti1beta2*ShiftAdvectionCoeff - 
-      PDupwindNthSymm1beta2*ShiftAdvectionCoeff*Abs(beta1L))*INV(ShiftGammaCoeff)*INV(theta),0) 
+      ShiftAdvectionCoeff*(beta1L*PDupwindNthAnti1beta2 + 
+      PDupwindNthSymm1beta2*Abs(beta1L)))*INV(ShiftGammaCoeff)*INV(theta),0) 
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetayL - 
-      beta2L*PDupwindNthAnti2beta2*ShiftAdvectionCoeff - 
-      PDupwindNthSymm2beta2*ShiftAdvectionCoeff*Abs(beta2L))*INV(ShiftGammaCoeff)*INV(theta),0) 
+      ShiftAdvectionCoeff*(beta2L*PDupwindNthAnti2beta2 + 
+      PDupwindNthSymm2beta2*Abs(beta2L)))*INV(ShiftGammaCoeff)*INV(theta),0) 
       + 2*(IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetayL - 
       beta1L*PDupwindNthAnti1beta2*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff)*INV(theta),0) 
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetayL - 
@@ -386,17 +389,17 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH const * restrict const cctkGH, int
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetayL - 
       PDupwindNthSymm3beta2*ShiftAdvectionCoeff*Abs(beta3L))*INV(ShiftGammaCoeff)*INV(theta),0)) 
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetayL - 
-      beta3L*PDupwindNthAnti3beta2*ShiftAdvectionCoeff - 
-      PDupwindNthSymm3beta2*ShiftAdvectionCoeff*Abs(beta3L))*INV(ShiftGammaCoeff)*INV(theta),0);
+      ShiftAdvectionCoeff*(beta3L*PDupwindNthAnti3beta2 + 
+      PDupwindNthSymm3beta2*Abs(beta3L)))*INV(ShiftGammaCoeff)*INV(theta),0);
     
     CCTK_REAL B3L = 12*IfThen(ShiftBCoeff*ShiftGammaCoeff != 
       0,dtbetazL*INV(ShiftGammaCoeff)*INV(theta),0) + 
       IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetazL - 
-      beta1L*PDupwindNthAnti1beta3*ShiftAdvectionCoeff - 
-      PDupwindNthSymm1beta3*ShiftAdvectionCoeff*Abs(beta1L))*INV(ShiftGammaCoeff)*INV(theta),0) 
+      ShiftAdvectionCoeff*(beta1L*PDupwindNthAnti1beta3 + 
+      PDupwindNthSymm1beta3*Abs(beta1L)))*INV(ShiftGammaCoeff)*INV(theta),0) 
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetazL - 
-      beta2L*PDupwindNthAnti2beta3*ShiftAdvectionCoeff - 
-      PDupwindNthSymm2beta3*ShiftAdvectionCoeff*Abs(beta2L))*INV(ShiftGammaCoeff)*INV(theta),0) 
+      ShiftAdvectionCoeff*(beta2L*PDupwindNthAnti2beta3 + 
+      PDupwindNthSymm2beta3*Abs(beta2L)))*INV(ShiftGammaCoeff)*INV(theta),0) 
       + 2*(IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetazL - 
       beta1L*PDupwindNthAnti1beta3*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff)*INV(theta),0) 
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetazL - 
@@ -410,8 +413,8 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH const * restrict const cctkGH, int
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetazL - 
       PDupwindNthSymm3beta3*ShiftAdvectionCoeff*Abs(beta3L))*INV(ShiftGammaCoeff)*INV(theta),0)) 
       + IfThen(ShiftBCoeff*ShiftGammaCoeff != 0,(dtbetazL - 
-      beta3L*PDupwindNthAnti3beta3*ShiftAdvectionCoeff - 
-      PDupwindNthSymm3beta3*ShiftAdvectionCoeff*Abs(beta3L))*INV(ShiftGammaCoeff)*INV(theta),0);
+      ShiftAdvectionCoeff*(beta3L*PDupwindNthAnti3beta3 + 
+      PDupwindNthSymm3beta3*Abs(beta3L)))*INV(ShiftGammaCoeff)*INV(theta),0);
     
     
     /* Copy local copies back to grid functions */
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBase.c b/ML_BSSN/src/ML_BSSN_convertToADMBase.c
index 95058c1..9aac628 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBase.c
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBase.c
@@ -38,6 +38,9 @@ void ML_BSSN_convertToADMBase_Body(cGH const * restrict const cctkGH, int const
     return;
   }
   
+  const char *groups[] = {"ADMBase::curv","ADMBase::lapse","ADMBase::metric","ADMBase::shift","ML_BSSN::ML_curv","ML_BSSN::ML_lapse","ML_BSSN::ML_log_confac","ML_BSSN::ML_metric","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_convertToADMBase", 10, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.c b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.c
index b8fe647..4122892 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.c
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.c
@@ -53,6 +53,9 @@ void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const cctkGH
     return;
   }
   
+  const char *groups[] = {"ADMBase::dtlapse","ADMBase::dtshift","grid::coordinates","Grid::coordinates","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtshift","ML_BSSN::ML_Gamma","ML_BSSN::ML_lapse","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_convertToADMBaseDtLapseShift", 10, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c
index 72c6ff2..98d402c 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c
@@ -53,6 +53,9 @@ void ML_BSSN_convertToADMBaseDtLapseShiftBoundary_Body(cGH const * restrict cons
     return;
   }
   
+  const char *groups[] = {"ADMBase::dtlapse","ADMBase::dtshift","grid::coordinates","Grid::coordinates","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtshift","ML_BSSN::ML_Gamma","ML_BSSN::ML_lapse","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_convertToADMBaseDtLapseShiftBoundary", 10, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.c b/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.c
index e2d008c..d4a0c76 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.c
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.c
@@ -38,6 +38,9 @@ void ML_BSSN_convertToADMBaseFakeDtLapseShift_Body(cGH const * restrict const cc
     return;
   }
   
+  const char *groups[] = {"ADMBase::dtlapse","ADMBase::dtshift","grid::coordinates","Grid::coordinates","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtshift","ML_BSSN::ML_Gamma","ML_BSSN::ML_lapse","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_convertToADMBaseFakeDtLapseShift", 10, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
diff --git a/ML_BSSN/src/ML_BSSN_enforce.c b/ML_BSSN/src/ML_BSSN_enforce.c
index d8b0d21..0e7b581 100644
--- a/ML_BSSN/src/ML_BSSN_enforce.c
+++ b/ML_BSSN/src/ML_BSSN_enforce.c
@@ -38,6 +38,9 @@ void ML_BSSN_enforce_Body(cGH const * restrict const cctkGH, int const dir, int
     return;
   }
   
+  const char *groups[] = {"ML_BSSN::ML_curv","ML_BSSN::ML_lapse","ML_BSSN::ML_metric"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_enforce", 3, groups);
+  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
diff --git a/ML_BSSN/src/make.code.defn b/ML_BSSN/src/make.code.defn
index 3f0d62e..7895da0 100644
--- a/ML_BSSN/src/make.code.defn
+++ b/ML_BSSN/src/make.code.defn
@@ -1,3 +1,3 @@
 # File produced by Kranc
 
-SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_BSSN_Minkowski.c ML_BSSN_convertFromADMBase.c ML_BSSN_convertFromADMBaseGamma.c ML_BSSN_RHS1.c ML_BSSN_RHS2.c ML_BSSN_RHSStaticBoundary.c ML_BSSN_RHSRadiativeBoundary.c ML_BSSN_enforce.c ML_BSSN_boundary.c ML_BSSN_convertToADMBase.c ML_BSSN_convertToADMBaseDtLapseShift.c ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c ML_BSSN_convertToADMBaseFakeDtLapseShift.c ML_BSSN_constraints.c ML_BSSN_constraints_boundary.c Boundaries.c 
+SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_BSSN_Minkowski.c ML_BSSN_convertFromADMBase.c ML_BSSN_convertFromADMBaseGamma.c ML_BSSN_RHS1.c ML_BSSN_RHS2.c ML_BSSN_RHSStaticBoundary.c ML_BSSN_enforce.c ML_BSSN_boundary.c ML_BSSN_convertToADMBase.c ML_BSSN_convertToADMBaseDtLapseShift.c ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c ML_BSSN_convertToADMBaseFakeDtLapseShift.c ML_BSSN_constraints1.c ML_BSSN_constraints2.c Boundaries.c