Use extrapolated solution as initial guess for outermost solve level.

1 files changed, 10 insertions, 16 deletions

diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c
index d9d1508..23051df 100644
--- a/src/maximal_slicing_axi_mg.c
+++ b/src/maximal_slicing_axi_mg.c
@@ -1339,7 +1339,7 @@ void msa_mg_solve(CCTK_ARGUMENTS)
     if (reflevel > 0) {
         if (timestep % 2) {
             /* outer-most level for this solve, use extrapolated lapse as the
-             * boundary condition */
+             * boundary condition and initial guess */
             double time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel];
             double fact0, fact1;
 
@@ -1357,24 +1357,18 @@ 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);
 
-            if (cp->solver->local_start[1] + cp->solver->local_size[1] == cp->solver->domain_size) {
 #pragma omp parallel for
-                for (int j = 0; j < cp->solver->fd_stencil; j++) {
-                    for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->local_size[0] + j; i++) {
-                        const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cp->offset_left[1] + cp->solver->local_size[1] - 1 + j);
-                        lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx];
-                    }
-                }
-            }
-            if (cp->solver->local_start[0] + cp->solver->local_size[0] == cp->solver->domain_size) {
+        for (size_t i = 0; i < grid_size; i++)
+            lapse_mg[i] = fact0 * lapse_prev0[i] + fact1 * lapse_prev1[i];
+
 #pragma omp parallel for
-                for (int j = 0; j < cp->solver->fd_stencil; j++) {
-                    for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->local_size[1] + j; i++) {
-                        const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cp->offset_left[0] + cp->solver->local_size[0] - 1 + j, 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->local_size[1]; j++)
+            for (int i = 0; i < cp->solver->local_size[0]; i++) {
+                const ptrdiff_t idx_src = CPINDEX(cp, i + cp->offset_left[0], cp->y_idx, j + cp->offset_left[1]);
+                const int idx_dst = j * cp->solver->u_stride + i;
+                cp->solver->u[idx_dst] = lapse_mg[idx_src] - 1.0;
             }
+
         } else {
             /* use the solution from the coarser level as the initial guess */
             CoordPatch *cp1 = get_coord_patch(ms, reflevel - 1);