From 7f4b62ef082ee18090d00b3cad3bce3c89f5dc71 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sun, 3 Mar 2019 14:04:35 +0100 Subject: Implement C/r falloff boundary condition. API bump. --- ell_grid_solve.c | 70 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ libmg2d.v | 2 +- meson.build | 2 +- mg2d.c | 1 + mg2d_boundary.h | 5 ++++ 5 files changed, 78 insertions(+), 2 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++) { diff --git a/libmg2d.v b/libmg2d.v index f4dd528..c01a69f 100644 --- a/libmg2d.v +++ b/libmg2d.v @@ -1,4 +1,4 @@ -LIBMG2D_9 { +LIBMG2D_10 { global: mg2d_*; local: *; }; diff --git a/meson.build b/meson.build index 3e639c7..48d80ef 100644 --- a/meson.build +++ b/meson.build @@ -22,7 +22,7 @@ liblapacke = cc.find_library('lapacke') dep_tp = declare_dependency(link_args : '-lthreadpool') -deps = [dep_tp, liblapacke] +deps = [dep_tp, libm, liblapacke] cdata = configuration_data() cdata.set10('ARCH_X86', false) diff --git a/mg2d.c b/mg2d.c index 07e162e..0a17793 100644 --- a/mg2d.c +++ b/mg2d.c @@ -457,6 +457,7 @@ static int mg_levels_init(MG2DContext *ctx) bnd_zero(bdst, ctx->fd_stencil, cur->solver->domain_size[0]); break; case MG2D_BC_TYPE_REFLECT: + case MG2D_BC_TYPE_FALLOFF: bnd_zero(bdst, ctx->fd_stencil, cur->solver->domain_size[0]); break; default: diff --git a/mg2d_boundary.h b/mg2d_boundary.h index 910249f..752661d 100644 --- a/mg2d_boundary.h +++ b/mg2d_boundary.h @@ -35,6 +35,11 @@ enum MG2DBCType { * The unknown function is mirrored across the boundary */ MG2D_BC_TYPE_REFLECT, + + /** + * The unknown function falls off as C / x on the boundary. + */ + MG2D_BC_TYPE_FALLOFF, }; /** -- cgit v1.2.3