Smooth the solution at the outer boundary.

2 files changed, 52 insertions, 20 deletions

diff --git a/param.ccl b/param.ccl
index 1e8e2b5..bbbca32 100644
--- a/param.ccl
+++ b/param.ccl
@@ -27,6 +27,11 @@ CCTK_INT ccz4 "Use CCZ4 or BSSN" STEERABLE=recover
     0:1 :: ""
 } 0
 
+CCTK_REAL outer_smooth_fact "" STEERABLE=recover
+{
+    : :: ""
+} 0.5
+
 RESTRICTED:
 CCTK_INT fd_stencil "finite differencing stencil"
 {
diff --git a/src/qms.c b/src/qms.c
index 4ed07fd..b34c590 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -82,6 +82,7 @@ typedef struct QMSMGContext {
     int nb_relax_pre;
     double tol_residual_base;
     double cfl_factor;
+    double outer_smooth_fact;
     int adaptive_step;
     int solve_level;
     int solve_level_max;
@@ -1108,15 +1109,47 @@ static void mirror(CoordPatch *cp, double *dst)
     }
 }
 
-static void solution_to_grid(CoordPatch *cp, MG2DContext *solver, double *dst)
+static double smooth_step(double x)
 {
+    if (x <= 0.0)
+        return 0.0;
+    else if (x >= 1.0)
+        return 1.0;
+    return 1.0 - exp(-36.0 * pow(x, 6.0));
+}
+
+static void solution_to_grid(QMSMGContext *ms, CoordPatch *cp, MG2DContext *solver, double *dst)
+{
+    /* on the coarsest level, smooth the solution at the outer boundary */
+    if (cp->level == ms->solve_level_max) {
+        const double bnd_loc = (solver->domain_size - 1) * solver->step[0];
+        const double transition_start = bnd_loc * ms->outer_smooth_fact;
+
 #pragma omp parallel for
-    for (int j = 0; j < solver->local_size[1]; j++)
-        for (int i = 0; i < solver->local_size[0]; i++) {
-            const ptrdiff_t idx_dst = CPINDEX(cp, i + cp->offset_left[0], cp->y_idx, j + cp->offset_left[1]);
-            const int idx_src = j * solver->u_stride + i;
-            dst[idx_dst] = solver->u[idx_src];
+        for (int j = 0; j < solver->local_size[1]; j++) {
+            const double z_val = (j + solver->local_start[1]) * solver->step[1];
+
+            for (int i = 0; i < solver->local_size[0]; i++) {
+                const double x_val = (i + solver->local_start[0]) * solver->step[0];
+                const double r = sqrt(SQR(x_val) + SQR(z_val));
+                const double dist = (r - transition_start) / (bnd_loc - transition_start);
+
+                const double smooth_fact = 1.0 - smooth_step(dist);
+
+                const ptrdiff_t idx_dst = CPINDEX(cp, i + cp->offset_left[0], cp->y_idx, j + cp->offset_left[1]);
+                const int idx_src = j * solver->u_stride + i;
+                dst[idx_dst] = solver->u[idx_src] * smooth_fact;
+            }
         }
+    } else {
+#pragma omp parallel for
+        for (int j = 0; j < solver->local_size[1]; j++)
+            for (int i = 0; i < solver->local_size[0]; i++) {
+                const ptrdiff_t idx_dst = CPINDEX(cp, i + cp->offset_left[0], cp->y_idx, j + cp->offset_left[1]);
+                const int idx_src = j * solver->u_stride + i;
+                dst[idx_dst] = solver->u[idx_src];
+            }
+    }
 
     /* fill in the axial symmetry ghostpoints by mirroring */
     mirror(cp, dst);
@@ -1149,7 +1182,8 @@ static int rect_intersect(Rect *dst, const Rect *src1, const Rect *src2)
     return 0;
 }
 
-static int guess_from_coarser(QMSMGContext *ms, CoordPatch *cp, CoordPatch *cp_coarse)
+static int guess_from_coarser(QMSMGContext *ms, CoordPatch *cp, CoordPatch *cp_coarse,
+                              double *dst)
 {
     int ret;
 
@@ -1259,15 +1293,6 @@ static int guess_from_coarser(QMSMGContext *ms, CoordPatch *cp, CoordPatch *cp_c
     return ret;
 }
 
-static double smooth_step(double x)
-{
-    if (x <= 0.0)
-        return 0.0;
-    else if (x >= 1.0)
-        return 1.0;
-    return 1.0 - exp(-36.0 * pow(x, 6.0));
-}
-
 void qms_mg_solve(CCTK_ARGUMENTS)
 {
     QMSMGContext *ms = qms_context;
@@ -1333,7 +1358,7 @@ void qms_mg_solve(CCTK_ARGUMENTS)
             /* use the solution from the coarser level as the initial guess */
             CoordPatch *cp_coarse = get_coord_patch(ms, reflevel - 1);
 
-            ret = guess_from_coarser(ms, cp, cp_coarse);
+            ret = guess_from_coarser(ms, cp, cp_coarse, W_val);
             if (ret < 0)
                 CCTK_WARN(0, "Error setting the initial guess");
         }
@@ -1378,7 +1403,7 @@ skip_solve:
     start = gettime();
     if (reflevel >= ms->solve_level_max) {
         double *W_val_tl1 = CCTK_VarDataPtr(cctkGH, 1, "QuasiMaximalSlicingMG::W_val");
-        solution_to_grid(cp, cp->solver, W_val);
+        solution_to_grid(ms, cp, cp->solver, W_val);
         memcpy(W_val_tl1, W_val, grid_size * sizeof(*W_val_tl1));
     }
 
@@ -1573,7 +1598,8 @@ void qms_mg_eval(CCTK_ARGUMENTS)
 
 static int context_init(cGH *gh, int solve_level, int solve_level_max, int fd_stencil, int maxiter, int exact_size, int nb_cycles,
                         int nb_relax_pre, int nb_relax_post, double tol_residual_base,
-                        double cfl_factor, int adaptive_step, const char *loglevel_str, int boundary_offset,
+                        double cfl_factor, double outer_smooth_fact, int adaptive_step,
+                        const char *loglevel_str, int boundary_offset,
                         QMSMGContext **ctx)
 {
     QMSMGContext *qms;
@@ -1649,7 +1675,8 @@ void qms_mg_init(CCTK_ARGUMENTS)
     if (!qms_context) {
         ret = context_init(cctkGH, solve_level, solve_level_max, fd_stencil, maxiter, exact_size, nb_cycles,
                            nb_relax_pre, nb_relax_post, tol_residual_base,
-                           cfl_factor, adaptive_step, loglevel, boundary_offset, &qms_context);
+                           cfl_factor, outer_smooth_fact, adaptive_step,
+                           loglevel, boundary_offset, &qms_context);
         if (ret < 0)
             CCTK_WARN(0, "Error initializing the solver context");
     }