From 63b25eca76d8481a0fdad339464d1eb617826583 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sun, 28 Aug 2016 10:07:20 +0200 Subject: qfunc: disentangle the API from the global solver API --- eval.c | 36 +++++---------- init.c | 24 ++++------ internal.h | 15 ++----- qfunc.c | 146 +++++++++++++++++++++++++++++++++++++++++++------------------ qfunc.h | 31 +++++++++++++ solve.c | 4 +- 6 files changed, 159 insertions(+), 97 deletions(-) create mode 100644 qfunc.h diff --git a/eval.c b/eval.c index 6c70790..98fd0a7 100644 --- a/eval.c +++ b/eval.c @@ -1,5 +1,5 @@ /* - * Copyright 2014-2015 Anton Khirnov + * Copyright 2014-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 @@ -27,6 +27,7 @@ #include "brill_data.h" #include "internal.h" +#include "qfunc.h" int bd_eval_psi(const BDContext *bd, const double *rho, int nb_coords_rho, @@ -172,14 +173,14 @@ do { \ * derivatives as needed */ #define DO_EVAL_Q_WRAPPER(DO_EVAL, eval_drho, eval_dz, eval_d2rho, eval_d2z, eval_d2rhoz) \ do { \ - const double base = psi4 * exp(2 * s->q_func->q(bd, rrho, zz)); \ + const double base = psi4 * exp(2 * bdi_qfunc_eval(s->qfunc, rrho, zz, 0)); \ double drho, dz, d2rho, d2z, d2rho_z; \ \ - if (eval_drho) drho = s->q_func->dq_rho(bd, rrho, zz) + 2 * GRID(dpsi_rho) / ppsi; \ - if (eval_dz) dz = s->q_func->dq_z(bd, rrho, zz) + 2 * GRID(dpsi_z) / ppsi; \ - if (eval_d2rho) d2rho = s->q_func->d2q_rho(bd, rrho, zz) + 2 * GRID(d2psi_rho) / ppsi - 2 * SQR(GRID(dpsi_rho) / ppsi); \ - if (eval_d2z) d2z = s->q_func->d2q_z(bd, rrho, zz) + 2 * GRID(d2psi_z) / ppsi - 2 * SQR(GRID(dpsi_z) / ppsi); \ - if (eval_d2rhoz) d2rho_z = s->q_func->d2q_rho_z(bd, rrho, zz) + 2 * GRID(d2psi_rho_z) / ppsi - 2 * GRID(dpsi_rho) * GRID(dpsi_z) / psi2; \ + if (eval_drho) drho = bdi_qfunc_eval(s->qfunc, rrho, zz, 1) + 2 * GRID(dpsi_rho) / ppsi; \ + if (eval_dz) dz = bdi_qfunc_eval(s->qfunc, rrho, zz, 3) + 2 * GRID(dpsi_z) / ppsi; \ + if (eval_d2rho) d2rho = bdi_qfunc_eval(s->qfunc, rrho, zz, 2) + 2 * GRID(d2psi_rho) / ppsi - 2 * SQR(GRID(dpsi_rho) / ppsi); \ + if (eval_d2z) d2z = bdi_qfunc_eval(s->qfunc, rrho, zz, 6) + 2 * GRID(d2psi_z) / ppsi - 2 * SQR(GRID(dpsi_z) / ppsi); \ + if (eval_d2rhoz) d2rho_z = bdi_qfunc_eval(s->qfunc, rrho, zz, 4) + 2 * GRID(d2psi_rho_z) / ppsi - 2 * GRID(dpsi_rho) * GRID(dpsi_z) / psi2; \ \ DO_EVAL; \ } while (0) @@ -286,34 +287,19 @@ int bd_eval_q(const BDContext *bd, const double *rho, int nb_coords_rho, double *out, unsigned int out_stride) { BDPriv *s = bd->priv; - double (*eval)(const struct BDContext *bd, double rho, double z); + int diff_idx; if (diff_order[0] > 2 || diff_order[1] > 2) { bdi_log(bd, 0, "Derivatives of higher order than 2 are not supported.\n"); return -ENOSYS; } - switch (diff_order[0]) { - case 0: - switch (diff_order[1]) { - case 0: eval = s->q_func->q; break; - case 1: eval = s->q_func->dq_z; break; - case 2: eval = s->q_func->d2q_z; break; - } - break; - case 1: - switch (diff_order[1]) { - case 0: eval = s->q_func->dq_rho; break; - case 1: eval = s->q_func->d2q_rho_z; break; - } - break; - case 2: eval = s->q_func->d2q_rho; break; - } + diff_idx = diff_order[1] * 3 + diff_order[0]; for (int i = 0; i < nb_coords_z; i++) { double *dst = out + i * out_stride; for (int j = 0; j < nb_coords_rho; j++) - *dst++ = eval(bd, rho[j], z[i]); + *dst++ = bdi_qfunc_eval(s->qfunc, rho[j], z[i], diff_idx); } return 0; diff --git a/init.c b/init.c index 292c555..5f9e3e6 100644 --- a/init.c +++ b/init.c @@ -23,27 +23,17 @@ #include "brill_data.h" #include "internal.h" +#include "qfunc.h" static int brill_init_check_options(BDContext *bd) { BDPriv *s = bd->priv; + int ret; - switch (bd->q_func_type) { - case BD_Q_FUNC_GUNDLACH: - s->q_func = &bdi_q_func_gundlach; - break; - case BD_Q_FUNC_EPPLEY: - if (bd->eppley_n < 4) { - bdi_log(bd, 0, "Invalid n: %d < 4\n", bd->eppley_n); - return -EINVAL; - } - - s->q_func = &bdi_q_func_eppley; - break; - default: - bdi_log(bd, 0, "Unknown q function type: %d\n", bd->q_func_type); - return -EINVAL; - } + ret = bdi_qfunc_init(bd, &s->qfunc, bd->q_func_type, + bd->amplitude, bd->eppley_n); + if (!s->qfunc) + return ret; s->basis[0] = &bdi_sb_even_basis; s->basis[1] = &bdi_sb_even_basis; @@ -118,6 +108,8 @@ void bd_context_free(BDContext **pbd) s = bd->priv; + bdi_qfunc_free(&s->qfunc); + free(s->coeffs); free(bd->priv); diff --git a/internal.h b/internal.h index 19e6d7f..ae9932a 100644 --- a/internal.h +++ b/internal.h @@ -22,6 +22,8 @@ #include "brill_data.h" +#include "qfunc.h" + #define MAX(x, y) ((x) > (y) ? (x) : (y)) #define SQR(x) ((x) * (x)) @@ -50,18 +52,9 @@ typedef struct BasisSet { double (*colloc_point)(int order, int idx, double sf); } BasisSet; -typedef struct QFunc { - double (*q)(const BDContext *bd, double rho, double z); - double (*dq_rho)(const BDContext *bd, double rho, double z); - double (*dq_z)(const BDContext *bd, double rho, double z); - double (*d2q_rho)(const BDContext *bd, double rho, double z); - double (*d2q_z)(const BDContext *bd, double rho, double z); - double (*d2q_rho_z)(const BDContext *bd, double rho, double z); -} QFunc; - typedef struct BDPriv { const BasisSet *basis[2]; - const QFunc *q_func; + QFuncContext *qfunc; int nb_colloc_points[2]; int nb_coeffs[2]; @@ -75,8 +68,6 @@ typedef struct BDPriv { extern const BasisSet bdi_sb_even_basis; -extern const QFunc bdi_q_func_gundlach; -extern const QFunc bdi_q_func_eppley; int bdi_solve(BDContext *bd); diff --git a/qfunc.c b/qfunc.c index e7827df..82557ba 100644 --- a/qfunc.c +++ b/qfunc.c @@ -20,63 +20,74 @@ * definitions of the q functions in the exponential in the Brill data */ +#include #include +#include #include "brill_data.h" #include "internal.h" -static double q_gundlach(const BDContext *bd, double rho, double z) +typedef double (*QFunc)(const struct QFuncContext *s, double rho, double z); + +struct QFuncContext { + const QFunc *qfunc; + + double amplitude; + int n; +}; + +static double q_gundlach(const QFuncContext *s, double rho, double z) { - return bd->amplitude * SQR(rho) * exp(- (SQR(rho) + SQR(z))); + return s->amplitude * SQR(rho) * exp(-(SQR(rho) + SQR(z))); } -static double dq_rho_gundlach(const BDContext *bd, double rho, double z) +static double dq_rho_gundlach(const QFuncContext *s, double rho, double z) { - return bd->amplitude * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho)); + return s->amplitude * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho)); } -static double d2q_rho_gundlach(const BDContext *bd, double rho, double z) +static double d2q_rho_gundlach(const QFuncContext *s, double rho, double z) { double rho2 = SQR(rho); - return bd->amplitude * 2 * exp(-rho2 - SQR(z)) * ((1 - rho2) * (1 - 2 * rho2) - 2 * rho2); + return s->amplitude * 2 * exp(-rho2 - SQR(z)) * ((1 - rho2) * (1 - 2 * rho2) - 2 * rho2); } -static double d2q_rho_z_gundlach(const BDContext *bd, double rho, double z) +static double d2q_rho_z_gundlach(const QFuncContext *s, double rho, double z) { - return -bd->amplitude * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho)); + return -s->amplitude * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho)); } -static double dq_z_gundlach(const BDContext *bd, double rho, double z) +static double dq_z_gundlach(const QFuncContext *s, double rho, double z) { - return - bd->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z)); + return -s->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z)); } -static double d2q_z_gundlach(const BDContext *bd, double rho, double z) +static double d2q_z_gundlach(const QFuncContext *s, double rho, double z) { - return bd->amplitude * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1); + return s->amplitude * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1); } // q function from from PHYSICAL REVIEW D 88, 103009 (2013) // with σ and ρ_0 hardcoded to 1 for now -const QFunc bdi_q_func_gundlach = { - .q = q_gundlach, - .dq_rho = dq_rho_gundlach, - .dq_z = dq_z_gundlach, - .d2q_rho = d2q_rho_gundlach, - .d2q_z = d2q_z_gundlach, - .d2q_rho_z = d2q_rho_z_gundlach, +static const QFunc q_func_gundlach[9] = { + [0] = q_gundlach, + [1] = dq_rho_gundlach, + [3] = dq_z_gundlach, + [2] = d2q_rho_gundlach, + [6] = d2q_z_gundlach, + [4] = d2q_rho_z_gundlach, }; -static double q_eppley(const BDContext *bd, double rho, double z) +static double q_eppley(const QFuncContext *s, double rho, double z) { double r2 = SQR(rho) + SQR(z); - return bd->amplitude * SQR(rho) / (1.0 + pow(r2, bd->eppley_n / 2.0)); + return s->amplitude * SQR(rho) / (1.0 + pow(r2, s->n / 2.0)); } -static double dq_rho_eppley(const BDContext *bd, double rho, double z) +static double dq_rho_eppley(const QFuncContext *s, double rho, double z) { - double A = bd->amplitude; - double n = bd->eppley_n; + double A = s->amplitude; + double n = s->n; double rho2 = SQR(rho); double z2 = SQR(z); @@ -87,10 +98,10 @@ static double dq_rho_eppley(const BDContext *bd, double rho, double z) return A * rho * (2 * (1 + rnm2 * r2) - n * rho2 * rnm2) / SQR(1 + rnm2 * r2); } -static double dq_z_eppley(const BDContext *bd, double rho, double z) +static double dq_z_eppley(const QFuncContext *s, double rho, double z) { - double A = bd->amplitude; - double n = bd->eppley_n; + double A = s->amplitude; + double n = s->n; double rho2 = SQR(rho); double z2 = SQR(z); @@ -101,10 +112,10 @@ static double dq_z_eppley(const BDContext *bd, double rho, double z) return - A * n * rho2 * z * rnm2 / SQR(1 + rnm2 * r2); } -static double d2q_rho_eppley(const BDContext *bd, double rho, double z) +static double d2q_rho_eppley(const QFuncContext *s, double rho, double z) { - double A = bd->amplitude; - double n = bd->eppley_n; + double A = s->amplitude; + double n = s->n; double rho2 = SQR(rho); double z2 = SQR(z); @@ -117,10 +128,10 @@ static double d2q_rho_eppley(const BDContext *bd, double rho, double z) return A * (SQR(n) * SQR(rho2) * rnm4 * (rn - 1) - n * rho2 * rnm4 * (3 * rho2 + 5 * z2) * rnp1 + 2 * SQR(rnp1)) / (rnp1 * SQR(rnp1)); } -static double d2q_z_eppley(const BDContext *bd, double rho, double z) +static double d2q_z_eppley(const QFuncContext *s, double rho, double z) { - double A = bd->amplitude; - double n = bd->eppley_n; + double A = s->amplitude; + double n = s->n; double rho2 = SQR(rho); double z2 = SQR(z); @@ -133,10 +144,10 @@ static double d2q_z_eppley(const BDContext *bd, double rho, double z) return A * n * rho2 * rnm4 * (z2 - rho2 - n * z2 + rn * ((1 + n) * z2 - rho2)) / (rnp1 * SQR(rnp1)); } -static double d2q_rho_z_eppley(const BDContext *bd, double rho, double z) +static double d2q_rho_z_eppley(const QFuncContext *s, double rho, double z) { - double A = bd->amplitude; - double n = bd->eppley_n; + double A = s->amplitude; + double n = s->n; double rho2 = SQR(rho); double z2 = SQR(z); @@ -151,11 +162,60 @@ static double d2q_rho_z_eppley(const BDContext *bd, double rho, double z) // q function from from PHYSICAL REVIEW D 16, 1609 (1977) // with λ hardcoded to 1 for now -const QFunc bdi_q_func_eppley = { - .q = q_eppley, - .dq_rho = dq_rho_eppley, - .dq_z = dq_z_eppley, - .d2q_rho = d2q_rho_eppley, - .d2q_z = d2q_z_eppley, - .d2q_rho_z = d2q_rho_z_eppley, +static const QFunc q_func_eppley[9] = { + [0] = q_eppley, + [1] = dq_rho_eppley, + [3] = dq_z_eppley, + [2] = d2q_rho_eppley, + [6] = d2q_z_eppley, + [4] = d2q_rho_z_eppley, }; + +double bdi_qfunc_eval(QFuncContext *s, double rho, double z, int diff_idx) +{ + return s->qfunc[diff_idx](s, rho, z); +} + +void bdi_qfunc_free(QFuncContext **ps) +{ + QFuncContext *s = *ps; + free(s); + *ps = NULL; +} + +int bdi_qfunc_init(const BDContext *bd, QFuncContext **out, enum BDQFuncType type, + double amplitude, int n) +{ + QFuncContext *s = calloc(1, sizeof(*s)); + int ret = 0; + + if (!s) + return -ENOMEM; + + switch (type) { + case BD_Q_FUNC_GUNDLACH: + s->qfunc = q_func_gundlach; + break; + case BD_Q_FUNC_EPPLEY: + if (n < 4) { + bdi_log(bd, 0, "Invalid n: %d < 4\n", n); + ret = -EINVAL; + goto fail; + } + s->qfunc = q_func_eppley; + break; + default: + bdi_log(bd, 0, "Unknown q function type: %d\n", type); + ret = -EINVAL; + goto fail; + } + + s->amplitude = amplitude; + s->n = n; + + *out = s; + return 0; +fail: + free(s); + return ret; +} diff --git a/qfunc.h b/qfunc.h new file mode 100644 index 0000000..be3f60e --- /dev/null +++ b/qfunc.h @@ -0,0 +1,31 @@ +/* + * 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 . + */ + +#ifndef BRILL_DATA_QFUNC_H +#define BRILL_DATA_QFUNC_H + +#include "brill_data.h" + +typedef struct QFuncContext QFuncContext; + +int bdi_qfunc_init(const BDContext *bd, QFuncContext **out, enum BDQFuncType type, + double amplitude, int n); +void bdi_qfunc_free(QFuncContext **ps); + +double bdi_qfunc_eval(QFuncContext *s, double rho, double z, int diff_idx); + +#endif /* BRILL_DATA_QFUNC_H */ diff --git a/solve.c b/solve.c index 64dcb55..5803d87 100644 --- a/solve.c +++ b/solve.c @@ -32,6 +32,7 @@ #include "brill_data.h" #include "internal.h" +#include "qfunc.h" static int brill_construct_matrix(const BDContext *bd, double *mat, double *rhs) @@ -76,7 +77,8 @@ static int brill_construct_matrix(const BDContext *bd, double *mat, for (idx_grid_rho = 0; idx_grid_rho < s->nb_colloc_points[0]; idx_grid_rho++) { double x_val = s->basis[0]->colloc_point(s->nb_coeffs[0], idx_grid_rho, sf[0]); - double d2q = s->q_func->d2q_rho(bd, x_val, z_val) + s->q_func->d2q_z(bd, x_val, z_val); + double d2q = bdi_qfunc_eval(s->qfunc, x_val, z_val, 2) + + bdi_qfunc_eval(s->qfunc, x_val, z_val, 6); int idx_grid = idx_grid_z * s->nb_colloc_points[0] + idx_grid_rho; for (idx_coeff_z = 0; idx_coeff_z < s->nb_coeffs[1]; idx_coeff_z++) -- cgit v1.2.3