Initial commit.

6 files changed, 231 insertions, 0 deletions

diff --git a/configuration.ccl b/configuration.ccl
new file mode 100644
index 0000000..12661fa
--- /dev/null
+++ b/configuration.ccl
@@ -0,0 +1,2 @@
+# Configuration definition for thorn NullSurf
+
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..612f623
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,29 @@
+# Interface definition for thorn NullSurf
+implements: NullSurf
+
+INHERITS: ADMBase grid CoordBase MethodOfLines ML_BSSN Driver
+
+CCTK_INT FUNCTION MoLRegisterEvolved(CCTK_INT IN EvolvedIndex, CCTK_INT IN RHSIndex)
+
+REQUIRES FUNCTION MoLRegisterEvolved
+
+CCTK_POINTER FUNCTION                     \
+    VarDataPtrI                           \
+        (CCTK_POINTER_TO_CONST IN cctkGH, \
+         CCTK_INT IN map,                 \
+         CCTK_INT IN reflevel,            \
+         CCTK_INT IN component,           \
+         CCTK_INT IN timelevel,           \
+         CCTK_INT IN varindex)
+USES FUNCTION VarDataPtrI
+
+public:
+REAL null_surface     TYPE=GF TIMELEVELS=3 tags='tensortypealias="Scalar" tensorweight=0'
+{
+    F
+}
+
+REAL null_surface_rhs TYPE=GF TIMELEVELS=3 tags='tensortypealias="Scalar" tensorweight=0 Prolongation="None"'
+{
+    F_rhs
+}
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..0210523
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,12 @@
+# Parameter definitions for thorn NullSurf
+
+shares: MethodOfLines
+
+USES CCTK_INT MoL_Num_Evolved_Vars
+
+restricted:
+CCTK_INT mol_evolved "Number of evolved variables used by this thorn" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Evolved_Vars STEERABLE=RECOVER
+{
+  1:1 :: "Number of evolved variables used by this thorn"
+} 1
+
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..d6c38b3
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,39 @@
+# Schedule definitions for thorn NullSurf
+#
+SCHEDULE ns_initial IN ADMBase_InitialData {
+    LANG: C
+} "NullSurf initial data"
+
+SCHEDULE ns_eval_rhs IN MoL_CalcRHS {
+    LANG: C
+    READS: ML_BSSN::alpha(Interior)
+    READS: ML_BSSN::phi(Interior)
+    READS: ML_BSSN::gt11(Interior)
+    READS: ML_BSSN::gt12(Interior)
+    READS: ML_BSSN::gt13(Interior)
+    READS: ML_BSSN::gt22(Interior)
+    READS: ML_BSSN::gt23(Interior)
+    READS: ML_BSSN::gt33(Interior)
+    READS: ML_BSSN::beta1(Interior)
+    READS: ML_BSSN::beta2(Interior)
+    READS: ML_BSSN::beta3(Interior)
+    WRITES: NullSurf::F(Interior)
+} "NullSurf eval RHS"
+
+SCHEDULE ns_mol_register in MoL_Register {
+    LANG: C
+} "NullSurf register MoL variables"
+
+SCHEDULE ns_register_symmetries in SymmetryRegister {
+    LANG: C
+} "NullSurf register symmetry properties"
+
+schedule ns_select_bc in MoL_PostStep
+{
+  LANG: C
+  OPTIONS: level
+  SYNC: null_surface
+} "select boundary conditions"
+
+STORAGE: null_surface[3]
+STORAGE: null_surface_rhs[3]
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..8fa0281
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,7 @@
+# Main make.code.defn file for thorn MaximalSlicingAxi
+
+# Source files in this directory
+SRCS = nullsurf.c
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/nullsurf.c b/src/nullsurf.c
new file mode 100644
index 0000000..2b88437
--- /dev/null
+++ b/src/nullsurf.c
@@ -0,0 +1,142 @@
+
+#include <math.h>
+#include <stddef.h>
+#include <stdint.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "Symmetry.h"
+
+#define SQR(x) ((x) * (x))
+
+void ns_mol_register(CCTK_ARGUMENTS)
+{
+    MoLRegisterEvolved(CCTK_VarIndex("NullSurf::F"),  CCTK_VarIndex("NullSurf::F_rhs"));
+}
+
+void ns_initial(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    for (int j = 0; j < cctkGH->cctk_lsh[2]; j++) {
+        for (int i = 0; i < cctkGH->cctk_lsh[0]; i++) {
+            const int idx = CCTK_GFINDEX3D(cctkGH, i, 0, j);
+            const double r2 = SQR(x[idx]) + SQR(z[idx]);
+
+            F[idx] = r2;
+        }
+    }
+}
+
+#define FD4(arr, stride, h_inv) \
+        ((-1.0 * arr[2 * stride] + 8.0 * arr[1 * stride] - 8.0 * arr[-1 * stride] + 1.0 * arr[-2 * stride]) * h_inv / 12.0)
+#define FD24(arr, stride, h_inv) \
+        ((-1.0 * arr[2 * stride] + 16.0 * arr[1 * stride] - 30.0 * arr[0] + 16.0 * arr[-1 * stride] - 1.0 * arr[-2 * stride]) * SQR(h_inv) / 12.0)
+
+static double diff4_dx(const double *arr, ptrdiff_t stride, double h_inv)
+{
+    return FD4(arr, 1, h_inv);
+}
+static double diff4_dz(const double *arr, ptrdiff_t stride, double h_inv)
+{
+    return FD4(arr, stride, h_inv);
+}
+
+void ns_eval_rhs(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    const int ghosts_x = cctkGH->cctk_nghostzones[0];
+    const int ghosts_z = cctkGH->cctk_nghostzones[2];
+
+    const ptrdiff_t stride_z = CCTK_GFINDEX3D(cctkGH, 0, 0, 1);
+
+    const double dx = x[1] - x[0];
+    const double dz = z[stride_z] - z[0];
+    const double dx_inv = 1.0 / dx;
+    const double dz_inv = 1.0 / dz;
+
+    for (int idx_z = ghosts_z; idx_z < cctkGH->cctk_lsh[2] - ghosts_z; idx_z++) {
+        for (int idx_x = ghosts_x; idx_x < cctkGH->cctk_lsh[0] - ghosts_x; idx_x++) {
+            const int idx = CCTK_GFINDEX3D(cctkGH, idx_x, 0, idx_z);
+
+            const double dF[3] = {diff4_dx(&F[idx], 1, dx_inv), 0.0, diff4_dz(&F[idx], stride_z, dz_inv)};
+
+            const double alpha_val   = alpha[idx];
+            const double beta_val[3] = { beta1[idx], 0.0, beta3[idx] };
+            const double phi_val     = phi[idx];
+
+            const double gtxx = gt11[idx];
+            const double gtxy = gt12[idx];
+            const double gtxz = gt13[idx];
+            const double gtyy = gt22[idx];
+            const double gtyz = gt23[idx];
+            const double gtzz = gt33[idx];
+
+            const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz);
+
+            double gtu[3][3], g4u[4][4];
+            double shift_term;
+
+            // \tilde{γ}^{ij}
+            gtu[0][0] =  (gtyy * gtzz - SQR(gtyz)) / det;
+            gtu[1][1] =  (gtxx * gtzz - SQR(gtxz)) / det;
+            gtu[2][2] =  (gtxx * gtyy - SQR(gtxy)) / det;
+            gtu[0][1] = -(gtxy * gtzz - gtyz * gtxz) / det;
+            gtu[0][2] =  (gtxy * gtyz - gtyy * gtxz) / det;
+            gtu[1][2] = -(gtxx * gtyz - gtxy * gtxz) / det;
+            gtu[1][0] = gtu[0][1];
+            gtu[2][0] = gtu[0][2];
+            gtu[2][1] = gtu[1][2];
+
+            g4u[0][0] = - 1.0 / SQR(alpha_val);
+            for (int i = 0; i < 3; i++) {
+                const double val = beta_val[i] / SQR(alpha_val);
+                g4u[1 + i][0] = val;
+                g4u[0][1 + i] = val;
+            }
+
+            for (int i = 0; i < 3; i++)
+                for (int j = 0; j < 3; j++) {
+                    g4u[1 + i][1 + j] = SQR(phi_val) * gtu[i][j] - beta_val[i] * beta_val[j] / SQR(alpha_val);
+                }
+
+            shift_term = 0.0;
+            for (int i = 0; i < 3; i++)
+                shift_term += g4u[0][1 + i] * dF[i];
+
+            double dF2 = 0.0;
+            for (int i = 0; i < 3; i++)
+                for (int j = 0; j < 3; j++)
+                    dF2 += g4u[1 + i][1 + j] * dF[i] * dF[j];
+
+            F_rhs[idx] = (-shift_term + sqrt(SQR(shift_term) - g4u[0][0] * dF2)) / g4u[0][0];
+        }
+    }
+}
+
+void ns_register_symmetries(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  int sym[3] = { 1, 1, 1 };
+  int ret;
+
+  ret = SetCartSymVN(cctkGH, sym, "NullSurf::F");
+  if (ret != 0)
+      CCTK_WARN(0, "Error registering symmetries");
+}
+
+void ns_select_bc(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  int ret = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, "NullSurf::F", "none");
+  if (ret != 0)
+      CCTK_WARN(0, "Error registering boundaries");
+}