Extrapolate the solution to fill coarsest-level outer ghostpoints.

1 files changed, 17 insertions, 0 deletions

diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c
index 44fe1ba..56a7c94 100644
--- a/src/maximal_slicing_axi_mg.c
+++ b/src/maximal_slicing_axi_mg.c
@@ -500,6 +500,23 @@ static void solution_to_grid(CoordPatch *cp, double *dst)
             dst[idx_dst] = 1.0 + cp->solver->u[idx_src];
         }
 
+    if (cp->level == 0) {
+        /* on the coarsest level, extrapolate the outer boundary ghost points */
+        for (int idx_z = cp->offset_left[1]; idx_z < cp->grid_size[2] - cp->offset_right[1]; idx_z++)
+            for (int idx_x = cp->grid_size[0] - cp->offset_right[0]; idx_x < cp->grid_size[0]; idx_x++) {
+                const ptrdiff_t idx_dst = CPINDEX(cp, idx_x, cp->y_idx, idx_z);
+                dst[idx_dst] = 2.0 * dst[idx_dst - 1] - dst[idx_dst - 2];
+            }
+
+        for (int idx_x = cp->offset_left[0]; idx_x < cp->grid_size[0]; idx_x++)
+            for (int idx_z = cp->grid_size[2] - cp->offset_right[1]; idx_z < cp->grid_size[2]; idx_z++) {
+                const ptrdiff_t idx_dst = CPINDEX(cp, idx_x, cp->y_idx, idx_z);
+                const ptrdiff_t idx_src0 = CPINDEX(cp, idx_x, cp->y_idx, idx_z - 1);
+                const ptrdiff_t idx_src1 = CPINDEX(cp, idx_x, cp->y_idx, idx_z - 2);
+                dst[idx_dst] = 2.0 * dst[idx_src0] - dst[idx_src1];
+            }
+    }
+
     /* fill in the axial symmetry ghostpoints by mirroring */
     for (int idx_z = cp->offset_left[1]; idx_z < cp->grid_size[2]; idx_z++)
         for (int idx_x = 0; idx_x < cp->offset_left[0]; idx_x++) {