From 4633e6db474e7825979bef26e68628d492a63eb1 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sun, 28 Aug 2016 10:53:54 +0200 Subject: eval: refactor evaluating the metric --- eval_metric.c | 161 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 161 insertions(+) create mode 100644 eval_metric.c (limited to 'eval_metric.c') 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 + * + * 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 . + */ + +#include +#include +#include + +#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); +} -- cgit v1.2.3