mg2d: factor out the boundary condition-related API

9 files changed, 147 insertions, 153 deletions

diff --git a/ell_relax.c b/ell_relax.c
index 982a95a..4ee2747 100644
--- a/ell_relax.c
+++ b/ell_relax.c
@@ -28,6 +28,8 @@
 #include "cpu.h"
 #include "ell_relax.h"
 #include "log.h"
+#include "mg2d_boundary.h"
+#include "mg2d_boundary_internal.h"
 #include "mg2d_constants.h"
 
 static const double relax_factors[FD_STENCIL_MAX] = {
@@ -42,8 +44,6 @@ struct EllRelaxInternal {
     double *residual_base;
     double *diff_coeffs_base[MG2D_DIFF_COEFF_NB];
 
-    double *boundary_val_base[4];
-
     double *residual_max;
     size_t  residual_max_size;
 
@@ -298,25 +298,25 @@ static void boundaries_apply(EllRelaxContext *ctx)
         const ptrdiff_t dst_strides[] = { stride_boundary,
             (boundary_def[i].is_upper ? 1 : -1) * stride_offset };
 
-        switch (ctx->boundaries[i].type) {
+        switch (ctx->boundaries[i]->type) {
         case MG2D_BC_TYPE_FIXVAL:
-            boundaries_apply_fixval(dst, dst_strides, ctx->boundaries[i].val,
-                                    ctx->boundaries[i].val_stride, size_boundary);
+            boundaries_apply_fixval(dst, dst_strides, ctx->boundaries[i]->val,
+                                    ctx->boundaries[i]->val_stride, size_boundary);
             break;
         case MG2D_BC_TYPE_FIXDIFF:
-            boundaries_apply_fixdiff(dst, dst_strides, ctx->boundaries[i].val,
-                                     ctx->boundaries[i].val_stride, size_boundary);
+            boundaries_apply_fixdiff(dst, dst_strides, ctx->boundaries[i]->val,
+                                     ctx->boundaries[i]->val_stride, size_boundary);
             break;
         }
     }
 
     /* 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) {
+    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_0L].type == MG2D_BC_TYPE_FIXDIFF)
+        if (ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXDIFF)
             fact_z *= -1;
         else
             fact_x *= -1;
@@ -328,12 +328,12 @@ static void boundaries_apply(EllRelaxContext *ctx)
                 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) {
+    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 + ctx->domain_size[0] - 1;
         int fact_x = 1, fact_z = -1;
 
-        if (ctx->boundaries[MG2D_BOUNDARY_0L].type == MG2D_BC_TYPE_FIXDIFF)
+        if (ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXDIFF)
             fact_z *= -1;
         else
             fact_x *= -1;
@@ -345,12 +345,12 @@ static void boundaries_apply(EllRelaxContext *ctx)
                 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) {
+    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 + strides[1] * (ctx->domain_size[1] - 1);
         int fact_x = -1, fact_z = 1;
 
-        if (ctx->boundaries[MG2D_BOUNDARY_0U].type == MG2D_BC_TYPE_FIXDIFF)
+        if (ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXDIFF)
             fact_z *= -1;
         else
             fact_x *= -1;
@@ -362,12 +362,12 @@ static void boundaries_apply(EllRelaxContext *ctx)
                 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) {
+    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_0U].type == MG2D_BC_TYPE_FIXDIFF)
+        if (ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXDIFF)
             fact_z *= -1;
         else
             fact_x *= -1;
@@ -475,13 +475,13 @@ int mg2di_ell_relax_init(EllRelaxContext *ctx)
     priv->residual_calc_offset  = 0;
     priv->residual_calc_size[0] = ctx->domain_size[0];
     priv->residual_calc_size[1] = ctx->domain_size[1];
-    if (ctx->boundaries[MG2D_BOUNDARY_0L].type == MG2D_BC_TYPE_FIXVAL)
+    if (ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXVAL)
         priv->residual_calc_offset += ctx->residual_stride;
-    if (ctx->boundaries[MG2D_BOUNDARY_1L].type == MG2D_BC_TYPE_FIXVAL)
+    if (ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXVAL)
         priv->residual_calc_offset++;
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
-        EllRelaxBoundary *bnd = &ctx->boundaries[i];
+        MG2DBoundary *bnd = ctx->boundaries[i];
         if (bnd->type == MG2D_BC_TYPE_FIXDIFF) {
            for (int k = 0; k < ctx->domain_size[boundary_def[i].stride_idx]; k++)
                if (bnd->val[k] != 0.0) {
@@ -553,16 +553,9 @@ static int ell_relax_arrays_alloc(EllRelaxContext *ctx, const size_t domain_size
     priv->stride = stride;
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
-        const size_t boundary_len0    = domain_size[boundary_def[i].stride_idx];
-        const size_t boundary_len_max = boundary_len0 + 2 * (FD_STENCIL_MAX - 1);
-
-        ret = posix_memalign((void**)&priv->boundary_val_base[i], 32,
-                             sizeof(*priv->boundary_val_base[i]) * boundary_len_max * FD_STENCIL_MAX);
-        if (ret != 0)
-            return -ret;
-
-        ctx->boundaries[i].val        = priv->boundary_val_base[i] + FD_STENCIL_MAX - 1;
-        ctx->boundaries[i].val_stride = boundary_len_max;
+        ctx->boundaries[i] = mg2di_bc_alloc(domain_size[boundary_def[i].stride_idx]);
+        if (!ctx->boundaries[i])
+            return -ENOMEM;
     }
 
     return 0;
@@ -613,8 +606,8 @@ void mg2di_ell_relax_free(EllRelaxContext **pctx)
     free(ctx->priv->residual_max);
     for (int i = 0; i < ARRAY_ELEMS(ctx->priv->diff_coeffs_base); i++)
         free(ctx->priv->diff_coeffs_base[i]);
-    for (int i = 0; i < ARRAY_ELEMS(ctx->priv->boundary_val_base); i++)
-        free(ctx->priv->boundary_val_base[i]);
+    for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++)
+        mg2di_bc_free(&ctx->boundaries[i]);
 
     tp_free(&ctx->priv->tp_internal);
     free(ctx->priv);
diff --git a/ell_relax.h b/ell_relax.h
index bc5b5fb..6a9ff2c 100644
--- a/ell_relax.h
+++ b/ell_relax.h
@@ -41,37 +41,11 @@
 #include <threadpool.h>
 
 #include "log.h"
+#include "mg2d_boundary.h"
 #include "mg2d_constants.h"
 
 typedef struct EllRelaxInternal EllRelaxInternal;
 
-/**
- * Boundary condition definition.
- */
-typedef struct EllRelaxBoundary {
-    /**
-     * Type of the boundary condition.
-     */
-    enum MG2DBCType type;
-    /**
-     * For type = ELL_RELAX_BC_TYPE_FIXVAL:
-     *     Values of the unknown function on the boundary.
-     *     The number of boundary layers is equal to fd_stencil.
-     *     The first boundary layer has the number of points equal to the
-     *     corresponding domain_size. Each subsequent boundary layer has one
-     *     more boundary point at each end of the domain.
-     *
-     * Ignored otherwise.
-     */
-    double *val;
-    /**
-     * Number of elements between rows in val. I.e. if val[0] is the first
-     * boundary point, then val[val_stride - 1] is the first boundary point in
-     * the second row and so on.
-     */
-    size_t  val_stride;
-} EllRelaxBoundary;
-
 typedef struct EllRelaxContext {
     /**
      * Solver private data, not to be accessed in any way by the caller.
@@ -113,10 +87,10 @@ typedef struct EllRelaxContext {
     size_t fd_stencil;
 
     /**
-     * Boundary specification, indexed by EllRelaxBoundaryLoc.
+     * Boundary specification, indexed by MG2DBoundaryLoc.
      * To be filled by the caller before mg2di_ell_relax_init().
      */
-    EllRelaxBoundary boundaries[4];
+    MG2DBoundary *boundaries[4];
 
     /**
      * The time stepping factor in relaxation.
diff --git a/meson.build b/meson.build
index c75f204..edd3b6f 100644
--- a/meson.build
+++ b/meson.build
@@ -4,6 +4,7 @@ project('libmg2d', 'c',
 add_project_arguments('-D_XOPEN_SOURCE=700', language : 'c')
 
 lib_src = [
+    'boundary.c',
     'cpu.c',
     'ell_relax.c',
     'log.c',
diff --git a/mg2d.c b/mg2d.c
index af33450..3406030 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -30,6 +30,8 @@
 #include "ell_relax.h"
 #include "log.h"
 #include "mg2d.h"
+#include "mg2d_boundary.h"
+#include "mg2d_boundary_internal.h"
 #include "mg2d_constants.h"
 
 #define LEVELS_MAX 16
@@ -64,8 +66,6 @@ struct MG2DInternal {
 
     int64_t time_solve;
     int64_t count_solve;
-
-    double *boundaries_base[4];
 };
 
 static void log_callback(MG2DLogger *log, int level, const char *fmt, va_list vl)
@@ -374,7 +374,7 @@ finish:
     return 0;
 }
 
-static void bnd_zero(EllRelaxBoundary *bdst, size_t nb_rows, size_t domain_size)
+static void bnd_zero(MG2DBoundary *bdst, size_t nb_rows, size_t domain_size)
 {
     for (size_t i = 0; i < nb_rows; i++) {
         memset(bdst->val + i * bdst->val_stride - i, 0,
@@ -382,7 +382,7 @@ static void bnd_zero(EllRelaxBoundary *bdst, size_t nb_rows, size_t domain_size)
     }
 }
 
-static void bnd_copy(EllRelaxBoundary *bdst, const double *src, ptrdiff_t src_stride,
+static void bnd_copy(MG2DBoundary *bdst, const double *src, ptrdiff_t src_stride,
                      size_t nb_rows, size_t domain_size)
 {
     for (size_t i = 0; i < nb_rows; i++) {
@@ -431,7 +431,7 @@ static int mg_levels_init(MG2DContext *ctx)
 
         for (int i = 0; i < ARRAY_ELEMS(cur->solver->boundaries); i++) {
             MG2DBoundary *bsrc = ctx->boundaries[i];
-            EllRelaxBoundary *bdst = &cur->solver->boundaries[i];
+            MG2DBoundary *bdst = cur->solver->boundaries[i];
             bdst->type = bsrc->type;
             switch (bsrc->type) {
             case MG2D_BC_TYPE_FIXVAL:
@@ -647,18 +647,11 @@ MG2DContext *mg2d_solver_alloc(size_t domain_size)
     priv->cpuflags      = mg2di_cpu_flags_get();
 
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
-        const size_t boundary_len_max = domain_size + 2 * (FD_STENCIL_MAX - 1);
-
-        ctx->boundaries[i] = calloc(1, sizeof(*ctx->boundaries[i]));
-        if (!ctx->boundaries[i])
-            goto fail;
-
-        ret = posix_memalign((void**)&priv->boundaries_base[i], 32,
-                             sizeof(*ctx->boundaries[i]->val) * boundary_len_max * FD_STENCIL_MAX);
-        if (ret != 0)
+        ctx->boundaries[i] = mg2di_bc_alloc(domain_size);
+        if (!ctx->boundaries[i]) {
+            ret = -ENOMEM;
             goto fail;
-        ctx->boundaries[i]->val        = priv->boundaries_base[i] + FD_STENCIL_MAX - 1;
-        ctx->boundaries[i]->val_stride = boundary_len_max;
+        }
     }
 
     ret = mg_levels_alloc(ctx, domain_size);
@@ -707,10 +700,8 @@ void mg2d_solver_free(MG2DContext **pctx)
         ctx->priv->root = next;
     }
 
-    for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
-        free(ctx->priv->boundaries_base[i]);
-        free(ctx->boundaries[i]);
-    }
+    for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++)
+        mg2di_bc_free(&ctx->boundaries[i]);
 
     tp_free(&ctx->priv->tp);
     free(ctx->priv);
diff --git a/mg2d.h b/mg2d.h
index 6b2def6..234d8cf 100644
--- a/mg2d.h
+++ b/mg2d.h
@@ -22,37 +22,11 @@
 #include <stddef.h>
 #include <stdint.h>
 
+#include "mg2d_boundary.h"
 #include "mg2d_constants.h"
 
 typedef struct MG2DInternal MG2DInternal;
 
-/**
- * Boundary condition definition.
- */
-typedef struct MG2DBoundary {
-    /**
-     * Type of the boundary condition.
-     */
-    enum MG2DBCType type;
-    /**
-     * For type = MG2D_BC_TYPE_FIXVAL:
-     *     Values of the unknown function on the boundary.
-     *     The number of boundary layers is equal to fd_stencil.
-     *     The first boundary layer has the number of points equal to the
-     *     corresponding domain_size. Each subsequent boundary layer has one
-     *     more boundary point at each end of the domain.
-     *
-     * Ignored otherwise.
-     */
-    double *val;
-    /**
-     * Number of elements between rows in val. I.e. if val[0] is the first
-     * boundary point, then val[val_stride - 1] is the first boundary point in
-     * the second row and so on.
-     */
-    size_t  val_stride;
-} MG2DBoundary;
-
 typedef struct MG2DContext {
     /**
      * Solver private data, not to be accessed in any way by the caller.
diff --git a/mg2d_boundary.h b/mg2d_boundary.h
new file mode 100644
index 0000000..4a7dcee
--- /dev/null
+++ b/mg2d_boundary.h
@@ -0,0 +1,97 @@
+/*
+ * Boundary condition declaration.
+ * Copyright 2018 Anton Khirnov <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef MG2D_BOUNDARY_H
+#define MG2D_BOUNDARY_H
+
+#include <stddef.h>
+
+typedef struct MG2DBoundaryInternal MG2DBoundaryInternal;
+
+/**
+ * Type of the boundary condition on a given boundary.
+ */
+enum MG2DBCType {
+    /**
+     * The value of the unknown function is fixed on the boundary.
+     */
+    MG2D_BC_TYPE_FIXVAL,
+    /**
+     * The normal derivative of the unkown function is fixed on the boundary.
+     *
+     * TODO: the only supported value of the derivative is currently zero, i.e.
+     * periodic boundary conditions.
+     */
+    MG2D_BC_TYPE_FIXDIFF,
+};
+
+/**
+ * Location of the boundary.
+ */
+enum MG2DBoundaryLoc {
+    /**
+     * coord0 lower
+     */
+    MG2D_BOUNDARY_0L,
+    /**
+     * coord0 upper
+     */
+    MG2D_BOUNDARY_0U,
+    /**
+     * coord1 lower
+     */
+    MG2D_BOUNDARY_1L,
+    /**
+     * coord1 upper
+     */
+    MG2D_BOUNDARY_1U,
+};
+
+/**
+ * Boundary condition definition.
+ */
+typedef struct MG2DBoundary {
+    /**
+     * Type of the boundary condition.
+     */
+    enum MG2DBCType type;
+    /**
+     * For type = MG2D_BC_TYPE_FIXVAL:
+     *     Values of the unknown function on the boundary.
+     *     The number of boundary layers is equal to fd_stencil.
+     *     The first boundary layer has the number of points equal to the
+     *     corresponding domain_size. Each subsequent boundary layer has one
+     *     more boundary point at each end of the domain.
+     *
+     * Ignored otherwise.
+     */
+    double *val;
+    /**
+     * Number of elements between rows in val. I.e. if val[0] is the first
+     * boundary point, then val[val_stride - 1] is the first boundary point in
+     * the second row and so on.
+     */
+    size_t  val_stride;
+
+    /**
+     * Private data, not to be touched by callers.
+     */
+    MG2DBoundaryInternal *priv;
+} MG2DBoundary;
+
+#endif // MG2D_BOUNDARY_H
diff --git a/mg2d_constants.h b/mg2d_constants.h
index 1ab4bc6..aa3d5cb 100644
--- a/mg2d_constants.h
+++ b/mg2d_constants.h
@@ -24,45 +24,6 @@ enum MG2DError {
 };
 
 /**
- * Type of the boundary condition on a given boundary.
- */
-enum MG2DBCType {
-    /**
-     * The value of the unknown function is fixed on the boundary.
-     */
-    MG2D_BC_TYPE_FIXVAL,
-    /**
-     * The normal derivative of the unkown function is fixed on the boundary.
-     *
-     * TODO: the only supported value of the derivative is currently zero, i.e.
-     * periodic boundary conditions.
-     */
-    MG2D_BC_TYPE_FIXDIFF,
-};
-
-/**
- * Location of the boundary.
- */
-enum MG2DBoundaryLoc {
-    /**
-     * coord0 lower
-     */
-    MG2D_BOUNDARY_0L,
-    /**
-     * coord0 upper
-     */
-    MG2D_BOUNDARY_0U,
-    /**
-     * coord1 lower
-     */
-    MG2D_BOUNDARY_1L,
-    /**
-     * coord1 upper
-     */
-    MG2D_BOUNDARY_1U,
-};
-
-/**
  * Derivative specifier.
  */
 enum MG2DDiffCoeff {
diff --git a/mg2d_test.c b/mg2d_test.c
index c04e824..6bdbb00 100644
--- a/mg2d_test.c
+++ b/mg2d_test.c
@@ -5,6 +5,8 @@
 #include <string.h>
 
 #include "mg2d.h"
+#include "mg2d_boundary.h"
+#include "mg2d_constants.h"
 
 #define MAXITER 64
 #define TOL 1e-9
diff --git a/relax_test.c b/relax_test.c
index 1994201..e885d55 100644
--- a/relax_test.c
+++ b/relax_test.c
@@ -6,6 +6,7 @@
 
 #include "ell_relax.h"
 #include "log.h"
+#include "mg2d_boundary.h"
 #include "mg2d_constants.h"
 
 #if 1
@@ -84,7 +85,7 @@ int main(int argc, char **argv)
     ctx->fd_stencil = 2;
 
     {
-        EllRelaxBoundary *bnd = &ctx->boundaries[MG2D_BOUNDARY_0L];
+        MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_0L];
 
         bnd->type = MG2D_BC_TYPE_FIXVAL;
 
@@ -98,7 +99,7 @@ int main(int argc, char **argv)
         }
     }
     {
-        EllRelaxBoundary *bnd = &ctx->boundaries[MG2D_BOUNDARY_0U];
+        MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_0U];
 
         bnd->type = MG2D_BC_TYPE_FIXVAL;
 
@@ -112,7 +113,7 @@ int main(int argc, char **argv)
         }
     }
     {
-        EllRelaxBoundary *bnd = &ctx->boundaries[MG2D_BOUNDARY_1L];
+        MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_1L];
 
         bnd->type = MG2D_BC_TYPE_FIXVAL;
 
@@ -126,7 +127,7 @@ int main(int argc, char **argv)
         }
     }
     {
-        EllRelaxBoundary *bnd = &ctx->boundaries[MG2D_BOUNDARY_1U];
+        MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_1U];
 
         bnd->type = MG2D_BC_TYPE_FIXVAL;