eval: refactor evaluating the metric

1 files changed, 161 insertions, 0 deletions

diff --git a/eval_metric.c b/eval_metric.c
new file mode 100644
index 0000000..f499e32
--- /dev/null
+++ b/eval_metric.c
@@ -0,0 +1,161 @@
+/*
+ * 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/>.
+ */
+
+#include <math.h>
+#include <stddef.h>
+#include <string.h>
+
+#include "brill_data.h"
+#include "internal.h"
+
+#define GRID(arr) (arr[i * nb_coords_rho + j])
+
+
+/* φφ components */
+
+#define COMP phiphi
+
+#define DIFF 0
+#define DO_EVAL val = SQR(rrho) * psi4
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#define DIFF 3
+#define DO_EVAL val = 4 * SQR(rrho) * psi3 * GRID(psi[3])
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#define DIFF 6
+#define DO_EVAL val = 4 * SQR(rrho) * (GRID(psi[6]) * psi3 + 3 * psi2 * SQR(GRID(psi[3])))
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#define DIFF 1
+#define DO_EVAL val = 2 * rrho * psi4 + 4 * SQR(rrho) * psi3 * GRID(psi[1])
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#define DIFF 4
+#define DO_EVAL val =  8 * rrho * psi3 * GRID(psi[3]) +\
+                      12 * SQR(rrho) * psi2 * GRID(psi[3]) * GRID(psi[1]) +\
+                       4 * SQR(rrho) * psi3 * GRID(psi[4])
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#define DIFF 2
+#define DO_EVAL val =  4 * SQR(rrho) * psi3 * GRID(psi[2]) +\
+                      12 * SQR(rrho * ppsi * GRID(psi[1])) +\
+                      16 * rrho * psi3 * GRID(psi[1]) + 2 * psi4
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#undef COMP
+
+/* ρρ components */
+#define COMP rhorho
+
+#define DIFF 0
+#define DO_EVAL val = base
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#define DIFF 3
+#define DO_EVAL val = 2 * base * dqz
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#define DIFF 6
+#define DO_EVAL val = 4 * SQR(dqz) * base + 2 * base * dq2z
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#define DIFF 1
+#define DO_EVAL val = 2 * base * dqrho
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#define DIFF 4
+#define DO_EVAL val = 4 * base * dqrho * dqz + 2 * base * dq2rho_z
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#define DIFF 2
+#define DO_EVAL val = 4 * SQR(dqrho) * base + 2 * base * dq2rho
+#include "eval_metric_template.c"
+#undef DIFF
+#undef DO_EVAL
+
+#undef COMP
+
+typedef void (*EvalMetricFunc)(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]);
+
+static const EvalMetricFunc eval_metric_phiphi[9] = {
+    [0] = eval_metric_phiphi_0,
+    [1] = eval_metric_phiphi_1,
+    [2] = eval_metric_phiphi_2,
+    [3] = eval_metric_phiphi_3,
+    [4] = eval_metric_phiphi_4,
+    [6] = eval_metric_phiphi_6,
+};
+
+static const EvalMetricFunc eval_metric_rhorho[9] = {
+    [0] = eval_metric_rhorho_0,
+    [1] = eval_metric_rhorho_1,
+    [2] = eval_metric_rhorho_2,
+    [3] = eval_metric_rhorho_3,
+    [4] = eval_metric_rhorho_4,
+    [6] = eval_metric_rhorho_6,
+};
+
+void bdi_eval_metric(const BDContext *bd,
+                     const double *rho, int nb_coords_rho,
+                     const double *z,   int nb_coords_z,
+                     enum BDMetricComponent comp,
+                     const unsigned int diff_order[2],
+                     double *psi[9], double *q[9],
+                     double *out, ptrdiff_t out_stride)
+{
+    int diff_idx = diff_order[0] * 3 + diff_order[1];
+    const EvalMetricFunc *functab;
+
+    if (comp != BD_METRIC_COMPONENT_RHORHO &&
+        comp != BD_METRIC_COMPONENT_ZZ     &&
+        comp != BD_METRIC_COMPONENT_PHIPHI) {
+        memset(out, 0, out_stride * nb_coords_z * sizeof(*out));
+        return;
+    }
+
+    functab = (comp == BD_METRIC_COMPONENT_PHIPHI) ?
+              eval_metric_phiphi : eval_metric_rhorho;
+
+    functab[diff_idx](bd, out, out_stride, nb_coords_rho, rho,
+                      nb_coords_z, z, psi, q);
+}