Implement multiple surfaces.HEADmaster

4 files changed, 37 insertions, 16 deletions

diff --git a/interface.ccl b/interface.ccl
index b0f1325..18580c9 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -15,7 +15,7 @@ USES FUNCTION VarDataPtrI
 
 
 public:
-REAL photon_coord TYPE=array DIM=2 SIZE=nb_surfaces,photons_per_surface DISTRIB=constant
+REAL photon_coord TYPE=array DIM=2 SIZE=photons_per_surface,nb_surfaces DISTRIB=constant
 {
     pos_x,
     # stored only to send zeroes to interpolator
@@ -25,7 +25,7 @@ REAL photon_coord TYPE=array DIM=2 SIZE=nb_surfaces,photons_per_surface DISTRIB=
     mom_z
 }
 
-REAL photon_rhs TYPE=array DIM=3 SIZE=4,nb_surfaces,photons_per_surface DISTRIB=constant
+REAL photon_rhs TYPE=array DIM=3 SIZE=photons_per_surface,nb_surfaces,4 DISTRIB=constant
 {
     pos_x_rhs,
     pos_z_rhs,
@@ -33,6 +33,7 @@ REAL photon_rhs TYPE=array DIM=3 SIZE=4,nb_surfaces,photons_per_surface DISTRIB=
     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
+REAL origin_proper_time TYPE=scalar
+REAL origin_proper_time_rhs TYPE=array DIM=1 SIZE=4
diff --git a/param.ccl b/param.ccl
index af9189c..25286ad 100644
--- a/param.ccl
+++ b/param.ccl
@@ -10,6 +10,11 @@ CCTK_INT nb_surfaces "Number of null surfaces to trace"
     1: :: ""
 } 1
 
+CCTK_REAL launch_delta "Proper time between subsequent surface launches"
+{
+    0: :: ""
+} 0.25
+
 CCTK_INT solve_level "Refinement level to integrate photons on"
 {
     0: :: ""
diff --git a/schedule.ccl b/schedule.ccl
index e743105..9b32340 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -16,7 +16,8 @@ if (nb_surfaces > 0) {
 
     STORAGE: photon_coord
     STORAGE: photon_rhs
-    STORAGE: photon_times
     STORAGE: evol_next_step
     STORAGE: surfaces_active
+    STORAGE: origin_proper_time
+    STORAGE: origin_proper_time_rhs
 }
diff --git a/src/nullsurf.c b/src/nullsurf.c
index 11563ec..b97ee44 100644
--- a/src/nullsurf.c
+++ b/src/nullsurf.c
@@ -473,8 +473,16 @@ ns_evol_rk_substep(NSContext *ns, int substep, int surfaces_active)
     double *mom_x = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_x");
     double *mom_z = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_z");
 
+    const int origin_idx = CCTK_GFINDEX3D(ns->gh, ns->gh->cctk_nghostzones[0],
+                                          0,      ns->gh->cctk_nghostzones[2]);
+    double *tau_rhs = (double*)CCTK_VarDataPtr(ns->gh, 0, "NullSurf::origin_proper_time_rhs") + substep;
+    double  alpha   = ((double*)CCTK_VarDataPtr(ns->gh, 0, metric_vars[ALPHA]))[origin_idx];
+
     double *step_pos_x, *step_pos_z, *step_mom_x, *step_mom_z;
 
+    /* proper time at origin */
+    *tau_rhs = alpha;
+
     /* compute the coordinate/momenta to evaluate the RHS for */
     switch (substep) {
         case 0:
@@ -548,6 +556,11 @@ static void ns_evol_rk_finish(NSContext *ns, int surfaces_active)
     const double *mom_z_rhs2 = mom_z_rhs + 2 * rk_rhs_stride;
     const double *mom_z_rhs3 = mom_z_rhs + 3 * rk_rhs_stride;
 
+    double *tau_rhs = (double*)CCTK_VarDataPtr(ns->gh, 0, "NullSurf::origin_proper_time_rhs");
+    double *tau     = (double*)CCTK_VarDataPtr(ns->gh, 0, "NullSurf::origin_proper_time");
+
+    *tau += ns->dt * (tau_rhs[0] / 6. + tau_rhs[1] / 3. + tau_rhs[2] / 3. + tau_rhs[3] / 6.);
+
 #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.);
@@ -573,24 +586,28 @@ void ns_evol(CCTK_ARGUMENTS)
     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]);
 
+    /* do the last Runge-Kutta sub-step */
+    if (*evol_next_step == 3) {
+        ns_evol_rk_substep(ns, *evol_next_step, *surfaces_active);
+        ns_evol_rk_finish(ns, *surfaces_active);
+        *evol_next_step = 0;
+    }
+
     /* launch a new surface if its time has come */
-    // TODO: actually implement multiple surfaces
-    if (cctk_time == 0.0 && *surfaces_active == 0) {
+    if (*surfaces_active < nb_surfaces &&
+        *origin_proper_time >= *surfaces_active * launch_delta) {
         const int offset = *surfaces_active * photons_per_surface;
+        fprintf(stderr, "launching surface %d at t=%g/tau=%g\n",
+                *surfaces_active, cctk_time, *origin_proper_time);
         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) */
+    /* do the Runge-Kutta evolution step(s): either 0->1 or 1->2 */
     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) {
+    if (*evol_next_step == 2) {
         ns_evol_rk_substep(ns, *evol_next_step, *surfaces_active);
         *evol_next_step += 1;
     }
@@ -686,9 +703,6 @@ void ns_init(CCTK_ARGUMENTS)
     if (reflevel != solve_level)
         return;
 
-    for (int i = 0; i < nb_surfaces; i++)
-        photon_times[i] = -1;
-
     *surfaces_active = 0;
     *evol_next_step  = 0;
 }