From 763dfb979c06707d334532cb31a776fb67564ce7 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Wed, 1 Aug 2018 17:04:12 +0200 Subject: higher finite difference order --- ell_relax.c | 40 +++++++++++++++++++++++++++++++++++----- 1 file changed, 35 insertions(+), 5 deletions(-) diff --git a/ell_relax.c b/ell_relax.c index 9454e84..96038e8 100644 --- a/ell_relax.c +++ b/ell_relax.c @@ -25,7 +25,7 @@ #include "ell_relax.h" #include "log.h" -#define FD_STENCIL_MAX 1 +#define FD_STENCIL_MAX 2 struct EllRelaxInternal { ptrdiff_t stride; @@ -33,6 +33,8 @@ struct EllRelaxInternal { double *rhs_base; double *residual_base; double *diff_coeffs_base[ELL_RELAX_DIFF_COEFF_NB]; + + void (*derivatives_calc)(EllRelaxContext *ctx, double *dst, ptrdiff_t idx, ptrdiff_t stride); }; static const struct { @@ -57,8 +59,8 @@ static const struct { }, }; -static inline void -derivatives_calc(EllRelaxContext *ctx, double *dst, ptrdiff_t idx, ptrdiff_t stride) +static void +derivatives_calc_s1(EllRelaxContext *ctx, double *dst, ptrdiff_t idx, ptrdiff_t stride) { const double *u = ctx->u; const double dx = ctx->step[0]; @@ -74,6 +76,27 @@ derivatives_calc(EllRelaxContext *ctx, double *dst, ptrdiff_t idx, ptrdiff_t str dst[ELL_RELAX_DIFF_COEFF_11] = (u[idx + 1 + stride] - u[idx + stride - 1] - u[idx - stride + 1] + u[idx - stride - 1]) / (4.0 * dx * dy); } +static void +derivatives_calc_s2(EllRelaxContext *ctx, double *dst, ptrdiff_t idx, ptrdiff_t stride) +{ + const double *u = ctx->u; + const double dx = ctx->step[0]; + const double dy = ctx->step[1]; + + dst[ELL_RELAX_DIFF_COEFF_00] = u[idx]; + dst[ELL_RELAX_DIFF_COEFF_10] = (-1.0 * u[idx + 2] + 8.0 * u[idx + 1] - 8.0 * u[idx - 1] + 1.0 * u[idx - 2]) / (12.0 * dx); + dst[ELL_RELAX_DIFF_COEFF_01] = (-1.0 * u[idx + 2 * stride] + 8.0 * u[idx + stride] - 8.0 * u[idx - stride] + 1.0 * u[idx - 2 * stride]) / (12.0 * dy); + + dst[ELL_RELAX_DIFF_COEFF_20] = (-1.0 * u[idx + 2] + 16.0 * u[idx + 1] - 30.0 * u[idx] + 16.0 * u[idx - 1] - 1.0 * u[idx - 2]) / (12.0 * SQR(dx)); + dst[ELL_RELAX_DIFF_COEFF_02] = (-1.0 * u[idx + 2 * stride] + 16.0 * u[idx + stride] - 30.0 * u[idx] + 16.0 * u[idx - stride] - 1.0 * u[idx - 2 * stride]) / (12.0 * SQR(dy)); + + dst[ELL_RELAX_DIFF_COEFF_11] = ( 1.0 * u[idx + 2 + 2 * stride] - 8.0 * u[idx + 2 + stride] + 8.0 * u[idx + 2 - stride] - 1.0 * u[idx + 2 - 2 * stride] + -8.0 * u[idx + 1 + 2 * stride] + 64.0 * u[idx + 1 + stride] - 64.0 * u[idx + 1 - stride] + 8.0 * u[idx + 1 - 2 * stride] + +8.0 * u[idx - 1 + 2 * stride] - 64.0 * u[idx - 1 + stride] + 64.0 * u[idx - 1 - stride] - 8.0 * u[idx - 1 - 2 * stride] + -1.0 * u[idx - 2 + 2 * stride] + 8.0 * u[idx - 2 + stride] - 8.0 * u[idx - 2 - stride] + 1.0 * u[idx - 2 - 2 * stride]) + / (144.0 * dx * dy); +} + static void residual_calc(EllRelaxContext *ctx) { EllRelaxInternal *priv = ctx->priv; @@ -85,7 +108,7 @@ static void residual_calc(EllRelaxContext *ctx) double u[ELL_RELAX_DIFF_COEFF_NB]; double res; - derivatives_calc(ctx, u, idx, priv->stride); + priv->derivatives_calc(ctx, u, idx, priv->stride); res = -ctx->rhs[idx]; for (int i = 0; i < ARRAY_ELEMS(u); i++) @@ -229,7 +252,14 @@ int mg2di_ell_relax_init(EllRelaxContext *ctx) EllRelaxInternal *priv = ctx->priv; int ret; - if (ctx->fd_stencil == 0 || ctx->fd_stencil > FD_STENCIL_MAX) { + switch (ctx->fd_stencil) { + case 1: + priv->derivatives_calc = derivatives_calc_s1; + break; + case 2: + priv->derivatives_calc = derivatives_calc_s2; + break; + default: mg2di_log(&ctx->logger, 0, "Invalid finite difference stencil: %zd\n", ctx->fd_stencil); return -EINVAL; -- cgit v1.2.3