higher finite difference order

1 files 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;