Add runtime-configurable support for CCZ4.

6 files changed, 52 insertions, 22 deletions

diff --git a/param.ccl b/param.ccl
index 49169aa..35fa495 100644
--- a/param.ccl
+++ b/param.ccl
@@ -39,3 +39,8 @@ CCTK_REAL switchoff_time "" STEERABLE=recover
 BOOLEAN export_coeffs "Export the coefficients of the spectral expansion in alpha_coeffs" STEERABLE=recover
 {
 } "no"
+
+CCTK_INT ccz4 "Use CCZ4 or BSSN" STEERABLE=recover
+{
+    0:1 :: ""
+} 0
diff --git a/schedule.ccl b/schedule.ccl
index 3923b92..b4fcda2 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -3,6 +3,9 @@
 SCHEDULE quasimaximal_slicing_axi_eval IN ML_BSSN_evolCalcGroup BEFORE ML_BSSN_lapse_evol {
     LANG: C
 } "Quasimaximal slicing eval W"
+SCHEDULE quasimaximal_slicing_axi_eval IN ML_CCZ4_evolCalcGroup BEFORE ML_CCZ4_lapse_evol {
+    LANG: C
+} "Quasimaximal slicing eval W"
 
 #SCHEDULE quasimaximal_slicing_axi_solve IN ML_BSSN_evolCalcGroup BEFORE quasimaximal_slicing_axi_eval {
 #SCHEDULE quasimaximal_slicing_axi_solve IN MoL_PostStep AFTER ML_BSSN_ApplyBCs {
diff --git a/src/common.h b/src/common.h
index 44b0674..6eb579f 100644
--- a/src/common.h
+++ b/src/common.h
@@ -4,7 +4,6 @@
 #define HAVE_OPENCL 0
 #define QMS_VERIFY  0
 #define QMS_POLAR   0
-#define QMS_CCZ4    0
 
 #define SQR(x) ((x) * (x))
 #define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
diff --git a/src/qms.c b/src/qms.c
index 90e1f76..6faa62e 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -232,7 +232,7 @@ static int context_init(cGH *cctkGH)
     qms->max_radius = 60.0;
 
     ret = qms_solver_init(&qms->solver, cctkGH, basis_order_r, basis_order_z,
-                          outer_bound, filter_power, 0.0);
+                          outer_bound, filter_power, 0.0, ccz4);
     if (ret < 0)
         return ret;
 
@@ -460,16 +460,30 @@ void qms_init(CCTK_ARGUMENTS)
     if (!qms_context)
         context_init(cctkGH);
 
-    double *Kdot11 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot11");
-    double *Kdot22 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot22");
-    double *Kdot33 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot33");
-    double *Kdot12 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot12");
-    double *Kdot13 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot13");
-    double *Kdot23 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot23");
-
-    double *Xtdot1 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Xtdot1");
-    double *Xtdot2 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Xtdot2");
-    double *Xtdot3 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Xtdot3");
+    double *Kdot11, *Kdot22, *Kdot33, *Kdot12, *Kdot13, *Kdot23, *Xtdot1, *Xtdot2, *Xtdot3;
+    if (ccz4) {
+        Kdot11 = CCTK_VarDataPtr(cctkGH, 0, "ML_CCZ4::Kdot11");
+        Kdot22 = CCTK_VarDataPtr(cctkGH, 0, "ML_CCZ4::Kdot22");
+        Kdot33 = CCTK_VarDataPtr(cctkGH, 0, "ML_CCZ4::Kdot33");
+        Kdot12 = CCTK_VarDataPtr(cctkGH, 0, "ML_CCZ4::Kdot12");
+        Kdot13 = CCTK_VarDataPtr(cctkGH, 0, "ML_CCZ4::Kdot13");
+        Kdot23 = CCTK_VarDataPtr(cctkGH, 0, "ML_CCZ4::Kdot23");
+
+        Xtdot1 = CCTK_VarDataPtr(cctkGH, 0, "ML_CCZ4::Xtdot1");
+        Xtdot2 = CCTK_VarDataPtr(cctkGH, 0, "ML_CCZ4::Xtdot2");
+        Xtdot3 = CCTK_VarDataPtr(cctkGH, 0, "ML_CCZ4::Xtdot3");
+    } else {
+        Kdot11 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot11");
+        Kdot22 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot22");
+        Kdot33 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot33");
+        Kdot12 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot12");
+        Kdot13 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot13");
+        Kdot23 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Kdot23");
+
+        Xtdot1 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Xtdot1");
+        Xtdot2 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Xtdot2");
+        Xtdot3 = CCTK_VarDataPtr(cctkGH, 0, "ML_BSSN::Xtdot3");
+    }
 
     for (int k = 0; k < cctk_lsh[2]; k++)
         for (int j = 0; j < cctk_lsh[1]; j++)
diff --git a/src/qms_solve.c b/src/qms_solve.c
index f1e6ee8..3ce463b 100644
--- a/src/qms_solve.c
+++ b/src/qms_solve.c
@@ -188,8 +188,7 @@ struct QMSSolverPriv {
 };
 
 /* mapping between our indices and thorn names */
-static const char *metric_vars[] = {
-#if QMS_CCZ4
+static const char *metric_vars_ccz4[] = {
     [GTXX]  = "ML_CCZ4::gt11",
     [GTYY]  = "ML_CCZ4::gt22",
     [GTZZ]  = "ML_CCZ4::gt33",
@@ -214,7 +213,8 @@ static const char *metric_vars[] = {
     [KDOT_XY] = "ML_CCZ4::Kdot12",
     [KDOT_XZ] = "ML_CCZ4::Kdot13",
     [KDOT_YZ] = "ML_CCZ4::Kdot23",
-#else
+};
+static const char *metric_vars_bssn[] = {
     [GTXX]  = "ML_BSSN::gt11",
     [GTYY]  = "ML_BSSN::gt22",
     [GTZZ]  = "ML_BSSN::gt33",
@@ -240,7 +240,6 @@ static const char *metric_vars[] = {
     [KDOT_XY] = "ML_BSSN::Kdot12",
     [KDOT_XZ] = "ML_BSSN::Kdot13",
     [KDOT_YZ] = "ML_BSSN::Kdot23",
-#endif
 };
 
 /* mapping between the cactus grid values and interpolated values */
@@ -1093,7 +1092,8 @@ fail:
 int qms_solver_init(QMSSolver **pctx,
                     cGH *cctkGH,
                     int basis_order_r, int basis_order_z,
-                    double outer_bound, double filter_power, double input_filter_power)
+                    double outer_bound, double filter_power, double input_filter_power,
+                    int ccz4)
 {
     QMSSolver *ctx;
     QMSSolverPriv *s;
@@ -1219,10 +1219,18 @@ int qms_solver_init(QMSSolver **pctx,
         s->interp_value_codes[i] = CCTK_VARIABLE_REAL;
     }
 
-    for (int i = 0; i < ARRAY_ELEMS(metric_vars); i++) {
-        s->interp_vars_indices[i] = CCTK_VarIndex(metric_vars[i]);
-        if (s->interp_vars_indices[i] < 0)
-            CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, "Error getting the index of variable: %s\n", metric_vars[i]);
+    if (ccz4) {
+        for (int i = 0; i < ARRAY_ELEMS(metric_vars_ccz4); i++) {
+            s->interp_vars_indices[i] = CCTK_VarIndex(metric_vars_ccz4[i]);
+            if (s->interp_vars_indices[i] < 0)
+                CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, "Error getting the index of variable: %s\n", metric_vars_ccz4[i]);
+        }
+    } else {
+        for (int i = 0; i < ARRAY_ELEMS(metric_vars_bssn); i++) {
+            s->interp_vars_indices[i] = CCTK_VarIndex(metric_vars_bssn[i]);
+            if (s->interp_vars_indices[i] < 0)
+                CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, "Error getting the index of variable: %s\n", metric_vars_bssn[i]);
+        }
     }
 
     s->coord_system = CCTK_CoordSystemHandle("cart3d");
diff --git a/src/qms_solve.h b/src/qms_solve.h
index 3e9ca49..91520d9 100644
--- a/src/qms_solve.h
+++ b/src/qms_solve.h
@@ -41,7 +41,8 @@ typedef struct QMSSolver {
 int qms_solver_init(QMSSolver **ctx,
                     cGH *cctkGH,
                     int basis_order_r, int basis_order_z,
-                    double outer_bound, double filter_power, double input_filter_power);
+                    double outer_bound, double filter_power, double input_filter_power,
+                    int ccz4);
 
 void qms_solver_free(QMSSolver **ctx);