Stop filling the boundary values for coarsest level.

1 files changed, 34 insertions, 32 deletions

diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c
index 56a7c94..c83c2fe 100644
--- a/src/maximal_slicing_axi_mg.c
+++ b/src/maximal_slicing_axi_mg.c
@@ -633,7 +633,7 @@ void msa_mg_eval(CCTK_ARGUMENTS)
             }
 
         for (int j = 0; j < cp->solver->fd_stencil; j++) {
-            MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_0U];
+            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, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1] + j);
@@ -641,7 +641,7 @@ void msa_mg_eval(CCTK_ARGUMENTS)
             }
         }
         for (int j = 0; j < cp->solver->fd_stencil; j++) {
-            MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_1U];
+            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, cctkGH->cctk_lsh[0] - cp->offset_right[0] + j, cp->y_idx, cp->offset_left[1] + i);
@@ -736,48 +736,50 @@ void msa_mg_solve(CCTK_ARGUMENTS)
     ms->time_solve_fill_coeffs += gettime() - start;
 
     start = gettime();
-    if (reflevel == 0 || timestep % 2) {
-        /* outer-most level for this solve, use extrapolated lapse as the
-         * boundary condition */
-        double time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel];
-        double fact0 = (lapse_prev1_time[reflevel] - t) / time_interp_step;
-        double fact1 = (t - lapse_prev0_time[reflevel]) / time_interp_step;
-
-        LOGDEBUG( "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0);
+    if (reflevel > 0) {
+        if (timestep % 2) {
+            /* outer-most level for this solve, use extrapolated lapse as the
+             * boundary condition */
+            double time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel];
+            double fact0 = (lapse_prev1_time[reflevel] - t) / time_interp_step;
+            double fact1 = (t - lapse_prev0_time[reflevel]) / time_interp_step;
+
+            LOGDEBUG( "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0);
+
+            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 (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 (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, 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];
+                dst[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, 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];
+                dst[i] = lapse_mg[idx] - 1.0;
             }
         }
     }
-    /* if the condition above was false, then lapse_mg should be filled by
-     * prolongation from the coarser level */
-
-    /* fill the solver boundary conditions */
-    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 (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;
 
     /* do the elliptic solve */