From 373dbd15f216c2ce97396af37ea3f7cc5d39455c Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Tue, 9 Apr 2019 09:45:16 +0200
Subject: egs: add higher-order finite difference operators

---
 residual_calc_template.c | 90 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)
 create mode 100644 residual_calc_template.c

(limited to 'residual_calc_template.c')

diff --git a/residual_calc_template.c b/residual_calc_template.c
new file mode 100644
index 0000000..dcf7eb4
--- /dev/null
+++ b/residual_calc_template.c
@@ -0,0 +1,90 @@
+/*
+ * Residual calculation template
+ * Copyright 2019 Anton Khirnov <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#define JOIN2(a, b) a ## _ ## b
+#define FUNC2(name, stencil) JOIN2(name, stencil)
+#define FUNC(name) FUNC2(name, STENCIL)
+
+static void
+FUNC(derivatives_calc)(double *dst, const double *u, const double *fd_factors, ptrdiff_t stride)
+{
+    dst[MG2D_DIFF_COEFF_00] = u[0];
+    
+    {
+        double val = 0.0;
+        for (int i = -STENCIL; i <= STENCIL; i++)
+            val += fd_coeffs[STENCIL - 1][0][i + STENCIL] * u[i];
+        dst[MG2D_DIFF_COEFF_10] = val * fd_factors[MG2D_DIFF_COEFF_10];
+    }
+    {
+        double val = 0.0;
+        for (int i = -STENCIL; i <= STENCIL; i++)
+            val += fd_coeffs[STENCIL - 1][0][i + STENCIL] * u[i * stride];
+        dst[MG2D_DIFF_COEFF_01] = val * fd_factors[MG2D_DIFF_COEFF_01];
+    }
+    {
+        double val = 0.0;
+        for (int i = -STENCIL; i <= STENCIL; i++)
+            val += fd_coeffs[STENCIL - 1][1][i + STENCIL] * u[i];
+        dst[MG2D_DIFF_COEFF_20] = val * fd_factors[MG2D_DIFF_COEFF_20];
+    }
+    {
+        double val = 0.0;
+        for (int i = -STENCIL; i <= STENCIL; i++)
+            val += fd_coeffs[STENCIL - 1][1][i + STENCIL] * u[i * stride];
+        dst[MG2D_DIFF_COEFF_02] = val * fd_factors[MG2D_DIFF_COEFF_02];
+    }
+    {
+        double val = 0.0;
+        for (int j = -STENCIL; j <= STENCIL; j++) {
+            double tmp = 0.0;
+            for (int i = -STENCIL; i <= STENCIL; i++)
+                tmp += fd_coeffs[STENCIL - 1][0][i + STENCIL] * u[j * stride + i];
+            val += fd_coeffs[STENCIL - 1][0][j + STENCIL] * tmp;
+        }
+        dst[MG2D_DIFF_COEFF_11] = val * fd_factors[MG2D_DIFF_COEFF_11];
+    }
+}
+static void
+FUNC(residual_calc_line)(size_t linesize, double *dst, double *dst_max,
+                         ptrdiff_t stride, const double *u, const double *rhs,
+                         const double * const diff_coeffs[MG2D_DIFF_COEFF_NB],
+                         const double *fd_factors)
+{
+    double res_max = 0.0, res_abs;
+    for (size_t i = 0; i < linesize; i++) {
+        double u_vals[MG2D_DIFF_COEFF_NB];
+        double res;
+
+        FUNC(derivatives_calc)(u_vals, u + i, fd_factors, stride);
+
+        res = -rhs[i];
+        for (int j = 0; j < ARRAY_ELEMS(u_vals); j++)
+            res += u_vals[j] * diff_coeffs[j][i];
+        dst[i] = res;
+
+        res_abs = fabs(res);
+        res_max = MAX(res_max, res_abs);
+    }
+
+    *dst_max = MAX(*dst_max, res_max);
+}
+
+#undef FUNC
+#undef FUNC2
+#undef JOIN2
-- 
cgit v1.2.3