Smoothly zero the lapse source at the outer boundary.

1 files changed, 21 insertions, 2 deletions

diff --git a/src/qms.c b/src/qms.c
index 1005290..d782d0c 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -362,8 +362,8 @@ static CoordPatch *get_coord_patch(QMSMGContext *ms, int level)
 
     solver->boundaries[MG2D_BOUNDARY_0L]->type = MG2D_BC_TYPE_REFLECT;
     solver->boundaries[MG2D_BOUNDARY_1L]->type = MG2D_BC_TYPE_REFLECT;
-    solver->boundaries[MG2D_BOUNDARY_0U]->type = level ? MG2D_BC_TYPE_FIXVAL : MG2D_BC_TYPE_FALLOFF;
-    solver->boundaries[MG2D_BOUNDARY_1U]->type = level ? MG2D_BC_TYPE_FIXVAL : MG2D_BC_TYPE_FALLOFF;
+    solver->boundaries[MG2D_BOUNDARY_0U]->type = MG2D_BC_TYPE_FIXVAL;
+    solver->boundaries[MG2D_BOUNDARY_1U]->type = MG2D_BC_TYPE_FIXVAL;
 
     solver->maxiter        = ms->maxiter;
     solver->tol            = ms->tol_residual_base / (solver->step[0] * solver->step[1]);
@@ -1420,6 +1420,25 @@ void qms_mg_solve(CCTK_ARGUMENTS)
                     W_pred1[idx] = (W_val1[idx] * fact1 + W_val0[idx] * fact0) * (1.0 - bnd_fact) + coarse_val * bnd_fact;
                 }
             mirror(cp, W_pred1);
+        } else if (reflevel == 0) {
+            const int interior_size = cctk_lsh[0] - cp->offset_left[0] - cctk_nghostzones[0];
+
+            const double bnd_loc = x[CCTK_GFINDEX3D(ms->gh, cp->offset_left[0] + interior_size, 0, 0)];
+            const double transition_start = bnd_loc * 0.7;
+
+#pragma omp parallel for
+            for (int idx = 0; idx < grid_size; idx++) {
+                const double      x_val = x[idx];
+                const double      z_val = z[idx];
+                const double coarse_val = W_pred1[idx];
+
+                const double r = sqrt(SQR(x_val) + SQR(z_val));
+                const double dist = (r - transition_start) / (bnd_loc - transition_start);
+
+                const double bnd_fact = 1.0 - smooth_step(dist);
+
+                W_pred1[idx] = (W_val1[idx] * fact1 + W_val0[idx] * fact0) * bnd_fact;
+            }
         } else {
 #pragma omp parallel for
             for (int i = 0; i < grid_size; i++)