aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAnton Khirnov <anton@khirnov.net>2019-06-14 15:29:02 +0200
committerAnton Khirnov <anton@khirnov.net>2019-06-14 16:10:14 +0200
commita1f10a1338e77dd147e7c94f3f401d2703181bca (patch)
tree6ee70770635391904b4146a4ed0a9f6c691609cd
parentd61fcbfbdb1517c39d6d894f687f087a8c728b27 (diff)
rescalc: improve reflection boundary conditions
Make parameter names more clear/consistent, document them, implement missing 1U boundary.
-rw-r--r--ell_grid_solve.c20
-rw-r--r--residual_calc.c39
-rw-r--r--residual_calc.h9
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