Update McLachlan. ML_BSSN_Helper is now required at run time

Calculate the complete RHS in one loop, including matter if desired.
This requires different McLachlan builds for runs with and without
matter.

Make the number of active timelevels a run-time parameter.

Register the slicing in ML_BSSN_Helper.  This requires CoordGauge to be
activated at run time.

Add custom boundary conditions, to be used instead of Cactus boundary
conditions.

Make McLachlan fully compliant with the ADMBase initial data mechanism.

Calculate RHS and constrains in special groups provided by ML_BSSN_Helper.
ML_BSSN_Helper is now required at run time.

5 files changed, 424 insertions, 17 deletions

diff --git a/ML_ADM/param.ccl b/ML_ADM/param.ccl
index 6ca561f..df67e04 100644
--- a/ML_ADM/param.ccl
+++ b/ML_ADM/param.ccl
@@ -32,11 +32,11 @@ KEYWORD my_initial_data "my_initial_data"
 } "ADMBase"
 
 private:
-KEYWORD SpaceTime "SpaceTime"
+KEYWORD my_boundary_condition "my_boundary_condition"
 {
-  "Space" :: "Space"
-  "Space+Matter" :: "Space+Matter"
-} "Space"
+  "none" :: "none"
+  "Minkowski" :: "Minkowski"
+} "none"
 
 restricted:
 CCTK_INT ML_ADM_MaxNumEvolvedVars "Number of evolved variables used by this thorn" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Evolved_Vars
@@ -50,6 +50,12 @@ CCTK_INT ML_ADM_MaxNumConstrainedVars "Number of constrained variables used by t
   38:38 :: "Number of constrained variables used by this thorn"
 } 38
 
+private:
+CCTK_INT timelevels "Number of active timelevels"
+{
+  0:3 :: ""
+} 3
+
 restricted:
 CCTK_INT ML_ADM_Minkowski_calc_every "ML_ADM_Minkowski_calc_every"
 {
@@ -69,6 +75,12 @@ CCTK_INT ML_ADM_RHS_calc_every "ML_ADM_RHS_calc_every"
 } 1
 
 restricted:
+CCTK_INT ML_ADM_boundary_calc_every "ML_ADM_boundary_calc_every"
+{
+  *:* :: ""
+} 1
+
+restricted:
 CCTK_INT ML_ADM_convertToADMBase_calc_every "ML_ADM_convertToADMBase_calc_every"
 {
   *:* :: ""
@@ -81,6 +93,12 @@ CCTK_INT ML_ADM_constraints_calc_every "ML_ADM_constraints_calc_every"
 } 1
 
 restricted:
+CCTK_INT ML_ADM_constraints_boundary_calc_every "ML_ADM_constraints_boundary_calc_every"
+{
+  *:* :: ""
+} 1
+
+restricted:
 CCTK_INT ML_ADM_Minkowski_calc_offset "ML_ADM_Minkowski_calc_offset"
 {
   *:* :: ""
@@ -99,6 +117,12 @@ CCTK_INT ML_ADM_RHS_calc_offset "ML_ADM_RHS_calc_offset"
 } 0
 
 restricted:
+CCTK_INT ML_ADM_boundary_calc_offset "ML_ADM_boundary_calc_offset"
+{
+  *:* :: ""
+} 0
+
+restricted:
 CCTK_INT ML_ADM_convertToADMBase_calc_offset "ML_ADM_convertToADMBase_calc_offset"
 {
   *:* :: ""
@@ -110,6 +134,12 @@ CCTK_INT ML_ADM_constraints_calc_offset "ML_ADM_constraints_calc_offset"
   *:* :: ""
 } 0
 
+restricted:
+CCTK_INT ML_ADM_constraints_boundary_calc_offset "ML_ADM_constraints_boundary_calc_offset"
+{
+  *:* :: ""
+} 0
+
 private:
 KEYWORD K11_bound "Boundary condition to implement"
 {
diff --git a/ML_ADM/schedule.ccl b/ML_ADM/schedule.ccl
index 36a4cdd..5dbca42 100644
--- a/ML_ADM/schedule.ccl
+++ b/ML_ADM/schedule.ccl
@@ -16,13 +16,57 @@ STORAGE: ml_metricrhs[1]
 
 STORAGE: ml_shiftrhs[1]
 
-STORAGE: ml_curv[3]
+if (timelevels == 1)
+{
+  STORAGE: ml_curv[1]
+}
+if (timelevels == 2)
+{
+  STORAGE: ml_curv[2]
+}
+if (timelevels == 3)
+{
+  STORAGE: ml_curv[3]
+}
 
-STORAGE: ml_lapse[3]
+if (timelevels == 1)
+{
+  STORAGE: ml_lapse[1]
+}
+if (timelevels == 2)
+{
+  STORAGE: ml_lapse[2]
+}
+if (timelevels == 3)
+{
+  STORAGE: ml_lapse[3]
+}
 
-STORAGE: ml_metric[3]
+if (timelevels == 1)
+{
+  STORAGE: ml_metric[1]
+}
+if (timelevels == 2)
+{
+  STORAGE: ml_metric[2]
+}
+if (timelevels == 3)
+{
+  STORAGE: ml_metric[3]
+}
 
-STORAGE: ml_shift[3]
+if (timelevels == 1)
+{
+  STORAGE: ml_shift[1]
+}
+if (timelevels == 2)
+{
+  STORAGE: ml_shift[2]
+}
+if (timelevels == 3)
+{
+  STORAGE: ml_shift[3]
+}
 
 schedule ML_ADM_Startup at STARTUP
 {
@@ -48,7 +92,6 @@ if (CCTK_EQUALS(my_initial_data, "Minkowski"))
   schedule ML_ADM_Minkowski IN ADMBase_InitialData
   {
     LANG: C
-    
   } "ML_ADM_Minkowski"
 }
 
@@ -58,40 +101,51 @@ if (CCTK_EQUALS(my_initial_data, "ADMBase"))
   schedule ML_ADM_convertFromADMBase AT initial AFTER ADMBase_PostInitial
   {
     LANG: C
-    
   } "ML_ADM_convertFromADMBase"
 }
 
 schedule ML_ADM_RHS IN MoL_CalcRHS
 {
   LANG: C
-  
 } "ML_ADM_RHS"
 
 schedule ML_ADM_RHS AT analysis
 {
   LANG: C
-  
   SYNC: ml_curvrhs
   SYNC: ml_lapserhs
   SYNC: ml_metricrhs
   SYNC: ml_shiftrhs
 } "ML_ADM_RHS"
 
-schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER ML_ADM_ApplyBCs
+
+if (CCTK_EQUALS(my_boundary_condition, "Minkowski"))
+{
+  schedule ML_ADM_boundary IN MoL_PostStep
+  {
+    LANG: C
+  } "ML_ADM_boundary"
+}
+
+schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER (ML_ADM_ApplyBCs ML_ADM_boundary)
 {
   LANG: C
-  
 } "ML_ADM_convertToADMBase"
 
 schedule ML_ADM_constraints AT analysis
 {
   LANG: C
-  
   SYNC: Ham
   SYNC: mom
+  TRIGGERS: Ham
+  TRIGGERS: mom
 } "ML_ADM_constraints"
 
+schedule ML_ADM_constraints_boundary AT analysis AFTER ML_ADM_constraints
+{
+  LANG: C
+} "ML_ADM_constraints_boundary"
+
 schedule ML_ADM_ApplyBoundConds in MoL_PostStep
 {
   LANG: C
@@ -111,5 +165,4 @@ schedule ML_ADM_CheckBoundaries at BASEGRID
 schedule group ApplyBCs as ML_ADM_ApplyBCs in MoL_PostStep  after ML_ADM_ApplyBoundConds
 {
   # no language specified
-  
 } "Apply boundary conditions controlled by thorn Boundary"
diff --git a/ML_ADM/src/ML_ADM_boundary.c b/ML_ADM/src/ML_ADM_boundary.c
new file mode 100644
index 0000000..02c77aa
--- /dev/null
+++ b/ML_ADM/src/ML_ADM_boundary.c
@@ -0,0 +1,181 @@
+/*  File produced by user eschnett */
+/*  Produced with Mathematica Version 6.0 for Mac OS X x86 (32-bit) (April 20, 2007) */
+
+/*  Mathematica script written by Ian Hinder and Sascha Husa */
+
+#define KRANC_C
+
+#include <math.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_ADM_boundary_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[3], CCTK_REAL tangentA[3], CCTK_REAL tangentB[3], CCTK_INT min[3], CCTK_INT max[3], CCTK_INT n_subblock_gfs, CCTK_REAL *subblock_gfs[])
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  
+  /* Declare finite differencing variables */
+  CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE;
+  CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE;
+  CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE;
+  CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE;
+  
+  
+  /* Declare predefined quantities */
+  CCTK_REAL p1o12dx = INITVALUE;
+  CCTK_REAL p1o12dy = INITVALUE;
+  CCTK_REAL p1o12dz = INITVALUE;
+  CCTK_REAL p1o144dxdy = INITVALUE;
+  CCTK_REAL p1o144dxdz = INITVALUE;
+  CCTK_REAL p1o144dydz = INITVALUE;
+  CCTK_REAL pm1o12dx2 = INITVALUE;
+  CCTK_REAL pm1o12dy2 = INITVALUE;
+  CCTK_REAL pm1o12dz2 = INITVALUE;
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_ADM_boundary_Body");
+  }
+  
+  if (cctk_iteration % ML_ADM_boundary_calc_every != ML_ADM_boundary_calc_offset)
+  {
+    return;
+  }
+  
+  /* Include user-supplied include files */
+  
+  /* Initialise finite differencing variables */
+  dx = CCTK_DELTA_SPACE(0);
+  dy = CCTK_DELTA_SPACE(1);
+  dz = CCTK_DELTA_SPACE(2);
+  dxi = 1.0 / dx;
+  dyi = 1.0 / dy;
+  dzi = 1.0 / dz;
+  khalf = 0.5;
+  kthird = 1/3.0;
+  ktwothird = 2.0/3.0;
+  kfourthird = 4.0/3.0;
+  keightthird = 8.0/3.0;
+  hdxi = 0.5 * dxi;
+  hdyi = 0.5 * dyi;
+  hdzi = 0.5 * dzi;
+  
+  /* Initialize predefined quantities */
+  p1o12dx = INV(dx)/12.;
+  p1o12dy = INV(dy)/12.;
+  p1o12dz = INV(dz)/12.;
+  p1o144dxdy = (INV(dx)*INV(dy))/144.;
+  p1o144dxdz = (INV(dx)*INV(dz))/144.;
+  p1o144dydz = (INV(dy)*INV(dz))/144.;
+  pm1o12dx2 = -pow(dx,-2)/12.;
+  pm1o12dy2 = -pow(dy,-2)/12.;
+  pm1o12dz2 = -pow(dz,-2)/12.;
+  
+  /* Loop over the grid points */
+  _Pragma ("omp parallel")
+  LC_LOOP3 (ML_ADM_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 subblock_index = INITVALUE;
+    index = CCTK_GFINDEX3D(cctkGH,i,j,k);
+    subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2]));
+    
+    /* Declare shorthands */
+    
+    /* Declare local copies of grid functions */
+    CCTK_REAL alphaL = INITVALUE;
+    CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE;
+    CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE;
+    CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE;
+    /* Declare precomputed derivatives*/
+    
+    /* Declare derivatives */
+    
+    /* Assign local copies of grid functions */
+    
+    /* Assign local copies of subblock grid functions */
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives (new style) */
+    
+    /* Precompute derivatives (old style) */
+    
+    /* Calculate temporaries and grid functions */
+    g11L  =  1;
+    
+    g12L  =  0;
+    
+    g13L  =  0;
+    
+    g22L  =  1;
+    
+    g23L  =  0;
+    
+    g33L  =  1;
+    
+    K11L  =  0;
+    
+    K12L  =  0;
+    
+    K13L  =  0;
+    
+    K22L  =  0;
+    
+    K23L  =  0;
+    
+    K33L  =  0;
+    
+    alphaL  =  1;
+    
+    beta1L  =  0;
+    
+    beta2L  =  0;
+    
+    beta3L  =  0;
+    
+    
+    /* Copy local copies back to grid functions */
+    alpha[index] = alphaL;
+    beta1[index] = beta1L;
+    beta2[index] = beta2L;
+    beta3[index] = beta3L;
+    g11[index] = g11L;
+    g12[index] = g12L;
+    g13[index] = g13L;
+    g22[index] = g22L;
+    g23[index] = g23L;
+    g33[index] = g33L;
+    K11[index] = K11L;
+    K12[index] = K12L;
+    K13[index] = K13L;
+    K22[index] = K22L;
+    K23[index] = K23L;
+    K33[index] = K33L;
+    
+    /* Copy local copies back to subblock grid functions */
+  }
+  LC_ENDLOOP3 (ML_ADM_boundary);
+}
+
+void ML_ADM_boundary(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  GenericFD_LoopOverBoundary(cctkGH, &ML_ADM_boundary_Body);
+}
diff --git a/ML_ADM/src/ML_ADM_constraints_boundary.c b/ML_ADM/src/ML_ADM_constraints_boundary.c
new file mode 100644
index 0000000..7a7053d
--- /dev/null
+++ b/ML_ADM/src/ML_ADM_constraints_boundary.c
@@ -0,0 +1,143 @@
+/*  File produced by user eschnett */
+/*  Produced with Mathematica Version 6.0 for Mac OS X x86 (32-bit) (April 20, 2007) */
+
+/*  Mathematica script written by Ian Hinder and Sascha Husa */
+
+#define KRANC_C
+
+#include <math.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_ADM_constraints_boundary_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[3], CCTK_REAL tangentA[3], CCTK_REAL tangentB[3], CCTK_INT min[3], CCTK_INT max[3], CCTK_INT n_subblock_gfs, CCTK_REAL *subblock_gfs[])
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  
+  /* Declare finite differencing variables */
+  CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE;
+  CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE;
+  CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE;
+  CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE;
+  
+  
+  /* Declare predefined quantities */
+  CCTK_REAL p1o12dx = INITVALUE;
+  CCTK_REAL p1o12dy = INITVALUE;
+  CCTK_REAL p1o12dz = INITVALUE;
+  CCTK_REAL p1o144dxdy = INITVALUE;
+  CCTK_REAL p1o144dxdz = INITVALUE;
+  CCTK_REAL p1o144dydz = INITVALUE;
+  CCTK_REAL pm1o12dx2 = INITVALUE;
+  CCTK_REAL pm1o12dy2 = INITVALUE;
+  CCTK_REAL pm1o12dz2 = INITVALUE;
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_ADM_constraints_boundary_Body");
+  }
+  
+  if (cctk_iteration % ML_ADM_constraints_boundary_calc_every != ML_ADM_constraints_boundary_calc_offset)
+  {
+    return;
+  }
+  
+  /* Include user-supplied include files */
+  
+  /* Initialise finite differencing variables */
+  dx = CCTK_DELTA_SPACE(0);
+  dy = CCTK_DELTA_SPACE(1);
+  dz = CCTK_DELTA_SPACE(2);
+  dxi = 1.0 / dx;
+  dyi = 1.0 / dy;
+  dzi = 1.0 / dz;
+  khalf = 0.5;
+  kthird = 1/3.0;
+  ktwothird = 2.0/3.0;
+  kfourthird = 4.0/3.0;
+  keightthird = 8.0/3.0;
+  hdxi = 0.5 * dxi;
+  hdyi = 0.5 * dyi;
+  hdzi = 0.5 * dzi;
+  
+  /* Initialize predefined quantities */
+  p1o12dx = INV(dx)/12.;
+  p1o12dy = INV(dy)/12.;
+  p1o12dz = INV(dz)/12.;
+  p1o144dxdy = (INV(dx)*INV(dy))/144.;
+  p1o144dxdz = (INV(dx)*INV(dz))/144.;
+  p1o144dydz = (INV(dy)*INV(dz))/144.;
+  pm1o12dx2 = -pow(dx,-2)/12.;
+  pm1o12dy2 = -pow(dy,-2)/12.;
+  pm1o12dz2 = -pow(dz,-2)/12.;
+  
+  /* Loop over the grid points */
+  _Pragma ("omp parallel")
+  LC_LOOP3 (ML_ADM_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 subblock_index = INITVALUE;
+    index = CCTK_GFINDEX3D(cctkGH,i,j,k);
+    subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2]));
+    
+    /* Declare shorthands */
+    
+    /* Declare local copies of grid functions */
+    CCTK_REAL HL = INITVALUE;
+    CCTK_REAL M1L = INITVALUE, M2L = INITVALUE, M3L = INITVALUE;
+    /* Declare precomputed derivatives*/
+    
+    /* Declare derivatives */
+    
+    /* Assign local copies of grid functions */
+    
+    /* Assign local copies of subblock grid functions */
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives (new style) */
+    
+    /* Precompute derivatives (old style) */
+    
+    /* Calculate temporaries and grid functions */
+    HL  =  0;
+    
+    M1L  =  0;
+    
+    M2L  =  0;
+    
+    M3L  =  0;
+    
+    
+    /* Copy local copies back to grid functions */
+    H[index] = HL;
+    M1[index] = M1L;
+    M2[index] = M2L;
+    M3[index] = M3L;
+    
+    /* Copy local copies back to subblock grid functions */
+  }
+  LC_ENDLOOP3 (ML_ADM_constraints_boundary);
+}
+
+void ML_ADM_constraints_boundary(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  GenericFD_LoopOverBoundary(cctkGH, &ML_ADM_constraints_boundary_Body);
+}
diff --git a/ML_ADM/src/make.code.defn b/ML_ADM/src/make.code.defn
index 6e2a3f8..9762fc1 100644
--- a/ML_ADM/src/make.code.defn
+++ b/ML_ADM/src/make.code.defn
@@ -3,4 +3,4 @@
 
 # Mathematica script written by Ian Hinder and Sascha Husa
 
-SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_ADM_Minkowski.c ML_ADM_convertFromADMBase.c ML_ADM_RHS.c ML_ADM_convertToADMBase.c ML_ADM_constraints.c Boundaries.c 
+SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_ADM_Minkowski.c ML_ADM_convertFromADMBase.c ML_ADM_RHS.c ML_ADM_boundary.c ML_ADM_convertToADMBase.c ML_ADM_constraints.c ML_ADM_constraints_boundary.c Boundaries.c