Add CCZ4 Helper thorn.

11 files changed, 488 insertions, 0 deletions

diff --git a/ML_CCZ4_Helper/configuration.ccl b/ML_CCZ4_Helper/configuration.ccl
new file mode 100644
index 0000000..edaa5b7
--- /dev/null
+++ b/ML_CCZ4_Helper/configuration.ccl
@@ -0,0 +1 @@
+REQUIRES THORNS: CartGrid3D CoordGauge
diff --git a/ML_CCZ4_Helper/interface.ccl b/ML_CCZ4_Helper/interface.ccl
new file mode 100644
index 0000000..9b2864f
--- /dev/null
+++ b/ML_CCZ4_Helper/interface.ccl
@@ -0,0 +1,38 @@
+IMPLEMENTS: ML_CCZ4_Helper
+
+INHERITS: ADMBase CoordGauge ML_CCZ4
+
+USES INCLUDE: Slicing.h
+USES INCLUDE: Symmetry.h
+
+
+
+CCTK_INT FUNCTION MoLRegisterConstrainedGroup \
+    (CCTK_INT IN ConstrainedIndex)
+USES FUNCTION MoLRegisterConstrainedGroup
+
+CCTK_INT FUNCTION                         \
+    Boundary_SelectGroupForBC             \
+        (CCTK_POINTER_TO_CONST IN cctkGH, \
+          CCTK_INT IN faces,              \
+          CCTK_INT IN boundary_width,     \
+          CCTK_INT IN table_handle,       \
+          CCTK_STRING IN group_name,      \
+          CCTK_STRING IN bc_name)
+USES FUNCTION Boundary_SelectGroupForBC
+
+CCTK_INT FUNCTION                         \
+    ExtrapolateGammas                     \
+        (CCTK_POINTER_TO_CONST IN cctkGH, \
+         CCTK_REAL ARRAY INOUT var)
+USES FUNCTION ExtrapolateGammas
+
+CCTK_INT FUNCTION                         \
+    NewRad_Apply                          \
+        (CCTK_POINTER_TO_CONST IN cctkGH, \
+         CCTK_REAL ARRAY IN var,          \
+         CCTK_REAL ARRAY INOUT rhs,       \
+         CCTK_REAL IN var0,               \
+         CCTK_REAL IN v0,                 \
+         CCTK_INT IN radpower)
+USES FUNCTION NewRad_Apply
diff --git a/ML_CCZ4_Helper/param.ccl b/ML_CCZ4_Helper/param.ccl
new file mode 100644
index 0000000..8ede1a2
--- /dev/null
+++ b/ML_CCZ4_Helper/param.ccl
@@ -0,0 +1,39 @@
+SHARES: Boundary
+
+USES CCTK_INT radpower
+
+SHARES: MethodofLines
+
+USES CCTK_INT MoL_Num_Constrained_Vars
+
+SHARES: ADMBase
+
+USES KEYWORD evolution_method
+
+SHARES: GenericFD
+
+USES CCTK_INT boundary_width
+
+SHARES: ML_CCZ4
+
+USES CCTK_INT timelevels
+USES CCTK_INT rhs_timelevels
+USES KEYWORD calculate_ADMBase_variables_at
+
+USES CCTK_REAL harmonicF
+
+EXTENDS KEYWORD my_initial_boundary_condition
+{
+  "extrapolate-gammas" :: "extrapolate-gammas"
+}
+EXTENDS KEYWORD my_rhs_boundary_condition
+{
+  "NewRad" :: "NewRad"
+}
+
+PRIVATE:
+
+CCTK_INT ML_CCZ4_MaxNumConstrainedVars "Number of constrained variables used by this thorn" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Constrained_Vars
+{
+  20:20 :: "Number of constrained variables used by this thorn"
+} 20
diff --git a/ML_CCZ4_Helper/schedule.ccl b/ML_CCZ4_Helper/schedule.ccl
new file mode 100644
index 0000000..7e95150
--- /dev/null
+++ b/ML_CCZ4_Helper/schedule.ccl
@@ -0,0 +1,161 @@
+if (CCTK_EQUALS (evolution_method, "ML_CCZ4")) {
+  
+  if (timelevels == 1) {
+    STORAGE: ADMBase::metric[1]
+    STORAGE: ADMBase::curv[1]
+    STORAGE: ADMBase::lapse[1]
+    STORAGE: ADMBase::shift[1]
+    STORAGE: ADMBase::dtlapse[1]
+    STORAGE: ADMBase::dtshift[1]
+  } else if (timelevels == 2) {
+    STORAGE: ADMBase::metric[2]
+    STORAGE: ADMBase::curv[2]
+    STORAGE: ADMBase::lapse[2]
+    STORAGE: ADMBase::shift[2]
+    STORAGE: ADMBase::dtlapse[2]
+    STORAGE: ADMBase::dtshift[2]
+  } else if (timelevels == 3) {
+    STORAGE: ADMBase::metric[3]
+    STORAGE: ADMBase::curv[3]
+    STORAGE: ADMBase::lapse[3]
+    STORAGE: ADMBase::shift[3]
+    STORAGE: ADMBase::dtlapse[3]
+    STORAGE: ADMBase::dtshift[3]
+  } else if (timelevels == 4) {
+    #STORAGE: ADMBase::metric[4]
+    #STORAGE: ADMBase::curv[4]
+    #STORAGE: ADMBase::lapse[4]
+    #STORAGE: ADMBase::shift[4]
+    #STORAGE: ADMBase::dtlapse[4]
+    #STORAGE: ADMBase::dtshift[4]
+    STORAGE: ADMBase::metric[3]
+    STORAGE: ADMBase::curv[3]
+    STORAGE: ADMBase::lapse[3]
+    STORAGE: ADMBase::shift[3]
+    STORAGE: ADMBase::dtlapse[3]
+    STORAGE: ADMBase::dtshift[3]
+  }
+  
+  SCHEDULE ML_CCZ4_RegisterSlicing AT startup
+  {
+    LANG: C
+    OPTIONS: meta
+  } "Register slicing"
+  
+  SCHEDULE ML_CCZ4_SetGroupTags AT startup BEFORE Driver_Startup
+  {
+    LANG: C
+    OPTIONS: meta
+  } "Set checkpointing and prolongation group tags"
+  
+  SCHEDULE ML_CCZ4_RegisterConstrained IN MoL_Register
+  {
+    LANG: C
+    OPTIONS: meta
+  } "Register ADMBase variables as constrained"
+  
+  
+  
+  # Why would any of these two be necessary?
+  #SCHEDULE GROUP ML_CCZ4_evolCalcGroup AT postinitial AFTER MoL_PostStep
+  #{
+  #} "Calculate BSSN RHS"
+  #SCHEDULE GROUP MoL_CalcRHS AT postinitial AFTER MoL_PostStep
+  #{
+  #} "Evaluate RHS"
+  
+  SCHEDULE GROUP ML_CCZ4_evolCalcGroup IN MoL_CalcRHS
+  {
+  } "Calculate BSSN RHS"
+  
+  SCHEDULE GROUP ML_CCZ4_evolCalcGroup AT analysis
+  {
+    TRIGGERS: ML_CCZ4::ML_log_confacrhs
+    TRIGGERS: ML_CCZ4::ML_metricrhs
+    TRIGGERS: ML_CCZ4::ML_Gammarhs
+    TRIGGERS: ML_CCZ4::ML_trace_curvrhs
+    TRIGGERS: ML_CCZ4::ML_curvrhs
+    TRIGGERS: ML_CCZ4::ML_lapserhs
+    TRIGGERS: ML_CCZ4::ML_dtlapserhs
+    TRIGGERS: ML_CCZ4::ML_shiftrhs
+    TRIGGERS: ML_CCZ4::ML_dtshiftrhs
+  } "Calculate BSSN RHS"
+  
+  
+  
+  if (CCTK_EQUALS (my_initial_boundary_condition, "extrapolate-gammas"))
+  {
+    SCHEDULE ML_CCZ4_ExtrapolateGammas AT initial AFTER ML_CCZ4_convertFromADMBaseGamma
+    {
+      LANG: C
+      # We don't need to synchronise here because extrapolation is now filling
+      # ghost zones
+      #SYNC: ML_Gamma
+      #SYNC: ML_dtlapse
+      #SYNC: ML_dtshift
+    } "Extrapolate Gammas and time derivatives of lapse and shift"
+  }
+  
+  if (CCTK_EQUALS (my_rhs_boundary_condition, "NewRad"))
+  {
+    SCHEDULE ML_CCZ4_NewRad IN ML_CCZ4_evolCalcGroup AFTER ML_CCZ4_RHS
+    {
+      LANG: C
+      #SYNC: ML_curvrhs
+      #SYNC: ML_dtlapserhs
+      #SYNC: ML_dtshiftrhs
+      #SYNC: ML_Gammarhs
+      #SYNC: ML_lapserhs
+      #SYNC: ML_log_confacrhs
+      #SYNC: ML_metricrhs
+      #SYNC: ML_shiftrhs
+      #SYNC: ML_trace_curvrhs
+    } "Apply NewRad boundary conditions to RHS"
+  }
+  
+  
+  
+  SCHEDULE GROUP ML_CCZ4_convertToADMBaseGroup IN ML_CCZ4_convertToADMBaseGroupWrapper
+  {
+    #SYNC: ADMBase::metric
+    #SYNC: ADMBase::curv
+    #SYNC: ADMBase::lapse
+    #SYNC: ADMBase::shift
+    #SYNC: ADMBase::dtlapse
+    #SYNC: ADMBase::dtshift
+  } "Calculate ADM variables"
+  
+  if (CCTK_EQUALS (calculate_ADMBase_variables_at, "MoL_PostStep"))
+  {
+    SCHEDULE GROUP ML_CCZ4_convertToADMBaseGroupWrapper IN MoL_PostStep AFTER (ML_CCZ4_ApplyBCs ML_CCZ4_enforce) BEFORE (ADMBase_SetADMVars Whisky_PostStep)
+    {
+    } "Calculate ADM variables"
+  }
+  else if  (CCTK_EQUALS (calculate_ADMBase_variables_at, "CCTK_EVOL"))
+  {
+    SCHEDULE GROUP ML_CCZ4_convertToADMBaseGroupWrapper IN MoL_PseudoEvolution BEFORE (ADMBase_SetADMVars Whisky_PostStep)
+    {
+    } "Calculate ADM variables"
+  }
+  else if  (CCTK_EQUALS (calculate_ADMBase_variables_at, "CCTK_ANALYSIS"))
+  {
+    SCHEDULE GROUP ML_CCZ4_convertToADMBaseGroupWrapper AT post_recover_variables
+    {
+    } "Calculate ADM variables"
+    
+    SCHEDULE GROUP ML_CCZ4_convertToADMBaseGroupWrapper AT analysis BEFORE (ADMBase_SetADMVars Whisky_PostStep)
+    {
+    } "Calculate ADM variables"
+  }
+
+  SCHEDULE ML_CCZ4_SelectBCsADMBase IN ML_CCZ4_convertToADMBaseGroupWrapper AFTER ML_CCZ4_convertToADMBaseGroup
+  {
+    LANG: C
+    OPTIONS: level
+  } "Select boundary conditions for ADMBase variables"
+  
+  SCHEDULE GROUP ApplyBCs AS ML_CCZ4_ApplyBCsADMBase IN ML_CCZ4_convertToADMBaseGroupWrapper AFTER ML_CCZ4_SelectBCsADMBase
+  {
+  } "Apply boundary conditions to ADMBase variables"
+  
+}
diff --git a/ML_CCZ4_Helper/src/ExtrapolateGammas.c b/ML_CCZ4_Helper/src/ExtrapolateGammas.c
new file mode 100644
index 0000000..efa08dd
--- /dev/null
+++ b/ML_CCZ4_Helper/src/ExtrapolateGammas.c
@@ -0,0 +1,29 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+
+static void
+extrap (cGH const * restrict cctkGH,
+        CCTK_REAL * restrict var);
+
+void
+ML_CCZ4_ExtrapolateGammas (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  
+  extrap (cctkGH, Xt1);
+  extrap (cctkGH, Xt2);
+  extrap (cctkGH, Xt3);
+  
+  extrap (cctkGH, A);
+  
+  extrap (cctkGH, B1);
+  extrap (cctkGH, B2);
+  extrap (cctkGH, B3);
+}
+
+static void
+extrap (cGH const * restrict const cctkGH,
+        CCTK_REAL * restrict const var)
+{
+  ExtrapolateGammas (cctkGH, var);
+}
diff --git a/ML_CCZ4_Helper/src/NewRad.c b/ML_CCZ4_Helper/src/NewRad.c
new file mode 100644
index 0000000..756e5bf
--- /dev/null
+++ b/ML_CCZ4_Helper/src/NewRad.c
@@ -0,0 +1,69 @@
+#include <math.h>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+static void
+newrad (cGH const * restrict cctkGH,
+        CCTK_REAL const * restrict var,
+        CCTK_REAL * restrict rhs,
+        CCTK_REAL var0,
+        CCTK_REAL v0);
+
+void
+ML_CCZ4_NewRad (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_REAL const v0 = sqrt (harmonicF);
+  
+  newrad (cctkGH, phi  , phirhs  , 0.0, v0 );
+  
+  newrad (cctkGH, gt11 , gt11rhs , 1.0, 1.0);
+  newrad (cctkGH, gt12 , gt12rhs , 0.0, 1.0);
+  newrad (cctkGH, gt13 , gt13rhs , 0.0, 1.0);
+  newrad (cctkGH, gt22 , gt22rhs , 1.0, 1.0);
+  newrad (cctkGH, gt23 , gt23rhs , 0.0, 1.0);
+  newrad (cctkGH, gt33 , gt33rhs , 1.0, 1.0);
+  
+  newrad (cctkGH, Xt1  , Xt1rhs  , 0.0, 1.0);
+  newrad (cctkGH, Xt2  , Xt2rhs  , 0.0, 1.0);
+  newrad (cctkGH, Xt3  , Xt3rhs  , 0.0, 1.0);
+  
+  newrad (cctkGH, trK  , trKrhs  , 0.0, v0 );
+  
+  newrad (cctkGH, At11 , At11rhs , 0.0, 1.0);
+  newrad (cctkGH, At12 , At12rhs , 0.0, 1.0);
+  newrad (cctkGH, At13 , At13rhs , 0.0, 1.0);
+  newrad (cctkGH, At22 , At22rhs , 0.0, 1.0);
+  newrad (cctkGH, At23 , At23rhs , 0.0, 1.0);
+  newrad (cctkGH, At33 , At33rhs , 0.0, 1.0);
+  
+  newrad (cctkGH, alpha, alpharhs, 1.0, v0 );
+  
+  newrad (cctkGH, A    , Arhs    , 0.0, v0 );
+  
+  newrad (cctkGH, beta1, beta1rhs, 0.0, 1.0);
+  newrad (cctkGH, beta2, beta2rhs, 0.0, 1.0);
+  newrad (cctkGH, beta3, beta3rhs, 0.0, 1.0);
+  
+  newrad (cctkGH, B1   , B1rhs   , 0.0, 1.0);
+  newrad (cctkGH, B2   , B2rhs   , 0.0, 1.0);
+  newrad (cctkGH, B3   , B3rhs   , 0.0, 1.0);
+
+  newrad (cctkGH, Theta, Thetarhs, 0.0, 1.0);
+}
+
+static void
+newrad (cGH const * restrict const cctkGH,
+        CCTK_REAL const * restrict const var,
+        CCTK_REAL * restrict const rhs,
+        CCTK_REAL const var0,
+        CCTK_REAL const v0)
+{
+  DECLARE_CCTK_PARAMETERS;
+  
+  NewRad_Apply (cctkGH, var, rhs, var0, v0, radpower);
+}
diff --git a/ML_CCZ4_Helper/src/RegisterConstrained.c b/ML_CCZ4_Helper/src/RegisterConstrained.c
new file mode 100644
index 0000000..ab906a3
--- /dev/null
+++ b/ML_CCZ4_Helper/src/RegisterConstrained.c
@@ -0,0 +1,32 @@
+#include <assert.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+
+static void
+register_constrained (char const * restrict const gn);
+
+void
+ML_CCZ4_RegisterConstrained (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  
+  register_constrained ("ADMBase::metric");
+  register_constrained ("ADMBase::curv");
+  register_constrained ("ADMBase::lapse");
+  register_constrained ("ADMBase::shift");
+  register_constrained ("ADMBase::dtlapse");
+  register_constrained ("ADMBase::dtshift");
+}
+
+static void
+register_constrained (char const * restrict const gn)
+{
+  assert (gn);
+  
+  int const gi = CCTK_GroupIndex (gn);
+  int const ierr = MoLRegisterConstrainedGroup (gi);
+  assert (! ierr);
+}
diff --git a/ML_CCZ4_Helper/src/RegisterSlicing.c b/ML_CCZ4_Helper/src/RegisterSlicing.c
new file mode 100644
index 0000000..437e185
--- /dev/null
+++ b/ML_CCZ4_Helper/src/RegisterSlicing.c
@@ -0,0 +1,9 @@
+#include <cctk.h>
+#include <Slicing.h>
+
+int
+ML_CCZ4_RegisterSlicing (void)
+{
+  Einstein_RegisterSlicing ("ML_CCZ4");
+  return 0;
+}
diff --git a/ML_CCZ4_Helper/src/SelectBCsADMBase.c b/ML_CCZ4_Helper/src/SelectBCsADMBase.c
new file mode 100644
index 0000000..10c01f9
--- /dev/null
+++ b/ML_CCZ4_Helper/src/SelectBCsADMBase.c
@@ -0,0 +1,31 @@
+#include <assert.h>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+static void
+select_bcs (cGH const * restrict cctkGH, char const * restrict gn);
+
+void
+ML_CCZ4_SelectBCsADMBase (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  
+  select_bcs (cctkGH, "ADMBase::metric");
+  select_bcs (cctkGH, "ADMBase::curv");
+  select_bcs (cctkGH, "ADMBase::lapse");
+  select_bcs (cctkGH, "ADMBase::shift");
+  select_bcs (cctkGH, "ADMBase::dtlapse");
+  select_bcs (cctkGH, "ADMBase::dtshift");
+}
+
+static void
+select_bcs (cGH const * restrict const cctkGH, char const * restrict const gn)
+{
+  DECLARE_CCTK_PARAMETERS;
+  
+  int const ierr = Boundary_SelectGroupForBC
+    (cctkGH, CCTK_ALL_FACES, boundary_width, -1, gn, "none");
+  assert (! ierr);
+}
diff --git a/ML_CCZ4_Helper/src/SetGroupTags.c b/ML_CCZ4_Helper/src/SetGroupTags.c
new file mode 100644
index 0000000..00878ee
--- /dev/null
+++ b/ML_CCZ4_Helper/src/SetGroupTags.c
@@ -0,0 +1,77 @@
+#include <cctk.h>
+#include <cctk_Parameters.h>
+#include <util_Table.h>
+
+#include <assert.h>
+
+static void
+set_group_tags (int const checkpoint,
+                int const persistent,
+                int const prolongate,
+                char const * restrict const gn);
+
+int
+ML_CCZ4_SetGroupTags (void)
+{
+  DECLARE_CCTK_PARAMETERS;
+  
+  int const checkpoint = timelevels > 1;
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::metric");
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::curv");
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::lapse");
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::shift");
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::dtlapse");
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::dtshift");
+  
+  set_group_tags (checkpoint, checkpoint, 0, "ML_CCZ4::ML_cons_detg");
+  set_group_tags (checkpoint, checkpoint, 0, "ML_CCZ4::ML_cons_Gamma");
+  set_group_tags (checkpoint, checkpoint, 0, "ML_CCZ4::ML_cons_traceA");
+  set_group_tags (checkpoint, checkpoint, 0, "ML_CCZ4::ML_Ham");
+  set_group_tags (checkpoint, checkpoint, 0, "ML_CCZ4::ML_mom");
+  
+  int const rhs_checkpoint = rhs_timelevels > 1;
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_CCZ4::ML_log_confacrhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_CCZ4::ML_metricrhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_CCZ4::ML_Gammarhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_CCZ4::ML_trace_curvrhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_CCZ4::ML_curvrhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_CCZ4::ML_lapserhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_CCZ4::ML_dtlapserhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_CCZ4::ML_shiftrhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_CCZ4::ML_dtshiftrhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_CCZ4::ML_Thetarhs");
+  
+  return 0;
+}
+
+static void
+set_group_tags (int const checkpoint,
+                int const persistent,
+                int const prolongate,
+                char const * restrict const gn)
+{
+  assert (gn);
+  
+  int const gi = CCTK_GroupIndex (gn);
+  assert (gi >= 0);
+  
+  int const table = CCTK_GroupTagsTableI (gi);
+  assert (table >= 0);
+  
+  if (! checkpoint) {
+    int const ierr = Util_TableSetString (table, "no", "Checkpoint");
+    assert (! ierr);
+  }
+  
+  if (! persistent) {
+    int const ierr = Util_TableSetString (table, "no", "Persistent");
+    assert (! ierr);
+  }
+  
+  if (! prolongate) {
+    int const iret = Util_TableDeleteKey (table, "ProlongationParameter");
+    assert (iret == 0 || iret == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    int const ierr = Util_TableSetString (table, "none", "Prolongation");
+    assert (! ierr);
+  }
+}
diff --git a/ML_CCZ4_Helper/src/make.code.defn b/ML_CCZ4_Helper/src/make.code.defn
new file mode 100644
index 0000000..a3ac433
--- /dev/null
+++ b/ML_CCZ4_Helper/src/make.code.defn
@@ -0,0 +1,2 @@
+# -*-Makefile-*-
+SRCS = ExtrapolateGammas.c NewRad.c RegisterConstrained.c RegisterSlicing.c SelectBCsADMBase.c SetGroupTags.c