Implement higher-order finite differences.

2 files changed, 53 insertions, 23 deletions

diff --git a/param.ccl b/param.ccl
index e6b39e2..1718290 100644
--- a/param.ccl
+++ b/param.ccl
@@ -8,3 +8,7 @@ EXTENDS KEYWORD lapse_evolution_method
 }
 
 RESTRICTED:
+CCTK_INT fd_stencil "finite differencing stencil"
+{
+    1: :: ""
+} 1
diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c
index b5255cc..ae851db 100644
--- a/src/maximal_slicing_axi_mg.c
+++ b/src/maximal_slicing_axi_mg.c
@@ -67,6 +67,8 @@ typedef struct CoordPatch {
 typedef struct MSMGContext {
     cGH *gh;
 
+    int fd_stencil;
+
     CoordPatch *patches;
     int nb_patches;
 
@@ -203,7 +205,7 @@ static CoordPatch *get_coord_patch(MSMGContext *ms, int level)
     cp->solver->step[0] = a_x[1] - a_x[0];
     cp->solver->step[1] = cp->solver->step[0];
 
-    cp->solver->fd_stencil = 1;
+    cp->solver->fd_stencil = ms->fd_stencil;
 
     cp->solver->boundaries[MG2D_BOUNDARY_0L]->type = MG2D_BC_TYPE_FIXDIFF;
     cp->solver->boundaries[MG2D_BOUNDARY_1L]->type = MG2D_BC_TYPE_FIXDIFF;
@@ -216,8 +218,11 @@ static CoordPatch *get_coord_patch(MSMGContext *ms, int level)
 
     /* initialize boundary values to zero,
      * non-zero values on the outer boundaries of refined levels are filled in elsewhere */
-    for (int i = 0; i < 4; i++)
-        memset(cp->solver->boundaries[i]->val, 0, sizeof(double) * cp->solver->boundaries[i]->val_len);
+    for (int i = 0; i < 4; i++) {
+        for (size_t j = 0; j < cp->solver->fd_stencil; j++)
+            memset(cp->solver->boundaries[i]->val + j * cp->solver->boundaries[i]->val_stride - j,
+                   0, sizeof(*cp->solver->boundaries[i]->val) * (domain_size[0] + 2 * j));
+    }
 
     if (cp->level > 0) {
         cp->boundary_val_stride[0] = gh->cctk_lsh[0] - cp->offset_left[0];
@@ -276,7 +281,7 @@ static void print_stats(MSMGContext *ms)
     }
 }
 
-static int context_init(cGH *gh, MSMGContext **ctx)
+static int context_init(cGH *gh, int fd_stencil, MSMGContext **ctx)
 {
     MSMGContext *ms;
     int ret;
@@ -286,6 +291,7 @@ static int context_init(cGH *gh, MSMGContext **ctx)
         return -ENOMEM;
 
     ms->gh = gh;
+    ms->fd_stencil = fd_stencil;
 
     *ctx = ms;
 
@@ -537,13 +543,21 @@ void msa_mg_eval(CCTK_ARGUMENTS)
                 lapse_mg_eval[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx];
             }
 
-        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++) {
-            const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1]);
-            cp->solver->boundaries[MG2D_BOUNDARY_0U]->val[i] = lapse_mg_eval[idx] - 1.0;
+        for (int j = 0; j < cp->solver->fd_stencil; j++) {
+            MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_0U];
+            double *dst = bnd->val + j * bnd->val_stride;
+            for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) {
+                const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1] + j);
+                dst[i] = lapse_mg_eval[idx] - 1.0;
+            }
         }
-        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++) {
-            const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0], cp->y_idx, cp->offset_left[1] + i);
-            cp->solver->boundaries[MG2D_BOUNDARY_1U]->val[i] = lapse_mg_eval[idx] - 1.0;
+        for (int j = 0; j < cp->solver->fd_stencil; j++) {
+            MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_1U];
+            double *dst = bnd->val + j * bnd->val_stride;
+            for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) {
+                const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0] + j, cp->y_idx, cp->offset_left[1] + i);
+                dst[i] = lapse_mg_eval[idx] - 1.0;
+            }
         }
         ms->time_fine_boundaries += gettime() - start;
 
@@ -642,26 +656,38 @@ void msa_mg_solve(CCTK_ARGUMENTS)
 
         LOGDEBUG( "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0);
 
-        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++) {
-            const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1]);
-            lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx];
+        for (int j = 0; j < cp->solver->fd_stencil; j++) {
+            for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) {
+                const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1] + j);
+                lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx];
+            }
         }
-        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++) {
-            const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0], cp->y_idx, cp->offset_left[1] + i);
-            lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx];
+        for (int j = 0; j < cp->solver->fd_stencil; j++) {
+            for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) {
+                const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0] + j, cp->y_idx, cp->offset_left[1] + i);
+                lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx];
+            }
         }
     }
     /* if the condition above was false, then lapse_mg should be filled by
      * prolongation from the coarser level */
 
     /* fill the solver boundary conditions */
-    for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++) {
-        const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1]);
-        cp->solver->boundaries[MG2D_BOUNDARY_0U]->val[i] = lapse_mg[idx] - 1.0;
+    for (int j = 0; j < cp->solver->fd_stencil; j++) {
+        MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_0U];
+        double *dst = bnd->val + j * bnd->val_stride;
+        for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) {
+            const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1] + j);
+            dst[i] = lapse_mg[idx] - 1.0;
+        }
     }
-    for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++) {
-        const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0], cp->y_idx, cp->offset_left[1] + i);
-        cp->solver->boundaries[MG2D_BOUNDARY_1U]->val[i] = lapse_mg[idx] - 1.0;
+    for (int j = 0; j < cp->solver->fd_stencil; j++) {
+        MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_1U];
+        double *dst = bnd->val + j * bnd->val_stride;
+        for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) {
+            const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0] + j, cp->y_idx, cp->offset_left[1] + i);
+            dst[i] = lapse_mg[idx] - 1.0;
+        }
     }
     ms->time_solve_boundaries += gettime() - start;
 
@@ -790,7 +816,7 @@ void msa_mg_init(CCTK_ARGUMENTS)
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
 
-    context_init(cctkGH, &ms);
+    context_init(cctkGH, fd_stencil, &ms);
 }
 
 void msa_mg_inithist(CCTK_ARGUMENTS)