Implement C/r falloff boundary condition.

API bump.

1 files changed, 70 insertions, 0 deletions

diff --git a/ell_grid_solve.c b/ell_grid_solve.c
index 1952f95..7a883e0 100644
--- a/ell_grid_solve.c
+++ b/ell_grid_solve.c
@@ -140,6 +140,38 @@ static void boundaries_apply_reflect(double *dst, const ptrdiff_t dst_stride[2],
     }
 }
 
+/* FIXME we currently use 1st-order forward differences in all cases,
+ * since centered and/or higher-order FDs seem to be unstable
+ * this may be suboptimal and should be investigated closer */
+static const double falloff_bnd_fd_factors[][FD_STENCIL_MAX * 2 + 1] = {
+    {  1.0 / 2.0, -1.0 / 2.0, },
+    {  1.0 / 2.0, -1.0 / 2.0, },
+    //{ -1.0, 8.0, 0.0, -8.0, 1.0 },
+};
+
+static void boundaries_apply_falloff(double *dst, const ptrdiff_t dst_stride[2],
+                                     const double *src, ptrdiff_t src_stride,
+                                     size_t boundary_size, EGSContext *ctx)
+{
+    EGSInternal *priv = ctx->priv;
+
+    for (ptrdiff_t bnd_layer = 1; bnd_layer <= ctx->fd_stencil; bnd_layer++)
+        for (ptrdiff_t bnd_idx = -(ptrdiff_t)ctx->fd_stencil; bnd_idx < (ptrdiff_t)(boundary_size + ctx->fd_stencil); bnd_idx++) {
+            const double x = ctx->step[0] * (ctx->domain_size[0] - 1 + bnd_layer - ctx->fd_stencil);
+            const double y = ctx->step[0] * bnd_idx;
+            const double r = sqrt(SQR(x) + SQR(y));
+
+            double *row = dst + bnd_layer * dst_stride[1] + bnd_idx * dst_stride[0];
+
+            double val = row[-(ptrdiff_t)ctx->fd_stencil * dst_stride[1]] / r * (x / r);
+
+            for (int i = 1; i < ctx->fd_stencil * 2 + 1; i++)
+                val += dst[(bnd_layer - i) * dst_stride[1] + bnd_idx * dst_stride[0]] * falloff_bnd_fd_factors[ctx->fd_stencil - 1][i] / ctx->step[0];
+
+            *row = -val * (ctx->step[0] / falloff_bnd_fd_factors[ctx->fd_stencil - 1][0]);
+        }
+}
+
 static void boundaries_apply(EGSContext *ctx)
 {
     EGSInternal *priv = ctx->priv;
@@ -147,6 +179,7 @@ static void boundaries_apply(EGSContext *ctx)
     int64_t start;
 
     start = gettime();
+    for (int tmp = 0; tmp < 2; tmp++)
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
         const int ci = mg2d_bnd_coord_idx(i);
         const ptrdiff_t stride_boundary = strides[!ci];
@@ -166,6 +199,10 @@ static void boundaries_apply(EGSContext *ctx)
             boundaries_apply_reflect(dst, dst_strides, ctx->boundaries[i]->val,
                                      ctx->boundaries[i]->val_stride, size_boundary);
             break;
+        case MG2D_BC_TYPE_FALLOFF:
+            boundaries_apply_falloff(dst, dst_strides, ctx->boundaries[i]->val,
+                                     ctx->boundaries[i]->val_stride, size_boundary, ctx);
+            break;
         }
     }
 
@@ -426,6 +463,39 @@ static int solve_exact(EGSContext *ctx)
                             }
                     }
 
+                    // falloff
+                    for (int bnd_loc = 0; bnd_loc < ARRAY_ELEMS(ctx->boundaries); bnd_loc++) {
+                        MG2DBoundary *bnd = ctx->boundaries[bnd_loc];
+
+                        const int ci  = mg2d_bnd_coord_idx(bnd_loc);
+                        const int dir = mg2d_bnd_out_dir(bnd_loc);
+
+                        const ptrdiff_t dst_stride[2] = { 1, row_stride };
+                        const ptrdiff_t bnd_stride[2] = { dst_stride[!ci], dst_stride[ci] };
+                        const ptrdiff_t bnd_len[2]    = { ctx->domain_size[!ci], ctx->domain_size[ci] };
+
+                        double *dst = mat_row + mg2d_bnd_is_upper(bnd_loc) * ((bnd_len[1] - 1) * bnd_stride[1]);
+
+                        if (!is_bnd[bnd_loc] || bnd->type != MG2D_BC_TYPE_FALLOFF)
+                            continue;
+
+                        for (ptrdiff_t bnd_layer = ctx->fd_stencil; bnd_layer > 0; bnd_layer--)
+                            for (ptrdiff_t bnd_idx = -(ptrdiff_t)ctx->fd_stencil; bnd_idx < (ptrdiff_t)(bnd_len[0] + ctx->fd_stencil); bnd_idx++) {
+                                const double x = ctx->step[0] * (ctx->domain_size[0] - 1 + bnd_layer - ctx->fd_stencil);
+                                const double y = ctx->step[0] * bnd_idx;
+                                const double r = sqrt(SQR(x) + SQR(y));
+                                double val = dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]];
+
+                                for (int i = 1; i < ctx->fd_stencil * 2 + 1; i++) {
+                                    dst[dir * (bnd_layer - i) * bnd_stride[1] + bnd_idx * bnd_stride[0]] -=
+                                        val * falloff_bnd_fd_factors[ctx->fd_stencil - 1][i] / falloff_bnd_fd_factors[ctx->fd_stencil - 1][0];
+                                }
+                                dst[dir * (bnd_layer - (ptrdiff_t)ctx->fd_stencil) * bnd_stride[1] + bnd_idx * bnd_stride[0]] -= val * (x / SQR(r)) * (ctx->step[0] / falloff_bnd_fd_factors[ctx->fd_stencil - 1][0]);
+
+                                dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]] = 0.0;
+                            }
+                    }
+
                     /* copy the interior values */
                     for (ptrdiff_t idx1_col = 0; idx1_col < ctx->domain_size[1]; idx1_col++)
                         for (ptrdiff_t idx0_col = 0; idx0_col < ctx->domain_size[0]; idx0_col++) {