eval: refactor evaluating the metric

1 files changed, 68 insertions, 0 deletions

diff --git a/eval_metric_template.c b/eval_metric_template.c
new file mode 100644
index 0000000..8443785
--- /dev/null
+++ b/eval_metric_template.c
@@ -0,0 +1,68 @@
+/*
+ * Copyright 2016 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/>.
+ */
+
+/**
+ * @file
+ * the template for the loop over the grid points
+*/
+
+#define FUNC3(a, b, c) a ## _ ## b ## _ ## c
+#define FUNC2(a, b, c) FUNC3(a, b, c)
+#define FUNC(x) FUNC2(x, COMP, DIFF)
+
+#define DIFF_RHO (DIFF % 3)
+#define DIFF_Z   (DIFF / 3)
+
+static void FUNC(eval_metric)(const BDContext *bd, double *out, ptrdiff_t out_stride,
+                              int nb_coords_rho, const double *rho,
+                              int nb_coords_z,   const double *z,
+                              double *psi[9], double *q[9])
+{
+    BDPriv *s = bd->priv;
+
+    for (int i = 0; i < nb_coords_z; i++) {
+        const double zz = z[i];
+        double *dst = out + i * out_stride;
+
+        for (int j = 0; j < nb_coords_rho; j++) {
+            const int idx = i * nb_coords_rho + j;
+            const double rrho = rho[j];
+
+            const double ppsi = psi[0][idx];
+            const double psi2 = SQR(ppsi);
+            const double psi3 = psi2 * ppsi;
+            const double psi4 = SQR(psi2);
+
+            const double base     = psi4 * exp(2 * q[0][idx]);
+            const double dqrho    = DIFF_RHO > 0 ? q[1][idx] + 2 * psi[1][idx] / ppsi : 0;
+            const double dqz      = DIFF_Z   > 0 ? q[3][idx] + 2 * psi[3][idx] / ppsi : 0;
+            const double dq2rho   = DIFF_RHO > 1 ? q[2][idx] + 2 * psi[2][idx] / ppsi - 2 * SQR(psi[1][idx] / ppsi) : 0;
+            const double dq2z     = DIFF_Z > 1 ? q[6][idx] + 2 * psi[6][idx] / ppsi - 2 * SQR(psi[3][idx] / ppsi) : 0;
+            const double dq2rho_z = DIFF_RHO && DIFF_Z ? q[4][idx] + 2 * psi[4][idx] / ppsi - 2 * psi[1][idx] * psi[3][idx] / psi2 : 0;
+
+            double val;
+
+            DO_EVAL;
+
+            *dst++ = val;
+        }
+    }
+}
+
+#undef FUNC3
+#undef FUNC2
+#undef FUNC