From 2926c4ea6f9e34b8962a08dff372252537473bda Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sun, 28 Aug 2016 10:07:20 +0200 Subject: basis: disentangle the API from the global solver API --- basis.c | 87 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++----- basis.h | 40 +++++++++++++++++++++++++++++ eval.c | 19 +++++++------- init.c | 10 ++++++-- internal.h | 18 ++----------- solve.c | 13 +++++----- 6 files changed, 145 insertions(+), 42 deletions(-) create mode 100644 basis.h diff --git a/basis.c b/basis.c index 0735752..face638 100644 --- a/basis.c +++ b/basis.c @@ -21,20 +21,47 @@ */ #include +#include +#include "basis.h" #include "internal.h" -static double sb_even_eval(double coord, int idx, double sf) + +/* a set of basis functions */ +typedef struct BasisSet { + /* evaluate the idx-th basis function at the specified point*/ + double (*eval) (const BasisSetContext *s, double coord, int idx); + /* evaluate the first derivative of the idx-th basis function at the specified point*/ + double (*eval_diff1)(const BasisSetContext *s, double coord, int idx); + /* evaluate the second derivative of the idx-th basis function at the specified point*/ + double (*eval_diff2)(const BasisSetContext *s, double coord, int idx); + + void (*eval_multi)(const BasisSetContext *s, double coord, + double **out, int order_start, int order_end); + /** + * Get the idx-th collocation point for the specified order. + * idx runs from 0 to order - 1 (inclusive) + */ + double (*colloc_point)(const BasisSetContext *s, int order, int idx); +} BasisSet; + +struct BasisSetContext { + const struct BasisSet *bs; + double sf; +}; + +static double sb_even_eval(const BasisSetContext *s, double coord, int idx) { - double val = atan2(sf, coord); + double val = atan2(s->sf, coord); idx *= 2; // even only return sin((idx + 1) * val); } -static double sb_even_eval_diff1(double coord, int idx, double sf) +static double sb_even_eval_diff1(const BasisSetContext *s, double coord, int idx) { + const double sf = s->sf; double val = atan2(sf, coord); idx *= 2; // even only @@ -42,8 +69,9 @@ static double sb_even_eval_diff1(double coord, int idx, double sf) return - sf * (idx + 1) * cos((idx + 1) * val) / (SQR(sf) + SQR(coord)); } -static double sb_even_eval_diff2(double coord, int idx, double sf) +static double sb_even_eval_diff2(const BasisSetContext *s, double coord, int idx) { + const double sf = s->sf; double val = atan2(sf, coord); idx *= 2; // even only @@ -51,14 +79,14 @@ static double sb_even_eval_diff2(double coord, int idx, double sf) return sf * (idx + 1) * (2 * coord * cos((idx + 1) * val) - sf * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(sf) + SQR(coord)); } -static double sb_even_colloc_point(int order, int idx, double sf) +static double sb_even_colloc_point(const BasisSetContext *s, int order, int idx) { double t; idx = order - idx - 1; t = (idx + 2) * M_PI / (2 * order + 2); - return sf / tan(t); + return s->sf / tan(t); } /* @@ -66,9 +94,54 @@ static double sb_even_colloc_point(int order, int idx, double sf) * SB(x, n) = sin((n + 1) arccot(|x| / L)) * They are symmetric wrt origin and decay as odd powers of x in infinity. */ -const BasisSet bdi_sb_even_basis = { +static const BasisSet sb_even_basis = { .eval = sb_even_eval, .eval_diff1 = sb_even_eval_diff1, .eval_diff2 = sb_even_eval_diff2, .colloc_point = sb_even_colloc_point, }; + +double bdi_basis_eval(BasisSetContext *s, enum BSEvalType type, double coord, int order) +{ + double (*eval)(const BasisSetContext *s, double coord, int order); + + switch (type) { + case BS_EVAL_TYPE_VALUE: eval = s->bs->eval; break; + case BS_EVAL_TYPE_DIFF1: eval = s->bs->eval_diff1; break; + case BS_EVAL_TYPE_DIFF2: eval = s->bs->eval_diff2; break; + } + + return eval(s, coord, order); +} + +double bdi_basis_colloc_point(BasisSetContext *s, int idx, int order) +{ + return s->bs->colloc_point(s, order, idx); +} + +void bdi_basis_free(BasisSetContext **ps) +{ + BasisSetContext *s = *ps; + free(s); + *ps = NULL; +} + +BasisSetContext *bdi_basis_init(enum BasisFamily family, double sf) +{ + BasisSetContext *s = calloc(1, sizeof(*s)); + + if (!s) + return NULL; + + switch (family) { + case BASIS_FAMILY_SB_EVEN: + s->bs = &sb_even_basis; + s->sf = sf; + break; + default: + free(s); + return NULL; + } + + return s; +} diff --git a/basis.h b/basis.h new file mode 100644 index 0000000..6b6f3c4 --- /dev/null +++ b/basis.h @@ -0,0 +1,40 @@ +/* + * Copyright 2015 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_BASIS_H +#define BRILL_DATA_BASIS_H + +enum BSEvalType { + BS_EVAL_TYPE_VALUE, + BS_EVAL_TYPE_DIFF1, + BS_EVAL_TYPE_DIFF2, +}; + +enum BasisFamily { + BASIS_FAMILY_SB_EVEN, +}; + +typedef struct BasisSetContext BasisSetContext; + +BasisSetContext *bdi_basis_init(enum BasisFamily family, double sf); +void bdi_basis_free(BasisSetContext **s); + +double bdi_basis_eval(BasisSetContext *s, enum BSEvalType type, + double coord, int order); +double bdi_basis_colloc_point(BasisSetContext *s, int idx, int order); + +#endif /* BRILL_DATA_BASIS_H */ diff --git a/eval.c b/eval.c index 98fd0a7..c8b4b95 100644 --- a/eval.c +++ b/eval.c @@ -36,8 +36,7 @@ int bd_eval_psi(const BDContext *bd, double *psi, unsigned int psi_stride) { BDPriv *s = bd->priv; - const double *sf = bd->basis_scale_factor; - double (*eval)(double coord, int idx, double sf); + enum BSEvalType type; double *basis_val_rho = NULL, *basis_val_z = NULL; @@ -62,26 +61,26 @@ int bd_eval_psi(const BDContext *bd, } switch (diff_order[0]) { - case 0: eval = s->basis[0]->eval; break; - case 1: eval = s->basis[0]->eval_diff1; break; - case 2: eval = s->basis[0]->eval_diff2; break; + case 0: type = BS_EVAL_TYPE_VALUE; break; + case 1: type = BS_EVAL_TYPE_DIFF1; break; + case 2: type = BS_EVAL_TYPE_DIFF2; break; } for (int i = 0; i < nb_coords_rho; i++) { double rrho = rho[i]; for (int j = 0; j < s->nb_coeffs[0]; j++) - basis_val_rho[i * s->nb_coeffs[0] + j] = eval(rrho, j, sf[0]); + basis_val_rho[i * s->nb_coeffs[0] + j] = bdi_basis_eval(s->basis[0], type, rrho, j); } switch (diff_order[1]) { - case 0: eval = s->basis[1]->eval; break; - case 1: eval = s->basis[1]->eval_diff1; break; - case 2: eval = s->basis[1]->eval_diff2; break; + case 0: type = BS_EVAL_TYPE_VALUE; break; + case 1: type = BS_EVAL_TYPE_DIFF1; break; + case 2: type = BS_EVAL_TYPE_DIFF2; break; } for (int i = 0; i < nb_coords_z; i++) { double zz = z[i]; for (int j = 0; j < s->nb_coeffs[1]; j++) - basis_val_z[i * s->nb_coeffs[1] + j] = eval(zz, j, sf[1]); + basis_val_z[i * s->nb_coeffs[1] + j] = bdi_basis_eval(s->basis[1], type, zz, j); } for (int i = 0; i < nb_coords_z; i++) { diff --git a/init.c b/init.c index 5f9e3e6..fa15325 100644 --- a/init.c +++ b/init.c @@ -22,6 +22,7 @@ #include #include "brill_data.h" +#include "basis.h" #include "internal.h" #include "qfunc.h" @@ -35,8 +36,10 @@ static int brill_init_check_options(BDContext *bd) if (!s->qfunc) return ret; - s->basis[0] = &bdi_sb_even_basis; - s->basis[1] = &bdi_sb_even_basis; + s->basis[0] = bdi_basis_init(BASIS_FAMILY_SB_EVEN, bd->basis_scale_factor[0]); + s->basis[1] = bdi_basis_init(BASIS_FAMILY_SB_EVEN, bd->basis_scale_factor[0]); + if (!s->basis[0] || !s->basis[1]) + return -ENOMEM; s->nb_colloc_points[0] = bd->nb_coeffs[0]; s->nb_colloc_points[1] = bd->nb_coeffs[1]; @@ -108,6 +111,9 @@ void bd_context_free(BDContext **pbd) s = bd->priv; + bdi_basis_free(&s->basis[0]); + bdi_basis_free(&s->basis[1]); + bdi_qfunc_free(&s->qfunc); free(s->coeffs); diff --git a/internal.h b/internal.h index ae9932a..78b3741 100644 --- a/internal.h +++ b/internal.h @@ -22,6 +22,7 @@ #include "brill_data.h" +#include "basis.h" #include "qfunc.h" #define MAX(x, y) ((x) > (y) ? (x) : (y)) @@ -37,23 +38,9 @@ */ #define EPS 1E-08 -/* a set of basis functions */ -typedef struct BasisSet { - /* evaluate the idx-th basis function at the specified point*/ - double (*eval) (double coord, int idx, double sf); - /* evaluate the first derivative of the idx-th basis function at the specified point*/ - double (*eval_diff1)(double coord, int idx, double sf); - /* evaluate the second derivative of the idx-th basis function at the specified point*/ - double (*eval_diff2)(double coord, int idx, double sf); - /** - * Get the idx-th collocation point for the specified order. - * idx runs from 0 to order - 1 (inclusive) - */ - double (*colloc_point)(int order, int idx, double sf); -} BasisSet; typedef struct BDPriv { - const BasisSet *basis[2]; + BasisSetContext *basis[2]; QFuncContext *qfunc; int nb_colloc_points[2]; @@ -66,7 +53,6 @@ typedef struct BDPriv { #define NB_COEFFS(s) (s->nb_coeffs[0] * s->nb_coeffs[1]) #define NB_COLLOC_POINTS(s) (s->nb_colloc_points[0] * s->nb_colloc_points[1]) -extern const BasisSet bdi_sb_even_basis; int bdi_solve(BDContext *bd); diff --git a/solve.c b/solve.c index 5803d87..fe47832 100644 --- a/solve.c +++ b/solve.c @@ -38,7 +38,6 @@ static int brill_construct_matrix(const BDContext *bd, double *mat, double *rhs) { BDPriv *s = bd->priv; - const double *sf = bd->basis_scale_factor; double *basis_val, *basis_dval, *basis_d2val; int idx_coeff_rho, idx_coeff_z, idx_grid_rho, idx_grid_z; @@ -56,11 +55,11 @@ static int brill_construct_matrix(const BDContext *bd, double *mat, } } for (int j = 0; j < s->nb_colloc_points[i]; j++) { - double coord = s->basis[i]->colloc_point(s->nb_coeffs[i], j, sf[i]); + double coord = bdi_basis_colloc_point(s->basis[i], j, s->nb_coeffs[i]); for (int k = 0; k < s->nb_coeffs[i]; k++) { - basis[i][0][j * s->nb_coeffs[i] + k] = s->basis[i]->eval (coord, k, sf[i]); - basis[i][1][j * s->nb_coeffs[i] + k] = s->basis[i]->eval_diff1(coord, k, sf[i]); - basis[i][2][j * s->nb_coeffs[i] + k] = s->basis[i]->eval_diff2(coord, k, sf[i]); + basis[i][0][j * s->nb_coeffs[i] + k] = bdi_basis_eval(s->basis[i], BS_EVAL_TYPE_VALUE, coord, k); + basis[i][1][j * s->nb_coeffs[i] + k] = bdi_basis_eval(s->basis[i], BS_EVAL_TYPE_DIFF1, coord, k); + basis[i][2][j * s->nb_coeffs[i] + k] = bdi_basis_eval(s->basis[i], BS_EVAL_TYPE_DIFF2, coord, k); } } } @@ -73,10 +72,10 @@ static int brill_construct_matrix(const BDContext *bd, double *mat, #define D2BASIS_Z (basis[1][2][idx_grid_z * s->nb_coeffs[1] + idx_coeff_z]) for (idx_grid_z = 0; idx_grid_z < s->nb_colloc_points[1]; idx_grid_z++) { - double z_val = s->basis[1]->colloc_point(s->nb_coeffs[1], idx_grid_z, sf[1]); + double z_val = bdi_basis_colloc_point(s->basis[1], idx_grid_z, s->nb_coeffs[1]); 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 x_val = bdi_basis_colloc_point(s->basis[0], idx_grid_rho, s->nb_coeffs[0]); 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; -- cgit v1.2.3