From be15318094225dc75d1f3c4a7984e770d4e128ec Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 25 Jan 2019 14:11:01 +0100 Subject: ell_relax: implement correct handling of the boundary corners --- ell_relax.c | 60 ++++++++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 48 insertions(+), 12 deletions(-) diff --git a/ell_relax.c b/ell_relax.c index 57c4112..891fc27 100644 --- a/ell_relax.c +++ b/ell_relax.c @@ -310,36 +310,72 @@ static void boundaries_apply(EllRelaxContext *ctx) } /* fill in the corner ghosts */ - { + if (ctx->boundaries[ELL_RELAX_BOUNDARY_0L].type == ELL_RELAX_BC_TYPE_FIXDIFF || + ctx->boundaries[ELL_RELAX_BOUNDARY_1L].type == ELL_RELAX_BC_TYPE_FIXDIFF) { double *dst = ctx->u; + int fact_x = -1, fact_z = -1; + + if (ctx->boundaries[ELL_RELAX_BOUNDARY_0L].type == ELL_RELAX_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 = j * strides[1] + i; - dst[-idx] = dst[idx]; + 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[ELL_RELAX_BOUNDARY_0L].type == ELL_RELAX_BC_TYPE_FIXDIFF || + ctx->boundaries[ELL_RELAX_BOUNDARY_1U].type == ELL_RELAX_BC_TYPE_FIXDIFF) { double *dst = ctx->u + ctx->domain_size[0] - 1; + int fact_x = 1, fact_z = -1; + + if (ctx->boundaries[ELL_RELAX_BOUNDARY_0L].type == ELL_RELAX_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 = j * strides[1] - i; - dst[-idx] = dst[idx]; + 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[ELL_RELAX_BOUNDARY_1L].type == ELL_RELAX_BC_TYPE_FIXDIFF || + ctx->boundaries[ELL_RELAX_BOUNDARY_0U].type == ELL_RELAX_BC_TYPE_FIXDIFF) { double *dst = ctx->u + strides[1] * (ctx->domain_size[1] - 1); + int fact_x = -1, fact_z = 1; + + if (ctx->boundaries[ELL_RELAX_BOUNDARY_0U].type == ELL_RELAX_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 = -j * strides[1] + i; - dst[-idx] = dst[idx]; + 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[ELL_RELAX_BOUNDARY_0U].type == ELL_RELAX_BC_TYPE_FIXDIFF || + ctx->boundaries[ELL_RELAX_BOUNDARY_1U].type == ELL_RELAX_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[ELL_RELAX_BOUNDARY_0U].type == ELL_RELAX_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 = -j * strides[1] - i; - dst[-idx] = dst[idx]; + 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]; } } ctx->time_boundaries += gettime() - start; -- cgit v1.2.3