ell_grid_solve: use the new boundary API to simplify code

1 files changed, 105 insertions, 196 deletions

diff --git a/ell_grid_solve.c b/ell_grid_solve.c
index 75960e7..0cefdfd 100644
--- a/ell_grid_solve.c
+++ b/ell_grid_solve.c
@@ -77,28 +77,6 @@ struct EGSInternal {
     TPContext *tp_internal;
 };
 
-static const struct {
-    unsigned int stride_idx;
-    unsigned int is_upper;
-} boundary_def[] = {
-    [MG2D_BOUNDARY_0L] = {
-        .stride_idx = 1,
-        .is_upper   = 0,
-    },
-    [MG2D_BOUNDARY_0U] = {
-        .stride_idx = 1,
-        .is_upper   = 1,
-    },
-    [MG2D_BOUNDARY_1L] = {
-        .stride_idx = 0,
-        .is_upper   = 0
-    },
-    [MG2D_BOUNDARY_1U] = {
-        .stride_idx = 0,
-        .is_upper   = 1,
-    },
-};
-
 static const double fd_denoms[][MG2D_DIFF_COEFF_NB] = {
     {
         [MG2D_DIFF_COEFF_00] = 1.0,
@@ -170,15 +148,14 @@ static void boundaries_apply(EGSContext *ctx)
 
     start = gettime();
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
-        const int si = boundary_def[i].stride_idx;
-        const ptrdiff_t stride_boundary = strides[si];
-        const ptrdiff_t stride_offset   = strides[!si];
-        const size_t size_boundary      = ctx->domain_size[si];
-        const size_t size_offset        = ctx->domain_size[!si];
+        const int ci = mg2d_bnd_coord_idx(i);
+        const ptrdiff_t stride_boundary = strides[!ci];
+        const ptrdiff_t stride_offset   = strides[ci];
+        const size_t size_boundary      = ctx->domain_size[!ci];
+        const size_t size_offset        = ctx->domain_size[ci];
 
-        double *dst = ctx->u + boundary_def[i].is_upper * ((size_offset - 1) * stride_offset);
-        const ptrdiff_t dst_strides[] = { stride_boundary,
-            (boundary_def[i].is_upper ? 1 : -1) * stride_offset };
+        double *dst = ctx->u + mg2d_bnd_is_upper(i) * ((size_offset - 1) * stride_offset);
+        const ptrdiff_t dst_strides[] = { stride_boundary, mg2d_bnd_out_dir(i) * stride_offset };
 
         switch (ctx->boundaries[i]->type) {
         case MG2D_BC_TYPE_FIXVAL:
@@ -193,73 +170,35 @@ static void boundaries_apply(EGSContext *ctx)
     }
 
     /* fill in the corner ghosts */
-    if (ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXDIFF ||
-        ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXDIFF) {
-        double *dst = ctx->u;
-        int fact_x = -1, fact_z = -1;
-
-        if (ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_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_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[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXDIFF ||
-        ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXDIFF) {
-        double *dst = ctx->u + ctx->domain_size[0] - 1;
-        int fact_x = 1, fact_z = -1;
-
-        if (ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_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_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[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXDIFF ||
-        ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_BC_TYPE_FIXDIFF) {
-        double *dst = ctx->u + strides[1] * (ctx->domain_size[1] - 1);
-        int fact_x = -1, fact_z = 1;
-
-        if (ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_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_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[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXDIFF ||
-        ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_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[MG2D_BOUNDARY_1U]->type == MG2D_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_dst = j * strides[1] + i;
-                const ptrdiff_t idx_src = fact_z * j * strides[1] + fact_x * i;
-                dst[idx_dst] = dst[idx_src];
-            }
+    for (int pos_y = 0; pos_y < 2; pos_y++) {
+        enum MG2DBoundaryLoc loc_y = mg2d_bnd_id(1, pos_y);
+        MG2DBoundary        *bnd_y = ctx->boundaries[loc_y];
+        const int            dir_y = mg2d_bnd_out_dir(loc_y);
+
+        for (int pos_x = 0; pos_x < 2; pos_x++) {
+            enum MG2DBoundaryLoc loc_x = mg2d_bnd_id(0, pos_x);
+            MG2DBoundary        *bnd_x = ctx->boundaries[loc_x];
+            const int            dir_x = mg2d_bnd_out_dir(loc_x);
+            int fact_x = dir_x, fact_y = dir_y;
+
+            double *dst = ctx->u
+                          + mg2d_bnd_is_upper(loc_y) * ((ctx->domain_size[1] - 1) * priv->stride)
+                          + mg2d_bnd_is_upper(loc_x) * (ctx->domain_size[0] - 1);
+
+            if (bnd_y->type == MG2D_BC_TYPE_FIXDIFF)
+                fact_y *= -1;
+            else if (bnd_x->type == MG2D_BC_TYPE_FIXDIFF)
+                fact_x *= -1;
+            else
+                continue;
+
+            for (int j = 1; j <= FD_STENCIL_MAX; j++)
+                for (int i = 1; i <= FD_STENCIL_MAX; i++) {
+                    const ptrdiff_t idx_dst = dir_y  * j * strides[1] + dir_x  * i;
+                    const ptrdiff_t idx_src = fact_y * j * strides[1] + fact_x * i;
+                    dst[idx_dst] = dst[idx_src];
+                }
+        }
     }
     ctx->time_boundaries += gettime() - start;
     ctx->count_boundaries++;
@@ -297,19 +236,6 @@ static int solve_relax_step(EGSContext *ctx)
     return 0;
 }
 
-static int bnd_id(int x_lower, int x_upper, int y_lower, int y_upper)
-{
-    if (x_lower)
-        return MG2D_BOUNDARY_0L;
-    if (x_upper)
-        return MG2D_BOUNDARY_0U;
-    if (y_lower)
-        return MG2D_BOUNDARY_1L;
-    if (y_upper)
-        return MG2D_BOUNDARY_1U;
-    return -1;
-}
-
 static void fill_mat_s1(double *mat_row, double **diff_coeffs, double *fd_factors,
                         ptrdiff_t idx_src, ptrdiff_t row_stride, ptrdiff_t row_offset)
 {
@@ -397,15 +323,12 @@ static int solve_exact(EGSContext *ctx)
     memset(e->mat, 0, SQR(e->N) * sizeof(*e->mat));
 
     for (ptrdiff_t idx1 = 0; idx1 < ctx->domain_size[1]; idx1++) {
-        const int boundary_y_lower = idx1 < ctx->fd_stencil;
-        const int boundary_y_upper = idx1 >= ctx->domain_size[1] - ctx->fd_stencil;
-
-        for (ptrdiff_t idx0 = 0; idx0 < ctx->domain_size[0]; idx0++) {
-            const int boundary_x_lower = idx0 < ctx->fd_stencil;
-            const int boundary_x_upper = idx0 >= ctx->domain_size[0] - ctx->fd_stencil;
+        int is_bnd[4], boundary;
 
-            const int boundary   = boundary_y_lower || boundary_y_upper || boundary_x_lower || boundary_x_upper;
+        is_bnd[MG2D_BOUNDARY_1L] = idx1 < ctx->fd_stencil;
+        is_bnd[MG2D_BOUNDARY_1U] = idx1 >= ctx->domain_size[1] - ctx->fd_stencil;
 
+        for (ptrdiff_t idx0 = 0; idx0 < ctx->domain_size[0]; idx0++) {
             const ptrdiff_t idx_src     = idx1 * priv->stride        + idx0;
             const ptrdiff_t mat_row_idx = idx1 * ctx->domain_size[0] + idx0;
 
@@ -414,6 +337,11 @@ static int solve_exact(EGSContext *ctx)
 
             double *mat_row;
 
+            is_bnd[MG2D_BOUNDARY_0L] = idx0 < ctx->fd_stencil;
+            is_bnd[MG2D_BOUNDARY_0U] = idx0 >= ctx->domain_size[0] - ctx->fd_stencil;
+
+            boundary = is_bnd[0] || is_bnd[1] || is_bnd[2] || is_bnd[3];
+
             if (boundary) {
                 memset(e->scratch_line, 0, e->N_ghosts * sizeof(*e->scratch_line));
                 row_stride = ctx->domain_size[0] + 2 * ctx->fd_stencil;
@@ -430,92 +358,71 @@ static int solve_exact(EGSContext *ctx)
             e->rhs[mat_row_idx] = ctx->rhs[idx_src];
 
             if (boundary) {
-                const int bnd_fixval = (idx0 == 0                       && ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXVAL) ||
-                                       (idx1 == 0                       && ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXVAL) ||
-                                       (idx0 == ctx->domain_size[1] - 1 && ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXVAL) ||
-                                       (idx1 == ctx->domain_size[0] - 1 && ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_BC_TYPE_FIXVAL);
-
                 double *mat_dst = e->mat + e->N * mat_row_idx;
-                ptrdiff_t bnd_idx  = (boundary_x_lower || boundary_x_upper) ? idx0 : idx1;
 
-                if (bnd_fixval) {
-                    const MG2DBoundary *bnd = ctx->boundaries[bnd_id(boundary_x_lower, boundary_x_upper,
-                                                                     boundary_y_lower, boundary_y_upper)];
+                const MG2DBoundary *bnd_fixval = NULL;
 
-                    mat_dst[mat_row_idx] = 1.0;
-                    e->rhs[mat_row_idx]     = bnd->val[bnd_idx];
+                for (int bnd_loc = 0; bnd_loc < ARRAY_ELEMS(ctx->boundaries); bnd_loc++) {
+                    const int ci  = mg2d_bnd_coord_idx(bnd_loc);
+                    const ptrdiff_t idx[] = { idx0, idx1 };
+
+                    if (!(idx[ci] == (ctx->domain_size[ci] - 1) * mg2d_bnd_is_upper(bnd_loc) &&
+                          ctx->boundaries[bnd_loc]->type == MG2D_BC_TYPE_FIXVAL))
+                        continue;
+
+                    bnd_fixval = ctx->boundaries[bnd_loc];
+
+                    mat_dst[mat_row_idx]    = 1.0;
+                    e->rhs[mat_row_idx]     = bnd_fixval->val[idx[!ci]];
                     memset(e->scratch_line, 0, e->N_ghosts * sizeof(*e->scratch_line));
-                } else {
+                    break;
+                }
+
+                if (!bnd_fixval) {
                     /* apply the boundary conditions to eliminate the ghostpoint values */
                     // fixdiff
-                    if (boundary_y_lower && ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXDIFF) {
-                        for (ptrdiff_t idx1_bnd = 0; idx1_bnd < ctx->fd_stencil; idx1_bnd++)
-                            for (ptrdiff_t idx0_col = -(ptrdiff_t)ctx->fd_stencil; idx0_col < (ptrdiff_t)(ctx->domain_size[0] + ctx->fd_stencil); idx0_col++) {
-                                const ptrdiff_t idx_ghost = -(idx1_bnd + 1) * row_stride + idx0_col;
-                                mat_row[(idx1_bnd + 1) * row_stride + idx0_col] += mat_row[idx_ghost];
-                                mat_row[idx_ghost] = 0.0;
-                            }
-                    }
-                    if (boundary_y_upper && ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_BC_TYPE_FIXDIFF) {
-                        for (ptrdiff_t idx1_bnd = 0; idx1_bnd < ctx->fd_stencil; idx1_bnd++)
-                            for (ptrdiff_t idx0_col = -(ptrdiff_t)ctx->fd_stencil; idx0_col < (ptrdiff_t)(ctx->domain_size[0] + ctx->fd_stencil); idx0_col++) {
-                                const ptrdiff_t idx_ghost = (ctx->domain_size[1] - 1 + idx1_bnd + 1) * row_stride + idx0_col;
-                                mat_row[(ctx->domain_size[1] - 1 - (idx1_bnd + 1)) * row_stride + idx0_col] += mat_row[idx_ghost];
-                                mat_row[idx_ghost] = 0.0;
-                            }
-                    }
-                    if (boundary_x_lower && ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXDIFF) {
-                        for (ptrdiff_t idx1_col = -(ptrdiff_t)ctx->fd_stencil; idx1_col < (ptrdiff_t)(ctx->domain_size[1] + ctx->fd_stencil); idx1_col++)
-                            for (ptrdiff_t idx0_bnd = 0; idx0_bnd < ctx->fd_stencil; idx0_bnd++) {
-                                const ptrdiff_t idx_ghost = idx1_col * row_stride - (idx0_bnd + 1);
-                                mat_row[idx1_col * row_stride + idx0_bnd + 1] += mat_row[idx_ghost];
-                                mat_row[idx_ghost] = 0.0;
-                            }
-                    }
-                    if (boundary_x_upper && ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXDIFF) {
-                        for (ptrdiff_t idx1_col = -(ptrdiff_t)ctx->fd_stencil; idx1_col < (ptrdiff_t)(ctx->domain_size[1] + ctx->fd_stencil); idx1_col++)
-                            for (ptrdiff_t idx0_bnd = 0; idx0_bnd < ctx->fd_stencil; idx0_bnd++) {
-                                const ptrdiff_t idx_ghost = idx1_col * row_stride + ctx->domain_size[0] - 1 + idx0_bnd + 1;
-                                mat_row[idx1_col * row_stride + ctx->domain_size[0] - 1 - (idx0_bnd + 1)] += mat_row[idx_ghost];
-                                mat_row[idx_ghost] = 0.0;
+                    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_FIXDIFF)
+                            continue;
+
+                        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)(bnd_len[0] + ctx->fd_stencil); bnd_idx++) {
+                                dst[-dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]] += dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]];
+                                dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]] = 0.0;
                             }
                     }
 
                     // fixval
-                    if (boundary_y_lower && ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXVAL) {
-                        const MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_1L];
-                        for (ptrdiff_t idx1_bnd = 1; idx1_bnd < ctx->fd_stencil; idx1_bnd++)
-                            for (ptrdiff_t idx0_col = -((ptrdiff_t)ctx->fd_stencil - 1); idx0_col < (ptrdiff_t)(ctx->domain_size[0] + ctx->fd_stencil - 1); idx0_col++) {
-                                const ptrdiff_t idx_ghost = -idx1_bnd * row_stride + idx0_col;
-                                e->rhs[mat_row_idx] -= bnd->val[idx1_bnd * bnd->val_stride + idx0_col] * mat_row[idx_ghost];
-                                mat_row[idx_ghost] = 0.0;
-                            }
-                    }
-                    if (boundary_y_upper && ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_BC_TYPE_FIXVAL) {
-                        const MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_1U];
-                        for (ptrdiff_t idx1_bnd = 1; idx1_bnd < ctx->fd_stencil; idx1_bnd++)
-                            for (ptrdiff_t idx0_col = -((ptrdiff_t)ctx->fd_stencil - 1); idx0_col < (ptrdiff_t)(ctx->domain_size[0] + ctx->fd_stencil - 1); idx0_col++) {
-                                const ptrdiff_t idx_ghost = (ctx->domain_size[1] - 1 + idx1_bnd) * row_stride + idx0_col;
-                                e->rhs[mat_row_idx] -= bnd->val[idx1_bnd * bnd->val_stride + idx0_col] * mat_row[idx_ghost];
-                                mat_row[idx_ghost] = 0.0;
-                            }
-                    }
-                    if (boundary_x_lower && ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXVAL) {
-                        const MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_0L];
-                        for (ptrdiff_t idx1_col = -((ptrdiff_t)ctx->fd_stencil - 1); idx1_col < (ptrdiff_t)(ctx->domain_size[1] + ctx->fd_stencil - 1); idx1_col++)
-                            for (ptrdiff_t idx0_bnd = 1; idx0_bnd < ctx->fd_stencil; idx0_bnd++) {
-                                const ptrdiff_t idx_ghost = idx1_col * row_stride - idx0_bnd;
-                                e->rhs[mat_row_idx] -= bnd->val[idx0_bnd * bnd->val_stride + idx1_col] * mat_row[idx_ghost];
-                                mat_row[idx1_col * row_stride - idx0_bnd] = 0.0;
-                            }
-                    }
-                    if (boundary_x_upper && ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXVAL) {
-                        const MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_0U];
-                        for (ptrdiff_t idx1_col = -((ptrdiff_t)ctx->fd_stencil - 1); idx1_col < (ptrdiff_t)(ctx->domain_size[1] + ctx->fd_stencil - 1); idx1_col++)
-                            for (ptrdiff_t idx0_bnd = 1; idx0_bnd < ctx->fd_stencil; idx0_bnd++) {
-                                const ptrdiff_t idx_ghost = idx1_col * row_stride + ctx->domain_size[0] - 1 + idx0_bnd;
-                                e->rhs[mat_row_idx] -= bnd->val[idx0_bnd * bnd->val_stride + idx1_col] * mat_row[idx_ghost];
-                                mat_row[idx_ghost] = 0.0;
+                    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_FIXVAL)
+                            continue;
+
+                        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)(bnd_len[0] + ctx->fd_stencil); bnd_idx++) {
+                                e->rhs[mat_row_idx] -= bnd->val[bnd_layer * bnd->val_stride + bnd_idx] * dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]];
+                                dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]] = 0.0;
                             }
                     }
 
@@ -637,14 +544,15 @@ int mg2di_egs_init(EGSContext *ctx)
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
         MG2DBoundary *bnd = ctx->boundaries[i];
+        const int ci  = mg2d_bnd_coord_idx(i);
         if (bnd->type == MG2D_BC_TYPE_FIXDIFF) {
-           for (int k = 0; k < ctx->domain_size[boundary_def[i].stride_idx]; k++)
+           for (int k = 0; k < ctx->domain_size[!ci]; k++)
                if (bnd->val[k] != 0.0) {
                    mg2di_log(&ctx->logger, 0, "Only zero boundary derivative supported\n");
                    return -ENOSYS;
                }
         } else if (bnd->type == MG2D_BC_TYPE_FIXVAL) {
-            priv->residual_calc_size[!boundary_def[i].stride_idx]--;
+            priv->residual_calc_size[ci]--;
         }
     }
 
@@ -709,7 +617,8 @@ static int arrays_alloc(EGSContext *ctx, const size_t domain_size[2])
     priv->stride = stride;
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
-        ctx->boundaries[i] = mg2di_bc_alloc(domain_size[boundary_def[i].stride_idx]);
+        const int ci  = mg2d_bnd_coord_idx(i);
+        ctx->boundaries[i] = mg2di_bc_alloc(domain_size[!ci]);
         if (!ctx->boundaries[i])
             return -ENOMEM;
     }