Integrate individual photons rather than null surfaces.

4 files changed, 691 insertions, 143 deletions

diff --git a/interface.ccl b/interface.ccl
index 612f623..b0f1325 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -1,11 +1,7 @@
 # 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
+INHERITS: ADMBase grid CoordBase ML_BSSN Driver
 
 CCTK_POINTER FUNCTION                     \
     VarDataPtrI                           \
@@ -17,13 +13,26 @@ CCTK_POINTER FUNCTION                     \
          CCTK_INT IN varindex)
 USES FUNCTION VarDataPtrI
 
+
 public:
-REAL null_surface     TYPE=GF TIMELEVELS=3 tags='tensortypealias="Scalar" tensorweight=0'
+REAL photon_coord TYPE=array DIM=2 SIZE=nb_surfaces,photons_per_surface DISTRIB=constant
 {
-    F
+    pos_x,
+    # stored only to send zeroes to interpolator
+    pos_y,
+    pos_z,
+    mom_x,
+    mom_z
 }
 
-REAL null_surface_rhs TYPE=GF TIMELEVELS=3 tags='tensortypealias="Scalar" tensorweight=0 Prolongation="None"'
+REAL photon_rhs TYPE=array DIM=3 SIZE=4,nb_surfaces,photons_per_surface DISTRIB=constant
 {
-    F_rhs
+    pos_x_rhs,
+    pos_z_rhs,
+    mom_x_rhs,
+    mom_z_rhs
 }
+
+REAL photon_times TYPE=array DIM=1 SIZE=nb_surfaces DISTRIB=constant
+INT evol_next_step TYPE=scalar
+INT surfaces_active TYPE=scalar
diff --git a/param.ccl b/param.ccl
index 0210523..af9189c 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,12 +1,16 @@
 # Parameter definitions for thorn NullSurf
 
-shares: MethodOfLines
-
-USES CCTK_INT MoL_Num_Evolved_Vars
+CCTK_INT photons_per_surface "Number of photons per surface"
+{
+    2: :: ""
+} 64
 
-restricted:
-CCTK_INT mol_evolved "Number of evolved variables used by this thorn" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Evolved_Vars STEERABLE=RECOVER
+CCTK_INT nb_surfaces "Number of null surfaces to trace"
 {
-  1:1 :: "Number of evolved variables used by this thorn"
+    1: :: ""
 } 1
 
+CCTK_INT solve_level "Refinement level to integrate photons on"
+{
+    0: :: ""
+} 0
diff --git a/schedule.ccl b/schedule.ccl
index d6c38b3..e743105 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,39 +1,22 @@
 # 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"
+if (nb_surfaces > 0) {
+    SCHEDULE ns_evol IN CCTK_ANALYSIS {
+        LANG: C
+    } "NullSurf eval RHS"
 
-SCHEDULE ns_mol_register in MoL_Register {
-    LANG: C
-} "NullSurf register MoL variables"
+    SCHEDULE ns_evol_init IN CCTK_BASEGRID AFTER TemporalSpacings {
+        LANG: C
+    } ""
 
-SCHEDULE ns_register_symmetries in SymmetryRegister {
-    LANG: C
-} "NullSurf register symmetry properties"
+    SCHEDULE ns_init IN CCTK_INITIAL {
+        LANG: C
+    } ""
 
-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]
+    STORAGE: photon_coord
+    STORAGE: photon_rhs
+    STORAGE: photon_times
+    STORAGE: evol_next_step
+    STORAGE: surfaces_active
+}
diff --git a/src/nullsurf.c b/src/nullsurf.c
index 2b88437..755be14 100644
--- a/src/nullsurf.c
+++ b/src/nullsurf.c
@@ -1,142 +1,694 @@
 
+#include <errno.h>
+#include <float.h>
+#include <limits.h>
 #include <math.h>
 #include <stddef.h>
 #include <stdint.h>
+#include <stdlib.h>
+#include <string.h>
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
+#include "util_Table.h"
 #include "Symmetry.h"
 
 #define SQR(x) ((x) * (x))
+#define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
+#define MAX(x, y) ((x) > (y) ? (x) : (y))
+#define MIN(x, y) ((x) > (y) ? (y) : (x))
+#define ABS(x) ((x >= 0) ? (x) : -(x))
+#define ARRAY_ELEMS(arr) (sizeof(arr) / sizeof(*arr))
 
-void ns_mol_register(CCTK_ARGUMENTS)
+/* indices (in our code, not cactus structs) of the grid functions which we'll need to
+ * interpolate on the photon positions */
+enum MetricVars {
+    GTXX = 0,
+    GTYY,
+    GTZZ,
+    GTXY,
+    GTXZ,
+    GTYZ,
+    ATXX,
+    ATYY,
+    ATZZ,
+    ATXY,
+    ATXZ,
+    ATYZ,
+    PHI,
+    ALPHA,
+    TRK,
+    NB_METRIC_VARS,
+};
+
+/* indices of the interpolated values of the above grid functions and their derivatives */
+enum InterpMetricVars {
+    I_GTXX = 0,
+    I_GTYY,
+    I_GTZZ,
+    I_GTXY,
+    I_GTXZ,
+    I_GTYZ,
+    I_GTXX_DX,
+    I_GTYY_DX,
+    I_GTZZ_DX,
+    I_GTXZ_DX,
+    I_GTXX_DZ,
+    I_GTYY_DZ,
+    I_GTZZ_DZ,
+    I_GTXZ_DZ,
+    I_ATXX,
+    I_ATYY,
+    I_ATZZ,
+    I_ATXY,
+    I_ATXZ,
+    I_ATYZ,
+    I_PHI,
+    I_PHI_DX,
+    I_PHI_DZ,
+    I_ALPHA,
+    I_ALPHA_DX,
+    I_ALPHA_DZ,
+    I_TRK,
+    NB_INTERP_VARS,
+};
+
+typedef struct NSContext {
+    cGH *gh;
+
+    int reflevel;
+    double dt;
+
+    int nb_surfaces;
+    int photons_per_surface;
+
+    // interpolation parameters
+    int coord_system;
+    int interp_operator;
+    int interp_params;
+
+    int        interp_vars_indices[NB_METRIC_VARS];
+    double    *interp_values[NB_INTERP_VARS];
+    int        interp_value_codes[NB_INTERP_VARS];
+
+    double    *rk_scratch[4];
+} NSContext;
+
+static NSContext ns_ctx;
+
+/* mapping between our indices and thorn names */
+static const char *metric_vars[] = {
+    [GTXX]  = "ML_BSSN::gt11",
+    [GTYY]  = "ML_BSSN::gt22",
+    [GTZZ]  = "ML_BSSN::gt33",
+    [GTXY]  = "ML_BSSN::gt12",
+    [GTXZ]  = "ML_BSSN::gt13",
+    [GTYZ]  = "ML_BSSN::gt23",
+    [ATXX]  = "ML_BSSN::At11",
+    [ATYY]  = "ML_BSSN::At22",
+    [ATZZ]  = "ML_BSSN::At33",
+    [ATXY]  = "ML_BSSN::At12",
+    [ATXZ]  = "ML_BSSN::At13",
+    [ATYZ]  = "ML_BSSN::At23",
+    [PHI]   = "ML_BSSN::phi",
+    [ALPHA] = "ML_BSSN::alpha",
+    [TRK]   = "ML_BSSN::trK",
+};
+
+/* mapping between the cactus grid values and interpolated values */
+static const int interp_operation_indices[] = {
+    [I_GTXX]     = GTXX,
+    [I_GTYY]     = GTYY,
+    [I_GTZZ]     = GTZZ,
+    [I_GTXY]     = GTXY,
+    [I_GTXZ]     = GTXZ,
+    [I_GTYZ]     = GTYZ,
+    [I_GTXX_DX]  = GTXX,
+    [I_GTYY_DX]  = GTYY,
+    [I_GTZZ_DX]  = GTZZ,
+    [I_GTXZ_DX]  = GTXZ,
+    [I_GTXX_DZ]  = GTXX,
+    [I_GTYY_DZ]  = GTYY,
+    [I_GTZZ_DZ]  = GTZZ,
+    [I_GTXZ_DZ]  = GTXZ,
+    [I_ATXX]     = ATXX,
+    [I_ATYY]     = ATYY,
+    [I_ATZZ]     = ATZZ,
+    [I_ATXY]     = ATXY,
+    [I_ATXZ]     = ATXZ,
+    [I_ATYZ]     = ATYZ,
+    [I_PHI]      = PHI,
+    [I_PHI_DX]   = PHI,
+    [I_PHI_DZ]   = PHI,
+    [I_ALPHA]    = ALPHA,
+    [I_ALPHA_DX] = ALPHA,
+    [I_ALPHA_DZ] = ALPHA,
+    [I_TRK]      = TRK,
+};
+
+/* the operation (plain value or x/y/z-derivative) to apply during interpolation */
+static const CCTK_INT interp_operation_codes[] = {
+    [I_GTXX]     = 0,
+    [I_GTYY]     = 0,
+    [I_GTZZ]     = 0,
+    [I_GTXY]     = 0,
+    [I_GTXZ]     = 0,
+    [I_GTYZ]     = 0,
+    [I_GTXX_DX]  = 1,
+    [I_GTYY_DX]  = 1,
+    [I_GTZZ_DX]  = 1,
+    [I_GTXZ_DX]  = 1,
+    [I_GTXX_DZ]  = 3,
+    [I_GTYY_DZ]  = 3,
+    [I_GTZZ_DZ]  = 3,
+    [I_GTXZ_DZ]  = 3,
+    [I_ATXX]     = 0,
+    [I_ATYY]     = 0,
+    [I_ATZZ]     = 0,
+    [I_ATXY]     = 0,
+    [I_ATXZ]     = 0,
+    [I_ATYZ]     = 0,
+    [I_PHI]      = 0,
+    [I_PHI_DX]   = 1,
+    [I_PHI_DZ]   = 3,
+    [I_ALPHA]    = 0,
+    [I_ALPHA_DX] = 1,
+    [I_ALPHA_DZ] = 3,
+    [I_TRK]      = 0,
+};
+
+static int ctz(int a)
+{
+    int ret = 0;
+
+    if (!a)
+        return INT_MAX;
+
+    while (!(a & 1)) {
+        a >>= 1;
+        ret++;
+    }
+
+    return ret;
+}
+
+static void surface_launch(NSContext *ns, double *pos_x, double *pos_z, double *mom_x, double *mom_z)
 {
-    MoLRegisterEvolved(CCTK_VarIndex("NullSurf::F"),  CCTK_VarIndex("NullSurf::F_rhs"));
+    const int origin_idx = CCTK_GFINDEX3D(ns->gh, ns->gh->cctk_nghostzones[0], 0, ns->gh->cctk_nghostzones[2]);
+
+    double gtxx  = ((double*)CCTK_VarDataPtr(ns->gh, 0, metric_vars[GTXX]))[origin_idx];
+    double gtyy  = ((double*)CCTK_VarDataPtr(ns->gh, 0, metric_vars[GTYY]))[origin_idx];
+    double gtzz  = ((double*)CCTK_VarDataPtr(ns->gh, 0, metric_vars[GTZZ]))[origin_idx];
+    double gtxy  = ((double*)CCTK_VarDataPtr(ns->gh, 0, metric_vars[GTXY]))[origin_idx];
+    double gtxz  = ((double*)CCTK_VarDataPtr(ns->gh, 0, metric_vars[GTXZ]))[origin_idx];
+    double gtyz  = ((double*)CCTK_VarDataPtr(ns->gh, 0, metric_vars[GTYZ]))[origin_idx];
+    double phi   = ((double*)CCTK_VarDataPtr(ns->gh, 0, metric_vars[PHI]))[origin_idx];
+    double alpha = ((double*)CCTK_VarDataPtr(ns->gh, 0, metric_vars[ALPHA]))[origin_idx];
+
+    const double g[3][3] = { { gtxx / SQR(phi), gtxy / SQR(phi), gtxz / SQR(phi) },
+                             { gtxy / SQR(phi), gtyy / SQR(phi), gtyz / SQR(phi) },
+                             { gtxz / SQR(phi), gtyz / SQR(phi), gtzz / SQR(phi) }};  
+
+    /**
+     * send out photons along geodesics evenly spaced in angles θ between x- and z-axis
+     * p^μ = { 1, p^x, 0, p^z }
+     * g_μν p^μ p^ν = 0
+     * p^x = p^r sinθ
+     * p^z = p^r cosθ
+     */
+#pragma omp parallel for
+    for (int i = 0; i < ns->photons_per_surface; i++) {
+        const double theta = (double)i / (ns->photons_per_surface - 1) * M_PI / 2;
+        const double st = sin(theta);
+        const double ct = cos(theta);
+
+        const double a = g[0][0] * SQR(st) + g[2][2] * SQR(ct) + 2.0 * g[0][2] * st * ct;
+        const double b = 0.0; // assume shift is zero
+        const double c = -SQR(alpha);
+
+        const double p_r = (-b + sqrt(SQR(b) - 4 * a * c)) / (2 * a);
+
+        pos_x[i] = 0.0;
+        pos_z[i] = 0.0;
+
+        mom_x[i] = st * p_r;
+        mom_z[i] = ct * p_r;
+    }
 }
 
-void ns_initial(CCTK_ARGUMENTS)
+static void null_geodesic_rhs(size_t len,
+                              double *dpos_dt_x, double *dpos_dt_z,
+                              double *dmom_dt_x, double *dmom_dt_z,
+                              const double *pos_x, const double *pos_z,
+                              const double *mom_x, const double *mom_z,
+                              const double * const metric_vars[NB_INTERP_VARS])
 {
-    DECLARE_CCTK_ARGUMENTS;
-    DECLARE_CCTK_PARAMETERS;
+#pragma omp parallel for
+    for (int idx = 0; idx < len; idx++) {
+        const double x = pos_x[idx];
+
+        const int zaxis = x <= 1e-12;
+
+        const double gtxx = metric_vars[I_GTXX][idx];
+        const double gtyy = metric_vars[I_GTYY][idx];
+        const double gtzz = metric_vars[I_GTZZ][idx];
+        const double gtxy = metric_vars[I_GTXY][idx];
+        const double gtxz = metric_vars[I_GTXZ][idx];
+        const double gtyz = metric_vars[I_GTYZ][idx];
+
+        const double gt[3][3] = {{ gtxx, gtxy, gtxz },
+                                 { gtxy, gtyy, gtyz },
+                                 { gtxz, gtyz, gtzz }};
 
-    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]);
+        const double dx_gt11 = metric_vars[I_GTXX_DX][idx];
+        const double dx_gt22 = metric_vars[I_GTYY_DX][idx];
+        const double dx_gt33 = metric_vars[I_GTZZ_DX][idx];
+        const double dx_gt13 = metric_vars[I_GTXZ_DX][idx];
 
-            F[idx] = r2;
+        const double dz_gt11 = metric_vars[I_GTXX_DZ][idx];
+        const double dz_gt22 = metric_vars[I_GTYY_DZ][idx];
+        const double dz_gt33 = metric_vars[I_GTZZ_DZ][idx];
+        const double dz_gt13 = metric_vars[I_GTXZ_DZ][idx];
+
+        const double dgt[3][3][3] = {
+            {
+                { dx_gt11,     0.0, dx_gt13 },
+                {     0.0, dx_gt22,     0.0 },
+                { dx_gt13,     0.0, dx_gt33 },
+            },
+            {
+                {     0.0, zaxis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0 },
+                { zaxis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0, zaxis ? dx_gt13 : gtxz / x },
+                { 0.0, zaxis ? dx_gt13 : gtxz / x, 0.0 },
+            },
+            {
+                { dz_gt11,     0.0, dz_gt13 },
+                {     0.0, dz_gt22,     0.0 },
+                { dz_gt13,     0.0, dz_gt33 },
+            },
+        };
+
+        const double Atxx = metric_vars[I_ATXX][idx];
+        const double Atyy = metric_vars[I_ATYY][idx];
+        const double Atzz = metric_vars[I_ATZZ][idx];
+        const double Atxy = metric_vars[I_ATXY][idx];
+        const double Atxz = metric_vars[I_ATXZ][idx];
+        const double Atyz = metric_vars[I_ATYZ][idx];
+
+        const double At[3][3] = {{ Atxx, Atxy, Atxz },
+                                 { Atxy, Atyy, Atyz },
+                                 { Atxz, Atyz, Atzz }};
+
+        const double trK    = metric_vars[I_TRK][idx];
+
+        const double alpha     = metric_vars[I_ALPHA][idx];
+        const double dx_alpha  = metric_vars[I_ALPHA_DX][idx];
+        const double dz_alpha  = metric_vars[I_ALPHA_DZ][idx];
+
+        const double phi    = metric_vars[I_PHI][idx];
+        const double dx_phi = metric_vars[I_PHI_DX][idx];
+        const double dz_phi = metric_vars[I_PHI_DZ][idx];
+
+        const double dphi[3] = { dx_phi, 0.0, dz_phi };
+
+        double g[3][3], gu[3][3], gtu[3][3], dg[3][3][3];
+        double K[3][3];
+
+        double g4[4][4], gu4[4][4], dg4[4][4][4];
+        double G4_x[4][4], G4_z[4][4];
+        double mom4[4];
+
+        double tmp, shift_term;
+
+        const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz);
+
+        // \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];
+
+        // γ_{jk}/^{jk}
+        // K_{ij}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++) {
+                gu[j][k] = SQR(phi) * gtu[j][k];
+                g[j][k]  = gt[j][k] / SQR(phi);
+                K[j][k]  = At[j][k] / SQR(phi) + (1.0 / 3.0) * g[j][k] * trK;
+            }
+
+        // ∂_j γ_{kl}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++) {
+                    dg[j][k][l] = -2.0 * dphi[j] * gt[k][l] / (phi * SQR(phi)) + dgt[j][k][l] / SQR(phi);
+                }
+
+        // 4-metric
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++)
+                g4[1 + i][1 + j] = g[i][j];
+        for (int i = 0; i < 3; i++) {
+            g4[0][1 + i] = 0.0; // assume zero shift
+            g4[1 + i][0] = 0.0;
+        }
+        g4[0][0] = -SQR(alpha);
+
+        /* derivative of the 4-metric */
+        /* spatial part */
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++) {
+                for (int k = 0; k < 3; k++)
+                    dg4[1 + k][1 + i][1 + j] = dg[k][i][j];
+                dg4[0][1 + i][1 + j] = -2.0 * alpha * K[i][j]; // + Dbeta[i][j] + Dbeta[j][i]
+            }
+
+        /* temporal part, assume shift is zero */
+        for (int i = 0; i < 4; i++)
+            for (int j = 0; j < 3; j++) {
+                dg4[i][1 + j][0] = 0.0;
+                dg4[i][0][1 + j] = 0.0;
+            }
+        dg4[1][0][0] = -2 * alpha * dx_alpha;
+        dg4[2][0][0] = 0.0;
+        dg4[3][0][0] = -2 * alpha * dz_alpha;
+        dg4[0][0][0] = 1e6; // should only be used multiplied by zero
+
+        /* inverse 4-metric */
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++)
+                gu4[1 + i][1 + j] = gu[i][j];
+
+        for (int i = 0; i < 3; i++) {
+            gu4[1 + i][0] = 0.0;
+            gu4[0][1 + i] = 0.0;
         }
+        gu4[0][0] = -1. / SQR(alpha);
+
+        for (int i = 0; i < 4; i++)
+            for (int j = 0; j < 4; j++) {
+                double val_x = 0.0;
+                double val_z = 0.0;
+                for (int k = 0; k < 4; k++) {
+                    val_x += gu4[1][k] * (dg4[i][k][j] + dg4[j][k][i] - dg4[k][i][j]);
+                    val_z += gu4[3][k] * (dg4[i][k][j] + dg4[j][k][i] - dg4[k][i][j]);
+                }
+                G4_x[i][j] = 0.5 * val_x;
+                G4_z[i][j] = 0.5 * val_z;
+            }
+
+        mom4[1] = mom_x[idx];
+        mom4[2] = 0.0;
+        mom4[3] = mom_z[idx];
+
+        shift_term = 0.0;
+        for (int i = 0; i < 3; i++)
+            shift_term += g4[0][1 + i] * mom4[1 + i];
+
+        tmp = 0.0;
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++)
+                tmp += g4[1 + i][1 + j] * mom4[1 + i] * mom4[1 + j];
+        mom4[0] = 1.0 / g4[0][0] * (-shift_term - sqrt(SQR(shift_term) - g4[0][0] * tmp));
+
+        tmp = 0.0;
+        for (int i = 0; i < 4; i++)
+            for (int j = 0; j < 4; j++)
+                tmp += g4[i][j] * mom4[i] * mom4[j];
+        if (fabs(tmp) > 1e-10)
+            abort();
+
+        dpos_dt_x[idx] = mom_x[idx] / mom4[0];
+        dpos_dt_z[idx] = mom_z[idx] / mom4[0];
+
+        tmp = 0.0;
+        for (int i = 0; i < 4; i++)
+            for (int j = 0; j < 4; j++)
+                tmp += G4_x[i][j] * mom4[i] * mom4[j];
+        dmom_dt_x[idx] = -tmp / mom4[0];
+
+        tmp = 0.0;
+        for (int i = 0; i < 4; i++)
+            for (int j = 0; j < 4; j++)
+                tmp += G4_z[i][j] * mom4[i] * mom4[j];
+        dmom_dt_z[idx] = -tmp / mom4[0];
     }
 }
 
-#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 void interp_geometry(NSContext *ns, int surfaces_active,
+                            const double *pos_x, const double *pos_y, const double *pos_z)
+{
+    const void * const coords[] = { pos_x, pos_y, pos_z };
+    int ret;
+
+    ret = CCTK_InterpGridArrays(ns->gh, 3, ns->interp_operator, ns->interp_params,
+                                ns->coord_system, surfaces_active * ns->photons_per_surface,
+                                CCTK_VARIABLE_REAL, coords,
+                                ARRAY_ELEMS(ns->interp_vars_indices), ns->interp_vars_indices,
+                                ARRAY_ELEMS(ns->interp_values), ns->interp_value_codes,
+                                (void * const *)ns->interp_values);
+    if (ret < 0)
+        CCTK_WARN(0, "Error interpolating");
+}
 
-static double diff4_dx(const double *arr, ptrdiff_t stride, double h_inv)
+static void
+ns_evol_rk_substep(NSContext *ns, int substep, int surfaces_active)
 {
-    return FD4(arr, 1, h_inv);
+    const ptrdiff_t rk_rhs_stride = ns->nb_surfaces * ns->photons_per_surface;
+
+    double *pos_x_rhs = (double*)CCTK_VarDataPtr(ns->gh, 0, "NullSurf::pos_x_rhs") + rk_rhs_stride * substep;
+    double *pos_z_rhs = (double*)CCTK_VarDataPtr(ns->gh, 0, "NullSurf::pos_z_rhs") + rk_rhs_stride * substep;
+    double *mom_x_rhs = (double*)CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_x_rhs") + rk_rhs_stride * substep;
+    double *mom_z_rhs = (double*)CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_z_rhs") + rk_rhs_stride * substep;
+
+    double *pos_x = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::pos_x");
+    double *pos_y = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::pos_y");
+    double *pos_z = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::pos_z");
+    double *mom_x = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_x");
+    double *mom_z = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_z");
+
+    double *step_pos_x, *step_pos_z, *step_mom_x, *step_mom_z;
+
+    /* compute the coordinate/momenta to evaluate the RHS for */
+    switch (substep) {
+        case 0:
+            step_pos_x = pos_x;
+            step_pos_z = pos_z;
+            step_mom_x = mom_x;
+            step_mom_z = mom_z;
+            break;
+        case 1:
+        case 2:
+        case 3: {
+            double fact = ns->dt * (substep == 3 ? 1.0 : 0.5);
+
+            step_pos_x = ns->rk_scratch[0];
+            step_pos_z = ns->rk_scratch[1];
+            step_mom_x = ns->rk_scratch[2];
+            step_mom_z = ns->rk_scratch[3];
+
+#pragma omp parallel for
+            for (int i = 0; i < surfaces_active * ns->photons_per_surface; i++) {
+                step_pos_x[i] = pos_x[i] + fact * pos_x_rhs[i - rk_rhs_stride];
+                step_pos_z[i] = pos_z[i] + fact * pos_z_rhs[i - rk_rhs_stride];
+                step_mom_x[i] = mom_x[i] + fact * mom_x_rhs[i - rk_rhs_stride];
+                step_mom_z[i] = mom_z[i] + fact * mom_z_rhs[i - rk_rhs_stride];
+            }
+            break;
+        }
+    }
+
+    // interpolate the metric variables at this step's coordinates
+    interp_geometry(ns, surfaces_active, step_pos_x, pos_y, step_pos_z);
+
+    // eval RHS
+    null_geodesic_rhs(ns->photons_per_surface * surfaces_active,
+                      pos_x_rhs, pos_z_rhs, mom_x_rhs, mom_z_rhs,
+                      step_pos_x, step_pos_z, step_mom_x, step_mom_z,
+                      ns->interp_values);
 }
-static double diff4_dz(const double *arr, ptrdiff_t stride, double h_inv)
+
+static void ns_evol_rk_finish(NSContext *ns, int surfaces_active)
 {
-    return FD4(arr, stride, h_inv);
+    const ptrdiff_t rk_rhs_stride = ns->nb_surfaces * ns->photons_per_surface;
+
+    double *pos_x = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::pos_x");
+    double *pos_z = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::pos_z");
+    double *mom_x = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_x");
+    double *mom_z = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_z");
+
+    const double *pos_x_rhs = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::pos_x_rhs");
+    const double *pos_z_rhs = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::pos_z_rhs");
+    const double *mom_x_rhs = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_x_rhs");
+    const double *mom_z_rhs = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_z_rhs");
+
+    const double *pos_x_rhs0 = pos_x_rhs + 0 * rk_rhs_stride;
+    const double *pos_x_rhs1 = pos_x_rhs + 1 * rk_rhs_stride;
+    const double *pos_x_rhs2 = pos_x_rhs + 2 * rk_rhs_stride;
+    const double *pos_x_rhs3 = pos_x_rhs + 3 * rk_rhs_stride;
+
+    const double *pos_z_rhs0 = pos_z_rhs + 0 * rk_rhs_stride;
+    const double *pos_z_rhs1 = pos_z_rhs + 1 * rk_rhs_stride;
+    const double *pos_z_rhs2 = pos_z_rhs + 2 * rk_rhs_stride;
+    const double *pos_z_rhs3 = pos_z_rhs + 3 * rk_rhs_stride;
+
+    const double *mom_x_rhs0 = mom_x_rhs + 0 * rk_rhs_stride;
+    const double *mom_x_rhs1 = mom_x_rhs + 1 * rk_rhs_stride;
+    const double *mom_x_rhs2 = mom_x_rhs + 2 * rk_rhs_stride;
+    const double *mom_x_rhs3 = mom_x_rhs + 3 * rk_rhs_stride;
+
+    const double *mom_z_rhs0 = mom_z_rhs + 0 * rk_rhs_stride;
+    const double *mom_z_rhs1 = mom_z_rhs + 1 * rk_rhs_stride;
+    const double *mom_z_rhs2 = mom_z_rhs + 2 * rk_rhs_stride;
+    const double *mom_z_rhs3 = mom_z_rhs + 3 * rk_rhs_stride;
+
+#pragma omp parallel for
+    for (int i = 0; i < surfaces_active * ns->photons_per_surface; i++) {
+        pos_x[i] += ns->dt * (pos_x_rhs0[i] / 6. + pos_x_rhs1[i] / 3. + pos_x_rhs2[i] / 3. + pos_x_rhs3[i] / 6.);
+        pos_z[i] += ns->dt * (pos_z_rhs0[i] / 6. + pos_z_rhs1[i] / 3. + pos_z_rhs2[i] / 3. + pos_z_rhs3[i] / 6.);
+        mom_x[i] += ns->dt * (mom_x_rhs0[i] / 6. + mom_x_rhs1[i] / 3. + mom_x_rhs2[i] / 3. + mom_x_rhs3[i] / 6.);
+        mom_z[i] += ns->dt * (mom_z_rhs0[i] / 6. + mom_z_rhs1[i] / 3. + mom_z_rhs2[i] / 3. + mom_z_rhs3[i] / 6.);
+    }
 }
 
-void ns_eval_rhs(CCTK_ARGUMENTS)
+void ns_evol(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;
-            }
+    NSContext *ns = &ns_ctx;
 
-            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);
-                }
+    const int reflevel = ctz(cctk_levfac[0]);
 
-            shift_term = 0.0;
-            for (int i = 0; i < 3; i++)
-                shift_term += g4u[0][1 + i] * dF[i];
+    /* check if we are on the right reflevel */
+    if (reflevel != solve_level)
+        return;
 
-            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];
+    fprintf(stderr, "ns evol step %d at time %g; pos[zaxis] = %g; pos[xaxis] = %g\n",
+            *evol_next_step, cctk_time, pos_z[0], pos_x[photons_per_surface - 1]);
 
-            F_rhs[idx] = (-shift_term + sqrt(SQR(shift_term) - g4u[0][0] * dF2)) / g4u[0][0];
-        }
+    /* launch a new surface if its time has come */
+    // TODO: actually implement multiple surfaces
+    if (cctk_time == 0.0 && *surfaces_active == 0) {
+        const int offset = *surfaces_active * photons_per_surface;
+        surface_launch(ns, pos_x + offset, pos_z + offset, mom_x + offset, mom_z + offset);
+
+        *surfaces_active += 1;
+    }
+
+    /* do the Runge-Kutta evolution step(s) */
+    ns_evol_rk_substep(ns, *evol_next_step, *surfaces_active);
+    *evol_next_step += 1;
+
+    if (*evol_next_step == 4) {
+        ns_evol_rk_finish(ns, *surfaces_active);
+        *evol_next_step = 0;
+    }
+    if (*evol_next_step == 0 || *evol_next_step == 2) {
+        ns_evol_rk_substep(ns, *evol_next_step, *surfaces_active);
+        *evol_next_step += 1;
     }
 }
 
-void ns_register_symmetries(CCTK_ARGUMENTS)
+static int context_init(NSContext *ns, cGH *gh, int reflevel,
+                        int nb_surfaces, int photons_per_surface)
 {
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
+    int ret;
+
+    ns->gh       = gh;
+    ns->reflevel = reflevel;
+
+    /* we do one null surface evolution step per two cactus time steps */
+    ns->dt = 2.0 * gh->cctk_delta_time / (1 << reflevel);
+
+    ns->coord_system = CCTK_CoordSystemHandle("cart3d");
+    if (ns->coord_system < 0)
+        CCTK_WARN(0, "Error getting the coordinate system");
+
+    ns->interp_operator = CCTK_InterpHandle("Lagrange polynomial interpolation (tensor product)");
+    if (ns->interp_operator < 0)
+        CCTK_WARN(0, "Error getting the interpolation operator");
+
+    ns->interp_params = Util_TableCreateFromString("order=4 want_global_mode=1");
+    if (ns->interp_params < 0)
+        CCTK_WARN(0, "Error creating interpolation parameters table");
+
+    ret = Util_TableSetIntArray(ns->interp_params, NB_INTERP_VARS,
+                                interp_operation_codes, "operation_codes");
+    if (ret < 0)
+        CCTK_WARN(0, "Error setting operation codes");
+
+    ret = Util_TableSetIntArray(ns->interp_params, NB_INTERP_VARS,
+                                interp_operation_indices, "operand_indices");
+    if (ret < 0)
+        CCTK_WARN(0, "Error setting operand indices");
+
+    for (int i = 0; i < ARRAY_ELEMS(metric_vars); i++) {
+        ns->interp_vars_indices[i] = CCTK_VarIndex(metric_vars[i]);
+        if (ns->interp_vars_indices[i] < 0)
+            CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, "Error getting the index of variable: %s\n", metric_vars[i]);
+    }
+
+    for (int i = 0; i < ARRAY_ELEMS(ns->interp_values); i++) {
+        ret = posix_memalign((void**)&ns->interp_values[i], 32,
+                             nb_surfaces * photons_per_surface * sizeof(*ns->interp_values[i]));
+        if (ret)
+            goto fail;
+        ns->interp_value_codes[i] = CCTK_VARIABLE_REAL;
+    }
+
+    for (int i = 0; i < ARRAY_ELEMS(ns->rk_scratch); i++) {
+        ret = posix_memalign((void**)&ns->rk_scratch[i], 32,
+                             nb_surfaces * photons_per_surface * sizeof(*ns->rk_scratch[i]));
+        if (ret)
+            goto fail;
+    }
 
-  int sym[3] = { 1, 1, 1 };
-  int ret;
+    ns->nb_surfaces         = nb_surfaces;
+    ns->photons_per_surface = photons_per_surface;
 
-  ret = SetCartSymVN(cctkGH, sym, "NullSurf::F");
-  if (ret != 0)
-      CCTK_WARN(0, "Error registering symmetries");
+    return 0;
+fail:
+    return -ENOMEM;
 }
 
-void ns_select_bc(CCTK_ARGUMENTS)
+void ns_evol_init(CCTK_ARGUMENTS)
 {
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    const int reflevel = ctz(cctk_levfac[0]);
+
+    if (reflevel != solve_level)
+        return;
+
+    memset(&ns_ctx, 0, sizeof(ns_ctx));
+
+    context_init(&ns_ctx, cctkGH, reflevel, nb_surfaces, photons_per_surface);
+
+    /* initialize y-values of the photon positions */
+    memset(pos_y, 0, sizeof(*pos_y) * nb_surfaces * photons_per_surface);
+}
+
+void ns_init(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    const int reflevel = ctz(cctk_levfac[0]);
+
+    if (reflevel != solve_level)
+        return;
+
+    for (int i = 0; i < nb_surfaces; i++)
+        photon_times[i] = -1;
 
-  int ret = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, "NullSurf::F", "none");
-  if (ret != 0)
-      CCTK_WARN(0, "Error registering boundaries");
+    *surfaces_active = 0;
+    *evol_next_step  = 0;
 }