Smoothly send W to zero at the boundary.

1 files changed, 47 insertions, 6 deletions

diff --git a/src/qms.c b/src/qms.c
index 5998fee..8ba8b1e 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -37,6 +37,9 @@ typedef struct CoordPatch {
 
     ptrdiff_t offset_left[2];
     int bnd_intercomp[2][2];
+
+    double *filter;
+    ptrdiff_t filter_stride;
 } CoordPatch;
 
 typedef struct QMSMGContext {
@@ -54,6 +57,9 @@ typedef struct QMSMGContext {
     int solve_level;
     int boundary_offset;
 
+    double filter_start;
+    double filter_end;
+
     CoordPatch *patches;
     int nb_patches;
 } QMSMGContext;
@@ -208,6 +214,30 @@ static CoordPatch *get_coord_patch(QMSMGContext *ms, int level)
         }
     }
 
+    if (level == ms->solve_level) {
+        ms->filter_end = solver->step[0] * (solver->domain_size - 1);
+        ms->filter_start = 0.7 * ms->filter_end;
+    }
+
+    if ((solver->domain_size - 1) * solver->step[0] >= ms->filter_start) {
+        cp->filter_stride = solver->domain_size;
+        cp->filter = calloc(solver->domain_size * cp->filter_stride, sizeof(*cp->filter));
+        if (!cp->filter)
+            CCTK_WARN(0, "Error allocating the filter");
+
+        for (int idx_z = 0; idx_z < solver->domain_size; idx_z++) {
+            const double z = solver->step[1] * idx_z;
+            for (int idx_x = 0; idx_x < solver->domain_size; idx_x++) {
+                const double  x = solver->step[0] * idx_x;
+                const double  r = sqrt(SQR(x) + SQR(z));
+
+                double fact = exp(-36.0 * (pow((r - ms->filter_start) / (ms->filter_end - ms->filter_start), 6.0)));
+                fact = r >= ms->filter_start ? fact : 1.0;
+
+                cp->filter[idx_z * cp->filter_stride + idx_x] = fact;
+            }
+        }
+    }
 
     ms->nb_patches++;
     return cp;
@@ -832,13 +862,24 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
 
 static void solution_to_grid(CoordPatch *cp, MG2DContext *solver, double *dst)
 {
+    if (cp->filter) {
 #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++)
+            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;
+                const int idx_filter = j * cp->filter_stride + i;
+                dst[idx_dst] = solver->u[idx_src] * cp->filter[idx_filter];
+            }
+    } 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 */
     if (!cp->bnd_intercomp[0][0]) {