ell_relax: implement correct handling of the boundary corners

1 files changed, 48 insertions, 12 deletions

diff --git a/ell_relax.c b/ell_relax.c
index 57c4112..891fc27 100644
--- a/ell_relax.c
+++ b/ell_relax.c
@@ -310,36 +310,72 @@ static void boundaries_apply(EllRelaxContext *ctx)
     }
 
     /* fill in the corner ghosts */
-    {
+    if (ctx->boundaries[ELL_RELAX_BOUNDARY_0L].type == ELL_RELAX_BC_TYPE_FIXDIFF ||
+        ctx->boundaries[ELL_RELAX_BOUNDARY_1L].type == ELL_RELAX_BC_TYPE_FIXDIFF) {
         double *dst = ctx->u;
+        int fact_x = -1, fact_z = -1;
+
+        if (ctx->boundaries[ELL_RELAX_BOUNDARY_0L].type == ELL_RELAX_BC_TYPE_FIXDIFF)
+            fact_z *= -1;
+        else
+            fact_x *= -1;
+
         for (int i = 1; i <= FD_STENCIL_MAX; i++)
             for (int j = 1; j <= FD_STENCIL_MAX; j++) {
-                const ptrdiff_t idx = j * strides[1] + i;
-                dst[-idx] = dst[idx];
+                const ptrdiff_t idx_dst = -j * strides[1] - i;
+                const ptrdiff_t idx_src = fact_z * j * strides[1] + fact_x * i;
+                dst[idx_dst] = dst[idx_src];
             }
     }
-    {
+    if (ctx->boundaries[ELL_RELAX_BOUNDARY_0L].type == ELL_RELAX_BC_TYPE_FIXDIFF ||
+        ctx->boundaries[ELL_RELAX_BOUNDARY_1U].type == ELL_RELAX_BC_TYPE_FIXDIFF) {
         double *dst = ctx->u + ctx->domain_size[0] - 1;
+        int fact_x = 1, fact_z = -1;
+
+        if (ctx->boundaries[ELL_RELAX_BOUNDARY_0L].type == ELL_RELAX_BC_TYPE_FIXDIFF)
+            fact_z *= -1;
+        else
+            fact_x *= -1;
+
         for (int i = 1; i <= FD_STENCIL_MAX; i++)
             for (int j = 1; j <= FD_STENCIL_MAX; j++) {
-                const ptrdiff_t idx = j * strides[1] - i;
-                dst[-idx] = dst[idx];
+                const ptrdiff_t idx_dst = -j * strides[1] + i;
+                const ptrdiff_t idx_src = fact_z * j * strides[1] + fact_x * i;
+                dst[idx_dst] = dst[idx_src];
             }
     }
-    {
+    if (ctx->boundaries[ELL_RELAX_BOUNDARY_1L].type == ELL_RELAX_BC_TYPE_FIXDIFF ||
+        ctx->boundaries[ELL_RELAX_BOUNDARY_0U].type == ELL_RELAX_BC_TYPE_FIXDIFF) {
         double *dst = ctx->u + strides[1] * (ctx->domain_size[1] - 1);
+        int fact_x = -1, fact_z = 1;
+
+        if (ctx->boundaries[ELL_RELAX_BOUNDARY_0U].type == ELL_RELAX_BC_TYPE_FIXDIFF)
+            fact_z *= -1;
+        else
+            fact_x *= -1;
+
         for (int i = 1; i <= FD_STENCIL_MAX; i++)
             for (int j = 1; j <= FD_STENCIL_MAX; j++) {
-                const ptrdiff_t idx = -j * strides[1] + i;
-                dst[-idx] = dst[idx];
+                const ptrdiff_t idx_dst = j * strides[1] - i;
+                const ptrdiff_t idx_src = fact_z * j * strides[1] + fact_x * i;
+                dst[idx_dst] = dst[idx_src];
             }
     }
-    {
+    if (ctx->boundaries[ELL_RELAX_BOUNDARY_0U].type == ELL_RELAX_BC_TYPE_FIXDIFF ||
+        ctx->boundaries[ELL_RELAX_BOUNDARY_1U].type == ELL_RELAX_BC_TYPE_FIXDIFF) {
         double *dst = ctx->u + strides[1] * (ctx->domain_size[1] - 1) + ctx->domain_size[0] - 1;
+        int fact_x = 1, fact_z = 1;
+
+        if (ctx->boundaries[ELL_RELAX_BOUNDARY_0U].type == ELL_RELAX_BC_TYPE_FIXDIFF)
+            fact_z *= -1;
+        else
+            fact_x *= -1;
+
         for (int i = 1; i <= FD_STENCIL_MAX; i++)
             for (int j = 1; j <= FD_STENCIL_MAX; j++) {
-                const ptrdiff_t idx = -j * strides[1] - i;
-                dst[-idx] = dst[idx];
+                const ptrdiff_t idx_dst = j * strides[1] + i;
+                const ptrdiff_t idx_src = fact_z * j * strides[1] + fact_x * i;
+                dst[idx_dst] = dst[idx_src];
             }
     }
     ctx->time_boundaries += gettime() - start;