From d76687285a032e25cbd258035321cb8d6d8e0330 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 25 Jan 2019 16:16:17 +0100 Subject: mg2d: factor out the boundary condition-related API --- ell_relax.c | 61 ++++++++++++++++------------------- ell_relax.h | 32 ++----------------- meson.build | 1 + mg2d.c | 31 +++++++----------- mg2d.h | 28 +--------------- mg2d_boundary.h | 97 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ mg2d_constants.h | 39 ----------------------- mg2d_test.c | 2 ++ relax_test.c | 9 +++--- 9 files changed, 147 insertions(+), 153 deletions(-) create mode 100644 mg2d_boundary.h 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 #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 #include +#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 + * + * 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 . + */ + +#ifndef MG2D_BOUNDARY_H +#define MG2D_BOUNDARY_H + +#include + +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 @@ -23,45 +23,6 @@ enum MG2DError { MG2D_ERR_DIVERGE = -0xff00, }; -/** - * 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. */ 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 #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; -- cgit v1.2.3