egs: only apply the fixval condition when it is non-zero and on init

Both u and residual are initialized to zero on alloc and the residual is
not computed on the fixval boundaries. The value of u will thus never
change there, so we do not need to impose the condition there except at
the very beginning.

1 files changed, 25 insertions, 9 deletions

diff --git a/ell_grid_solve.c b/ell_grid_solve.c
index 63aed59..4ad671a 100644
--- a/ell_grid_solve.c
+++ b/ell_grid_solve.c
@@ -73,6 +73,8 @@ struct EGSInternal {
     size_t residual_calc_size[2];
     double fd_factors[MG2D_DIFF_COEFF_NB];
 
+    int bnd_zero[4];
+
     ResidualCalcContext *rescalc;
 
     union {
@@ -176,7 +178,7 @@ static void boundaries_apply_falloff(double *dst, const ptrdiff_t dst_stride[2],
         }
 }
 
-static void boundaries_apply(EGSContext *ctx)
+static void boundaries_apply(EGSContext *ctx, int init)
 {
     static const enum MG2DBCType bnd_type_order[] = { MG2D_BC_TYPE_FIXVAL, MG2D_BC_TYPE_REFLECT, MG2D_BC_TYPE_FALLOFF };
     EGSInternal *priv = ctx->priv;
@@ -202,8 +204,10 @@ static void boundaries_apply(EGSContext *ctx)
 
         switch (bnd->type) {
         case MG2D_BC_TYPE_FIXVAL:
-            boundaries_apply_fixval(dst, dst_strides, bnd->val,
-                                    bnd->val_stride, size_boundary);
+            if (init && !priv->bnd_zero[i]) {
+                boundaries_apply_fixval(dst, dst_strides, bnd->val,
+                                        bnd->val_stride, size_boundary);
+            }
             break;
         case MG2D_BC_TYPE_REFLECT:
             boundaries_apply_reflect(dst, dst_strides, bnd->val,
@@ -278,7 +282,7 @@ static int solve_relax_step(EGSContext *ctx)
     r->time_correct += gettime() - start;
     r->count_correct++;
 
-    boundaries_apply(ctx);
+    boundaries_apply(ctx, 0);
     residual_calc(ctx);
 
     return 0;
@@ -566,7 +570,7 @@ static int solve_exact(EGSContext *ctx)
     ec->time_lin_solve += gettime() - start;
     ec->count_lin_solve++;
 
-    boundaries_apply(ctx);
+    boundaries_apply(ctx, 0);
     residual_calc(ctx);
 
     return 0;
@@ -648,10 +652,22 @@ int mg2di_egs_init(EGSContext *ctx)
     if (ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXVAL)
         priv->residual_calc_offset++;
 
-    for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
-        MG2DBoundary *bnd = ctx->boundaries[i];
-        const int ci  = mg2d_bnd_coord_idx(i);
+    for (int bnd_idx = 0; bnd_idx < ARRAY_ELEMS(ctx->boundaries); bnd_idx++) {
+        MG2DBoundary *bnd = ctx->boundaries[bnd_idx];
+        const int ci  = mg2d_bnd_coord_idx(bnd_idx);
+
         if (bnd->type == MG2D_BC_TYPE_FIXVAL) {
+            double maxval = 0.0;
+
+            for (int j = 0; j < ctx->fd_stencil; j++) {
+                for (ptrdiff_t i = -j; i < (ptrdiff_t)ctx->domain_size[!ci] + j; i++) {
+                    double val = fabs(bnd->val[j * bnd->val_stride + i]);
+                    maxval = MAX(maxval, val);
+                }
+            }
+
+            priv->bnd_zero[bnd_idx] = maxval == 0.0;
+
             priv->residual_calc_size[ci]--;
         }
     }
@@ -664,7 +680,7 @@ int mg2di_egs_init(EGSContext *ctx)
     if (ret < 0)
         return ret;
 
-    boundaries_apply(ctx);
+    boundaries_apply(ctx, 1);
     residual_calc(ctx);
 
     ctx->time_total += gettime() - start;