From a1f10a1338e77dd147e7c94f3f401d2703181bca Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 14 Jun 2019 15:29:02 +0200 Subject: rescalc: improve reflection boundary conditions Make parameter names more clear/consistent, document them, implement missing 1U boundary. --- ell_grid_solve.c | 20 +++++++++++--------- residual_calc.c | 39 +++++++++++++++++++++++++-------------- residual_calc.h | 9 ++++++++- 3 files changed, 44 insertions(+), 24 deletions(-) diff --git a/ell_grid_solve.c b/ell_grid_solve.c index 6f1a029..5047718 100644 --- a/ell_grid_solve.c +++ b/ell_grid_solve.c @@ -152,12 +152,18 @@ static void residual_calc(EGSContext *ctx, int export_res) { EGSInternal *priv = ctx->priv; const double *diff_coeffs[MG2D_DIFF_COEFF_NB]; + int reflect_flags = 0; mg2di_timer_start(&ctx->timer_res_calc); for (int i = 0; i < ARRAY_ELEMS(diff_coeffs); i++) diff_coeffs[i] = NDPTR2D(priv->diff_coeffs[i], priv->residual_calc_offset[1], priv->residual_calc_offset[0]); + for (int bnd_idx = 0; bnd_idx < 4; bnd_idx++) + if (ctx->boundaries[bnd_idx]->type == MG2D_BC_TYPE_REFLECT && + priv->dg->components[priv->local_component].bnd_is_outer[bnd_idx]) + reflect_flags |= 1 << bnd_idx; + mg2di_residual_calc(priv->rescalc, priv->residual_calc_size, &ctx->residual_max, NDPTR2D(export_res ? ctx->residual : priv->u_next, priv->residual_calc_offset[1], priv->residual_calc_offset[0]), export_res ? ctx->residual->stride[0] : priv->u_next->stride[0], @@ -167,18 +173,14 @@ static void residual_calc(EGSContext *ctx, int export_res) ctx->rhs->stride[0], diff_coeffs, priv->diff_coeffs[0]->stride[0], export_res ? 0.0 : 1.0, export_res ? 1.0 : priv->r.relax_factor, - export_res ? 0 : - ((ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_REFLECT) << 0) | - ((ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_REFLECT) << 1) | - ((ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_REFLECT) << 2)); + reflect_flags, FD_STENCIL_MAX); mg2di_timer_stop(&ctx->timer_res_calc); - if (!export_res) { - priv->reflect_disable[MG2D_BOUNDARY_0L] = 1; - priv->reflect_disable[MG2D_BOUNDARY_0U] = 1; - priv->reflect_disable[MG2D_BOUNDARY_1L] = 1; - } + priv->reflect_disable[MG2D_BOUNDARY_0L] = 1; + priv->reflect_disable[MG2D_BOUNDARY_0U] = 1; + priv->reflect_disable[MG2D_BOUNDARY_1L] = 1; + priv->reflect_disable[MG2D_BOUNDARY_1U] = 1; } static void boundaries_apply_fixval(double *dst, const ptrdiff_t dst_stride[2], diff --git a/residual_calc.c b/residual_calc.c index 3b83c63..6d58c10 100644 --- a/residual_calc.c +++ b/residual_calc.c @@ -27,11 +27,12 @@ #include "common.h" #include "cpu.h" +#include "mg2d_boundary.h" #include "mg2d_constants.h" #include "residual_calc.h" typedef struct ResidualCalcTask { - size_t line_size; + size_t size[2]; double *dst; ptrdiff_t dst_stride; @@ -48,7 +49,8 @@ typedef struct ResidualCalcTask { double u_mult; double res_mult; - int mirror_line; + int reflect; + size_t reflect_dist; } ResidualCalcTask; struct ResidualCalcInternal { @@ -255,28 +257,34 @@ static int residual_calc_task(void *arg, unsigned int job_idx, unsigned int thre diff_coeffs[i] = task->diff_coeffs[i] + job_idx * task->diff_coeffs_stride; if (task->u_mult == 0.0) { - priv->residual_calc_line(task->line_size, dst, + priv->residual_calc_line(task->size[0], dst, priv->residual_max + thread_idx * priv->calc_blocksize, task->u_stride, task->u + job_idx * task->u_stride, task->rhs + job_idx * task->rhs_stride, diff_coeffs, task->res_mult); } else { - priv->residual_add_line(task->line_size, dst, + priv->residual_add_line(task->size[0], dst, priv->residual_max + thread_idx * priv->calc_blocksize, task->u_stride, task->u + job_idx * task->u_stride, task->rhs + job_idx * task->rhs_stride, diff_coeffs, task->res_mult, task->u_mult); } - if (task->mirror_line & (1 << 0)) { - for (int i = 1; i <= FD_STENCIL_MAX; i++) + + if (task->reflect & (1 << MG2D_BOUNDARY_0L)) { + for (int i = 1; i <= task->reflect_dist; i++) dst[-i] = dst[i]; } - if (task->mirror_line & (1 << 1)) { - for (int i = 1; i <= FD_STENCIL_MAX; i++) - dst[task->line_size - 1 + i] = dst[task->line_size - 1 - i]; + if (task->reflect & (1 << MG2D_BOUNDARY_0U)) { + for (int i = 1; i <= task->reflect_dist; i++) + dst[task->size[0] - 1 + i] = dst[task->size[0] - 1 - i]; + } + if ((task->reflect & (1 << MG2D_BOUNDARY_1L)) && + job_idx > 0 && job_idx <= task->reflect_dist) { + memcpy(task->dst - job_idx * task->dst_stride, dst, sizeof(*dst) * task->size[0]); } - if ((task->mirror_line & (1 << 2)) && job_idx > 0 && job_idx <= FD_STENCIL_MAX) { - memcpy(task->dst - job_idx * task->dst_stride, dst, sizeof(*dst) * task->line_size); + if ((task->reflect & (1 << MG2D_BOUNDARY_1U)) && + job_idx >= task->size[1] - 1 - task->reflect_dist && job_idx < task->size[1] - 1) { + memcpy(task->dst + (2 * (task->size[1] - 1) - job_idx) * task->dst_stride, dst, sizeof(*dst) * task->size[0]); } return 0; @@ -289,7 +297,8 @@ int mg2di_residual_calc(ResidualCalcContext *ctx, size_t size[2], const double *rhs, ptrdiff_t rhs_stride, const double * const diff_coeffs[MG2D_DIFF_COEFF_NB], ptrdiff_t diff_coeffs_stride, - double u_mult, double res_mult, int mirror_line) + double u_mult, double res_mult, + int reflect, size_t reflect_dist) { ResidualCalcInternal *priv = ctx->priv; ResidualCalcTask *task = &priv->task; @@ -297,7 +306,8 @@ int mg2di_residual_calc(ResidualCalcContext *ctx, size_t size[2], memset(priv->residual_max, 0, sizeof(*priv->residual_max) * priv->residual_max_size); - task->line_size = size[0]; + task->size[0] = size[0]; + task->size[1] = size[1]; task->dst = dst; task->dst_stride = dst_stride; task->u = u; @@ -308,7 +318,8 @@ int mg2di_residual_calc(ResidualCalcContext *ctx, size_t size[2], task->diff_coeffs_stride = diff_coeffs_stride; task->u_mult = u_mult; task->res_mult = res_mult; - task->mirror_line = mirror_line; + task->reflect = reflect; + task->reflect_dist = reflect_dist; tp_execute(ctx->tp, size[1], residual_calc_task, priv); diff --git a/residual_calc.h b/residual_calc.h index 8c9b628..1498bc9 100644 --- a/residual_calc.h +++ b/residual_calc.h @@ -67,6 +67,12 @@ int mg2di_residual_calc_init(ResidualCalcContext *ctx); /** * Calculate the residual. + * + * @param reflect Flags indicating if dst should be extended by reflection. + * If the 'MG2D_BOUNDARY_x'th bit is set in reflect, then dst should be + * reflected reflect_dist points beyond its bounds in the corresponding + * direction. + * @param reflect_dist number of points to fill by reflection */ int mg2di_residual_calc(ResidualCalcContext *ctx, size_t size[2], double *residual_max, @@ -75,6 +81,7 @@ int mg2di_residual_calc(ResidualCalcContext *ctx, size_t size[2], const double *rhs, ptrdiff_t rhs_stride, const double * const diff_coeffs[MG2D_DIFF_COEFF_NB], ptrdiff_t diff_coeffs_stride, - double u_mult, double res_mult, int mirror_line); + double u_mult, double res_mult, + int reflect, size_t reflect_dist); #endif // MG2D_RESIDUAL_CALC_H -- cgit v1.2.3