From 41f29cbf3076db51b96f240d27d432ef31b8aa71 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Tue, 14 Apr 2015 16:00:24 +0200 Subject: A major rewrite. Split the code into multiple files, drop the GSL dependency, introduce configurable logging, random cleanups. --- Makefile | 13 +- README | 6 +- basis.c | 74 ++++++ brill_data.c | 843 ---------------------------------------------------------- brill_data.h | 89 ++++--- brill_data.py | 37 ++- eval.c | 320 ++++++++++++++++++++++ init.c | 133 +++++++++ internal.h | 83 ++++++ log.c | 41 +++ qfunc.c | 161 +++++++++++ solve.c | 230 ++++++++++++++++ 12 files changed, 1136 insertions(+), 894 deletions(-) create mode 100644 basis.c delete mode 100644 brill_data.c create mode 100644 eval.c create mode 100644 init.c create mode 100644 internal.h create mode 100644 log.c create mode 100644 qfunc.c create mode 100644 solve.c diff --git a/Makefile b/Makefile index aec6c0b..0efa6c0 100644 --- a/Makefile +++ b/Makefile @@ -1,14 +1,19 @@ -CFLAGS = -std=gnu99 -D_XOPEN_SOURCE=700 -fPIC -LDFLAGS = --version-script=libbrilldata.v -shared -lgsl -lm -llapacke +CFLAGS = -std=c99 -D_XOPEN_SOURCE=700 -fPIC -g +LDFLAGS = --version-script=libbrilldata.v -shared -lm -llapacke TARGET = libbrilldata.so -OBJECTS = brill_data.o +OBJECTS = basis.o \ + eval.o \ + init.o \ + log.o \ + qfunc.o \ + solve.o all: $(TARGET) $(TARGET): $(OBJECTS) - ld ${LDFLAGS} -o $@ $< + ld ${LDFLAGS} -o $@ $(OBJECTS) clean: -rm $(OBJECTS) $(TARGET) diff --git a/README b/README index 3125972..9f5492f 100644 --- a/README +++ b/README @@ -13,9 +13,9 @@ matrix is inverted with LAPACK. Building and installation ========================= -The library requires GSL and LAPACKE (the C interface to LAPACK) to be present -where the compiler and linker can find them. A C99-compliant compiler and a -POSIX environemtn are expected. +The library requires LAPACKE (the C interface to LAPACK) to be present where the +compiler and linker can find it. A C99-compliant compiler and a POSIX +environment are expected. Simply running 'make' will then build the shared library libbrilldata.so. That must be copied manually to where your linker will find it (or set diff --git a/basis.c b/basis.c new file mode 100644 index 0000000..0735752 --- /dev/null +++ b/basis.c @@ -0,0 +1,74 @@ +/* + * Copyright 2014-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 . + */ + +/** + * @file + * definitions of the basis functions for the spectral method + */ + +#include + +#include "internal.h" + +static double sb_even_eval(double coord, int idx, double sf) +{ + double val = atan2(sf, coord); + + idx *= 2; // even only + + return sin((idx + 1) * val); +} + +static double sb_even_eval_diff1(double coord, int idx, double sf) +{ + double val = atan2(sf, coord); + + idx *= 2; // even only + + return - sf * (idx + 1) * cos((idx + 1) * val) / (SQR(sf) + SQR(coord)); +} + +static double sb_even_eval_diff2(double coord, int idx, double sf) +{ + double val = atan2(sf, coord); + + idx *= 2; // even only + + 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) +{ + double t; + + idx = order - idx - 1; + + t = (idx + 2) * M_PI / (2 * order + 2); + return sf / tan(t); +} + +/* + * The basis of even (n = 2 * idx) SB functions (Boyd 2000, Ch 17.9) + * 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 = { + .eval = sb_even_eval, + .eval_diff1 = sb_even_eval_diff1, + .eval_diff2 = sb_even_eval_diff2, + .colloc_point = sb_even_colloc_point, +}; diff --git a/brill_data.c b/brill_data.c deleted file mode 100644 index 9c98ec4..0000000 --- a/brill_data.c +++ /dev/null @@ -1,843 +0,0 @@ -/* - * Copyright 2014 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 -#include -#include - -#include -#include - -#include - -#include "brill_data.h" - -#define ACC_TEST 1 - -#define MAX(x, y) ((x) > (y) ? (x) : (y)) -#define MIN(x, y) ((x) > (y) ? (y) : (x)) -#define SQR(x) ((x) * (x)) -#define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0) - -/* - * small number to avoid r=0 singularities - */ -#define EPS 1E-08 - -/* a set of basis functions */ -typedef struct BasisSet { - /* evaluate the idx-th basis function at the specified point*/ - long double (*eval) (long double coord, int idx, long double sf); - /* evaluate the first derivative of the idx-th basis function at the specified point*/ - long double (*eval_diff1)(long double coord, int idx, long double sf); - /* evaluate the second derivative of the idx-th basis function at the specified point*/ - long double (*eval_diff2)(long double coord, int idx, long double sf); - /** - * Get the idx-th collocation point for the specified order. - * idx runs from 0 to order - 1 (inclusive) - */ - long double (*colloc_point)(int order, int idx, long double sf); -} BasisSet; - -typedef struct QFunc { - double (*q)(struct BDContext *bd, double rho, double z); - double (*dq_rho)(struct BDContext *bd, double rho, double z); - double (*dq_z)(struct BDContext *bd, double rho, double z); - double (*d2q_rho)(struct BDContext *bd, double rho, double z); - double (*d2q_z)(struct BDContext *bd, double rho, double z); - double (*d2q_rho_z)(struct BDContext *bd, double rho, double z); -} QFunc; - -typedef struct BDPriv { - const BasisSet *basis; - const QFunc *q_func; - - int nb_colloc_points; - int nb_colloc_points_rho; - int nb_colloc_points_z; - - int colloc_grid_order_rho; - int colloc_grid_order_z; - int nb_coeffs; - - gsl_vector *coeffs; -} BDPriv; - -/* - * The basis of even (n = 2 * idx) SB functions (Boyd 2000, Ch 17.9) - * SB(x, n) = sin((n + 1) arccot(|x| / L)) - * They are symmetric wrt origin and decay as 1/x in infinity. - */ - -static long double sb_even_eval(long double coord, int idx, long double sf) -{ - long double val = atan2l(sf, coord); - - idx *= 2; // even only - - return sinl((idx + 1) * val); -} - -static long double sb_even_eval_diff1(long double coord, int idx, long double sf) -{ - long double val = atan2l(sf, coord); - - idx *= 2; // even only - - return - sf * (idx + 1) * SGN(coord) * cosl((idx + 1) * val) / (SQR(sf) + SQR(coord)); -} - -static long double sb_even_eval_diff2(long double coord, int idx, long double sf) -{ - long double val = atan2l(sf, coord); - - idx *= 2; // even only - - return sf * (idx + 1) * SGN(coord) * (2 * fabsl(coord) * cosl((idx + 1) * val) - sf * (idx + 1) * sinl((idx + 1) * val)) / SQR(SQR(sf) + SQR(coord)); -} - -static long double sb_even_colloc_point(int order, int idx, long double sf) -{ - long double t; - - idx = order - idx - 1; - - t = (idx + 2) * M_PI / (2 *(order + 2)); - return sf / tanl(t); -} - -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, -}; - -static int brill_construct_matrix(BDContext *bd, gsl_matrix *mat, - gsl_vector *rhs) -{ - BDPriv *s = bd->priv; - const long double sfr = bd->basis_scale_factor_rho, sfz = bd->basis_scale_factor_z; - - double *basis_val, *basis_dval, *basis_d2val; - int idx_coeff_rho, idx_coeff_z, idx_grid_rho, idx_grid_z; - - /* precompute the basis values we will need */ - // FIXME: assumes same number of points in ρ and z directions - basis_val = malloc(sizeof(*basis_val) * s->nb_colloc_points_rho * bd->nb_coeffs_rho); - basis_dval = malloc(sizeof(*basis_dval) * s->nb_colloc_points_rho * bd->nb_coeffs_rho); - basis_d2val = malloc(sizeof(*basis_d2val) * s->nb_colloc_points_rho * bd->nb_coeffs_rho); - for (int i = 0; i < s->nb_colloc_points_rho; i++) { - double coord = s->basis->colloc_point(s->colloc_grid_order_rho, i, sfr); - for (int j = 0; j < bd->nb_coeffs_rho; j++) { - basis_val [i * bd->nb_coeffs_rho + j] = s->basis->eval(coord, j, sfr); - basis_dval [i * bd->nb_coeffs_rho + j] = s->basis->eval_diff1(coord, j, sfr); - basis_d2val[i * bd->nb_coeffs_rho + j] = s->basis->eval_diff2(coord, j, sfr); - } - } - - for (idx_grid_z = 0; idx_grid_z < s->nb_colloc_points_z; idx_grid_z++) { - double z_physical = s->basis->colloc_point(s->colloc_grid_order_z, idx_grid_z, sfz); - - for (idx_grid_rho = 0; idx_grid_rho < s->nb_colloc_points_rho; idx_grid_rho++) { - double x_physical = s->basis->colloc_point(s->colloc_grid_order_rho, idx_grid_rho, sfr); - double d2q = s->q_func->d2q_rho(bd, x_physical, z_physical) + s->q_func->d2q_z(bd, x_physical, z_physical); - int idx_grid = idx_grid_z * s->nb_colloc_points_z + idx_grid_rho; - - for (idx_coeff_z = 0; idx_coeff_z < bd->nb_coeffs_z; idx_coeff_z++) - for (idx_coeff_rho = 0; idx_coeff_rho < bd->nb_coeffs_rho; idx_coeff_rho++) { - int idx_coeff = idx_coeff_z * bd->nb_coeffs_z + idx_coeff_rho; - - mat->data[idx_grid + mat->tda * idx_coeff] = basis_d2val[idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] + - basis_d2val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] * basis_val[idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] + - basis_val [idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] * 0.25 * d2q; - if (x_physical > EPS) - mat->data[idx_grid + mat->tda * idx_coeff] += basis_dval [idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] / x_physical; - else - mat->data[idx_grid + mat->tda * idx_coeff] += basis_d2val [idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z]; - } - rhs->data[idx_grid] = -0.25 * d2q; - } - } - - free(basis_val); - free(basis_dval); - free(basis_d2val); - - return 0; -} - -static int brill_solve_linear(gsl_matrix *mat, gsl_vector **px, gsl_vector **prhs) -{ - int *ipiv; - gsl_matrix *mat_f; - double cond, ferr, berr, rpivot; - - gsl_vector *x = *px; - gsl_vector *rhs = *prhs; - char equed = 'N'; - int ret = 0; - - ipiv = malloc(mat->size1 * sizeof(*ipiv)); - mat_f = gsl_matrix_alloc(mat->size1, mat->size2); - if (!ipiv || !mat_f) { - ret = -ENOMEM; - goto fail; - } - - LAPACKE_dgesvx(LAPACK_COL_MAJOR, 'N', 'N', mat->size1, 1, - mat->data, mat->tda, mat_f->data, mat_f->tda, ipiv, &equed, - NULL, NULL, rhs->data, rhs->size, x->data, x->size, - &cond, &ferr, &berr, &rpivot); - - fprintf(stderr, "LU factorization solution to a %zdx%zd matrix: " - "condition number %16.16g; forward error %16.16g backward error %16.16g\n", - mat->size1, mat->size2, cond, ferr, berr); -fail: - gsl_matrix_free(mat_f); - free(ipiv); - - return ret; -} - -static int brill_solve_svd(gsl_matrix *mat, gsl_vector **px, gsl_vector **prhs) -{ - double *sv; - int rank; - - gsl_vector *x = *px; - gsl_vector *rhs = *prhs; - - sv = malloc(sizeof(*sv) * x->size); - if (!sv) - return -ENOMEM; - - LAPACKE_dgelsd(LAPACK_COL_MAJOR, rhs->size, x->size, 1, mat->data, mat->tda, - rhs->data, rhs->size, sv, -1, &rank); - - fprintf(stderr, "Least squares SVD solution to a %zdx%zd matrix: " - "rank %d, condition number %16.16g\n", mat->size1, mat->size2, - rank, sv[x->size - 1] / sv[0]); - - gsl_vector_free(x); - - *px = rhs; - (*px)->size = mat->size1; - *prhs = NULL; - - free(sv); - - return 0; -} - -/* - * Solve the equation - * Δψ + ¼ ψ (∂²q/∂r² + ∂²q/∂z²) = 0 - * for the coefficients of spectral approximation of ψ: - * ψ(r, z) = 1 + ΣaᵢⱼTᵢ(r)Tⱼ(z) - * where i = { 0, ... , bd->nb_coeffs_rho }; - * j = { 0, ... , bd->nb_coeffs_z }; - * q(r, z) and Tᵢ(x) are defined by bd->q_func, bd->laplace_q_func and - * bd->basis. - * - * The cofficients are exported in bd->priv->coeffs - */ -static int brill_solve(BDContext *bd) -{ - BDPriv *s = bd->priv; - gsl_matrix *mat = NULL; - gsl_vector *coeffs = NULL, *rhs = NULL; - -#if ACC_TEST - gsl_vector *rhs_bkp = NULL; - gsl_matrix *mat_bkp = NULL; -#endif - - int ret = 0; - - /* allocate and fill the pseudospectral matrix */ - mat = gsl_matrix_alloc (s->nb_coeffs, s->nb_colloc_points); // inverted order for lapack - coeffs = gsl_vector_alloc (s->nb_coeffs); - rhs = gsl_vector_calloc(s->nb_colloc_points); - if (!mat || !coeffs || !rhs) { - ret = -ENOMEM; - goto fail; - } - - /* fill the matrix */ - ret = brill_construct_matrix(bd, mat, rhs); - if (ret < 0) - goto fail; - -#if ACC_TEST - /* make backups of the matrix and RHS, since they might be destroyed later */ - mat_bkp = gsl_matrix_alloc(mat->size2, mat->size1); - rhs_bkp = gsl_vector_alloc(rhs->size); - if (!mat_bkp || !rhs_bkp) { - ret = -ENOMEM; - goto fail; - } - - gsl_vector_memcpy(rhs_bkp, rhs); - gsl_matrix_transpose_memcpy(mat_bkp, mat); -#endif - - /* solve for the coeffs */ - if (s->nb_colloc_points > s->nb_coeffs) - ret = brill_solve_svd(mat, &coeffs, &rhs); - else - ret = brill_solve_linear(mat, &coeffs, &rhs); - - /* export the result */ - s->coeffs = coeffs; - -#if ACC_TEST - { - double min, max, rmin, rmax; - - gsl_vector_minmax(rhs_bkp, &rmin, &rmax); - gsl_blas_dgemv(CblasNoTrans, 1.0, mat_bkp, coeffs, -1.0, rhs_bkp); - gsl_vector_minmax(rhs_bkp, &min, &max); - - fprintf(stderr, "max(|A·x - rhs|) / max(|rhs|): %16.16g\n", - MAX(fabs(min), fabs(max)) / MAX(fabs(rmin), fabs(rmax))); - } -#endif - -fail: - -#if ACC_TEST - gsl_matrix_free(mat_bkp); - gsl_vector_free(rhs_bkp); -#endif - - if (ret < 0) - gsl_vector_free(coeffs); - - gsl_vector_free(rhs); - gsl_matrix_free(mat); - - return ret; -} - -// q function form from PHYSICAL REVIEW D 88, 103009 (2013) -// with σ and ρ_0 hardcoded to 1 for now -static double q_gundlach(BDContext *bd, double rho, double z) -{ - return bd->amplitude * SQR(rho) * exp(- (SQR(rho) + SQR(z))); -} - -static double dq_rho_gundlach(BDContext *bd, double rho, double z) -{ - return bd->amplitude * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho)); -} - -static double d2q_rho_gundlach(BDContext *bd, double rho, double z) -{ - double rho2 = SQR(rho); - return bd->amplitude * 2 * exp(-rho2 - SQR(z)) * ((1 - rho2) * (1 - 2 * rho2) - 2 * rho2); -} - -static double d2q_rho_z_gundlach(BDContext *bd, double rho, double z) -{ - return -bd->amplitude * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho)); -} - -static double dq_z_gundlach(BDContext *bd, double rho, double z) -{ - return - bd->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z)); -} - -static double d2q_z_gundlach(BDContext *bd, double rho, double z) -{ - return bd->amplitude * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1); -} - -static const QFunc 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 double q_eppley(BDContext *bd, double rho, double z) -{ - double r2 = SQR(rho) + SQR(z); - return bd->amplitude * SQR(rho) / (1.0 + pow(r2, bd->eppley_n / 2.0)); -} - -static double dq_rho_eppley(BDContext *bd, double rho, double z) -{ - double A = bd->amplitude; - double n = bd->eppley_n; - double rho2 = SQR(rho); - double z2 = SQR(z); - - double r2 = rho2 + z2; - double r = sqrt(r2); - double rnm2 = pow(r, n - 2); - - return A * rho * (2 * (1 + rnm2 * r2) - n * rho2 * rnm2) / SQR(1 + rnm2 * r2); -} - -static double dq_z_eppley(BDContext *bd, double rho, double z) -{ - double A = bd->amplitude; - double n = bd->eppley_n; - double rho2 = SQR(rho); - double z2 = SQR(z); - - double r2 = rho2 + z2; - double r = sqrt(r2); - double rnm2 = pow(r, n - 2); - - return - A * n * rho2 * z * rnm2 / SQR(1 + rnm2 * r2); -} - -static double d2q_rho_eppley(BDContext *bd, double rho, double z) -{ - double A = bd->amplitude; - double n = bd->eppley_n; - double rho2 = SQR(rho); - double z2 = SQR(z); - - double r2 = rho2 + z2; - double r = sqrt(r2); - double rnm4 = pow(r, n - 4); - double rn = rnm4 * SQR(r2); - double rnp1 = 1 + rn; - - 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(BDContext *bd, double rho, double z) -{ - double A = bd->amplitude; - double n = bd->eppley_n; - double rho2 = SQR(rho); - double z2 = SQR(z); - - double r2 = rho2 + z2; - double r = sqrt(r2); - double rnm4 = pow(r, n - 4); - double rn = rnm4 * SQR(r2); - double rnp1 = 1 + rn; - - return A * n * rho2 * rnm4 * (z2 - rho2 - n * z2 + rn * ((1 + n) * z2 - rho2)) / (rnp1 * SQR(rnp1)); -} - -static double d2q_rho_z_eppley(BDContext *bd, double rho, double z) -{ - double A = bd->amplitude; - double n = bd->eppley_n; - double rho2 = SQR(rho); - double z2 = SQR(z); - - double r2 = rho2 + z2; - double r = sqrt(r2); - double rnm4 = pow(r, n - 4); - double rn = rnm4 * SQR(r2); - double rnp1 = 1 + rn; - - return A * n * rho * z * rnm4 * (rn * (n * rho2 - 2 * z2) - 2 * z2 - n * rho2) / (rnp1 * SQR(rnp1)); -} - -static const QFunc 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 double d2q_eppley(BDContext *bd, double rho, double z) -{ - double rho2 = SQR(rho); - double z2 = SQR(z); - double rho4 = SQR(rho2); - double z4 = SQR(z2); - double z6 = z4 * z2; - - double r2 = rho2 + z2; - double r = sqrt(r2); - - return bd->amplitude * (2.0 + r2 * (7 * rho4 * rho4 + 23 * rho4 * rho2 * z2 + 27 * rho4 * z4 + rho2 * (13 * z6 - 41 * r) + 2 * z2 * (z6 + 2 * r))) / pow(1.0 + r2 * r2 * r, 3); -} - -static int brill_init_check_options(BDContext *bd) -{ - BDPriv *s = bd->priv; - - switch (bd->q_func_type) { - case BD_Q_FUNC_GUNDLACH: - s->q_func = &q_func_gundlach; - break; - case BD_Q_FUNC_EPPLEY: - if (bd->eppley_n < 4) { - fprintf(stderr, "Invalid n: %d < 4\n", bd->eppley_n); - return -EINVAL; - } - - s->q_func = &q_func_eppley; - break; - default: - fprintf(stderr, "Unknown q function type: %d\n", bd->q_func_type); - return -EINVAL; - } - - if (!bd->nb_coeffs_z) - bd->nb_coeffs_z = bd->nb_coeffs_rho; - - if (bd->nb_coeffs_rho != bd->nb_coeffs_z) { - fprintf(stderr, "Different ρ and z basis orders are not supported.\n"); - return -EINVAL; - } - - s->basis = &sb_even_basis; - - s->nb_coeffs = bd->nb_coeffs_rho * bd->nb_coeffs_z; - - s->nb_colloc_points_rho = bd->nb_coeffs_rho + bd->overdet_rho; - s->nb_colloc_points_z = bd->nb_coeffs_z + bd->overdet_z; - - s->nb_colloc_points = s->nb_colloc_points_rho * s->nb_colloc_points_z; - - s->colloc_grid_order_rho = bd->nb_coeffs_rho + bd->colloc_grid_offset_rho; - s->colloc_grid_order_z = bd->nb_coeffs_z + bd->colloc_grid_offset_z; - - return 0; -} - -int bd_solve(BDContext *bd) -{ - BDPriv *s = bd->priv; - int ret; - - if (bd->psi_minus1_coeffs) { - fprintf(stderr, "Solve called multiple times on the same context.\n"); - return -EINVAL; - } - - ret = brill_init_check_options(bd); - if (ret < 0) - return ret; - - ret = brill_solve(bd); - if (ret < 0) - return ret; - - bd->psi_minus1_coeffs = s->coeffs->data; - bd->stride = bd->nb_coeffs_rho; - - return 0; -} - -BDContext *bd_context_alloc(void) -{ - BDContext *bd = calloc(1, sizeof(*bd)); - - if (!bd) - return NULL; - - bd->q_func_type = BD_Q_FUNC_GUNDLACH; - bd->amplitude = 1.0; - bd->eppley_n = 5; - - bd->nb_coeffs_rho = 80; - bd->nb_coeffs_z = 0; - - bd->colloc_grid_offset_rho = 3; - bd->colloc_grid_offset_z = 3; - - bd->basis_scale_factor_rho = 3.0; - bd->basis_scale_factor_z = 3.0; - - bd->priv = calloc(1, sizeof(BDPriv)); - if (!bd->priv) { - free(bd); - return NULL; - } - - return bd; -} - -void bd_context_free(BDContext **pbd) -{ - BDContext *bd = *pbd; - BDPriv *s; - - if (!bd) - return; - - s = bd->priv; - - gsl_vector_free(s->coeffs); - - free(bd->priv); - free(bd); - *pbd = NULL; -} - -int bd_eval_psi(BDContext *bd, const double *rho, int nb_coords_rho, - const double *z, int nb_coords_z, - const unsigned int diff_order[2], - double *psi) -{ - BDPriv *s = bd->priv; - const long double sfr = bd->basis_scale_factor_rho, sfz = bd->basis_scale_factor_z; - long double (*eval)(long double coord, int idx, long double sf); - - long double *basis_val_rho, *basis_val_z; - - long double add = 0.0; - - if (diff_order[0] > 2 || diff_order[1] > 2) { - fprintf(stderr, "Derivatives of higher order than 2 are not supported.\n"); - return -ENOSYS; - } - - if (diff_order[0] == 0 && diff_order[1] == 0) - add = 1.0; - - /* precompute the basis values on the grid points */ - basis_val_rho = malloc(sizeof(*basis_val_rho) * bd->nb_coeffs_rho * nb_coords_rho); - basis_val_z = malloc(sizeof(*basis_val_z) * bd->nb_coeffs_z * nb_coords_z); - - switch (diff_order[0]) { - case 0: eval = s->basis->eval; break; - case 1: eval = s->basis->eval_diff1; break; - case 2: eval = s->basis->eval_diff2; break; - } - - for (int i = 0; i < nb_coords_rho; i++) { - double rrho = rho[i]; - for (int j = 0; j < bd->nb_coeffs_rho; j++) - basis_val_rho[i * bd->nb_coeffs_rho + j] = eval(rrho, j, sfr); - } - - switch (diff_order[1]) { - case 0: eval = s->basis->eval; break; - case 1: eval = s->basis->eval_diff1; break; - case 2: eval = s->basis->eval_diff2; break; - } - for (int i = 0; i < nb_coords_z; i++) { - double zz = z[i]; - for (int j = 0; j < bd->nb_coeffs_z; j++) - basis_val_z[i * bd->nb_coeffs_z + j] = eval(zz, j, sfz); - } - - for (int i = 0; i < nb_coords_z; i++) - for (int j = 0; j < nb_coords_rho; j++) { - long double ppsi = add; - - for (int m = 0; m < bd->nb_coeffs_z; m++) - for (int n = 0; n < bd->nb_coeffs_rho; n++) - ppsi += s->coeffs->data[m * bd->nb_coeffs_rho + n] * basis_val_rho[j * bd->nb_coeffs_rho + n] * basis_val_z[i * bd->nb_coeffs_z + m]; - - psi[i * nb_coords_rho + j] = ppsi; - } - - free(basis_val_rho); - free(basis_val_z); - - return 0; - -} - -int bd_eval_metric(BDContext *bd, const double *rho, int nb_coords_rho, - const double *z, int nb_coords_z, - const unsigned int comp[2], const unsigned int diff_order[2], - double *out) -{ - BDPriv *s = bd->priv; - const int nb_coords = nb_coords_rho * nb_coords_z; - double *psi = NULL, *dpsi_rho = NULL, *dpsi_z = NULL, *d2psi_rho = NULL, *d2psi_rho_z = NULL, *d2psi_z = NULL; - int ret = 0; - - /* check the parameters for validity */ - if (comp[0] > 2 || comp[1] > 2) { - fprintf(stderr, "Invalid component indices: %d%d\n", comp[0], comp[1]); - return -EINVAL; - } - - if (comp[0] != comp[1]) { - memset(out, 0, nb_coords * sizeof(*out)); - return 0; - } - - if (diff_order[0] + diff_order[1] > 2) { - fprintf(stderr, "At most second order derivatives are supported.\n"); - return -ENOSYS; - } - - /* evaluate the conformal factor and its derivatives if necessary */ -#define EVAL_PSI(arr, order) \ -do { \ - arr = malloc(nb_coords * sizeof(*arr)); \ - if (!arr) { \ - ret = -ENOMEM; \ - goto fail; \ - } \ - \ - ret = bd_eval_psi(bd, rho, nb_coords_rho, z, nb_coords_z, \ - order, arr); \ - if (ret < 0) \ - goto fail; \ -} while (0) - - EVAL_PSI(psi, ((unsigned int[2]){ 0, 0 })); - if (diff_order[0]) - EVAL_PSI(dpsi_rho, ((unsigned int[2]){ 1, 0 })); - if (diff_order[0] > 1) - EVAL_PSI(d2psi_rho, ((unsigned int[2]){ 2, 0 })); - if (diff_order[1]) - EVAL_PSI(dpsi_z, ((unsigned int[2]){ 0, 1 })); - if (diff_order[1] > 1) - EVAL_PSI(d2psi_z, ((unsigned int[2]){ 0, 2 })); - if (diff_order[0] && diff_order[1]) - EVAL_PSI(d2psi_rho_z, ((unsigned int[2]){ 1, 1 })); - -/* the template for the loop over the grid points */ -#define EVAL_LOOP(DO_EVAL) \ -do { \ - for (int i = 0; i < nb_coords_z; i++) { \ - const double zz = z[i]; \ - \ - for (int j = 0; j < nb_coords_rho; j++) { \ - const double rrho = rho[j]; \ - const double ppsi = psi[i * nb_coords_rho + j]; \ - const double psi2 = SQR(ppsi); \ - const double psi3 = psi2 * ppsi; \ - const double psi4 = SQR(psi2); \ - double val; \ - \ - DO_EVAL; \ - \ - out[i * nb_coords_rho + j] = val; \ - } \ - } \ -} while (0) - -#define GRID(arr) (arr[i * nb_coords_rho + j]) - - if (comp[0] == 2) { - // γ_φφ - switch (diff_order[0]) { - case 0: - switch (diff_order[1]) { - case 0: EVAL_LOOP(val = SQR(rrho) * psi4); break; - case 1: EVAL_LOOP(val = 4 * SQR(rrho) * psi3 * GRID(dpsi_z)); break; // ∂_z - case 2: EVAL_LOOP(val = 4 * SQR(rrho) * (GRID(d2psi_z) * psi3 + 3 * psi2 * SQR(GRID(dpsi_z)))); break; // ∂_z ∂_z - } - break; - case 1: - switch (diff_order[1]) { - case 0: EVAL_LOOP(val = 2 * rrho * psi4 + 4 * SQR(rrho) * psi3 * GRID(dpsi_rho)); break; // ∂_ρ - case 1: EVAL_LOOP(val = 8 * rrho * psi3 * GRID(dpsi_z) + 12 * SQR(rrho) * psi2 * GRID(dpsi_z) * GRID(dpsi_rho) + 4 * SQR(rrho) * psi3 * GRID(d2psi_rho_z)); break; // ∂_ρ ∂_z - } - break; - case 2: EVAL_LOOP(val = 4 * SQR(rrho) * psi3 * GRID(d2psi_rho) + 12 * SQR(rrho * ppsi * GRID(dpsi_rho)) + 16 * rrho * psi3 * GRID(dpsi_rho) + 2 * psi4); break; // ∂_ρ ∂_ρ - } - } else { - // γ_ρρ / γ_zz - -/* a wrapper around the actual evaluation expression that provides the q function and its - * 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)); \ - 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; \ - \ - DO_EVAL; \ -} while (0) - - switch (diff_order[0]) { - case 0: - switch (diff_order[1]) { - case 0: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = base, 0, 0, 0, 0, 0)); break; - case 1: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = 2 * base * dz, 0, 1, 0, 0, 0)); break; // ∂_z - case 2: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = 4 * SQR(dz) * base + 2 * base * d2z, 0, 1, 0, 1, 0)); break; // ∂_z ∂_z - } - break; - case 1: - switch (diff_order[1]) { - case 0: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = 2 * base * drho, 1, 0, 0, 0, 0)); break; // ∂_ρ - case 1: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = 4 * base * drho * dz + 2 * base * d2rho_z, 1, 1, 0, 0, 1)); break; // ∂_ρ ∂_z - } - break; - case 2: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = 4 * SQR(drho) * base + 2 * base * d2rho, 1, 0, 1, 0, 0)); break; // ∂_ρ ∂_ρ - } - } - -fail: - free(psi); - free(dpsi_rho); - free(dpsi_z); - free(d2psi_rho); - free(d2psi_rho_z); - free(d2psi_z); - - return ret; -} - -int bd_eval_q(BDContext *bd, const double *rho, int nb_coords_rho, - const double *z, int nb_coords_z, const unsigned int diff_order[2], - double *out) -{ - BDPriv *s = bd->priv; - double (*eval)(struct BDContext *bd, double rho, double z); - - if (diff_order[0] > 2 || diff_order[1] > 2) { - fprintf(stderr, "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; - } - - for (int i = 0; i < nb_coords_z; i++) - for (int j = 0; j < nb_coords_rho; j++) - out[i * nb_coords_rho + j] = eval(bd, rho[j], z[i]); - - return 0; -} diff --git a/brill_data.h b/brill_data.h index 37b338a..e928db9 100644 --- a/brill_data.h +++ b/brill_data.h @@ -18,8 +18,9 @@ #ifndef BRILL_DATA_H #define BRILL_DATA_H -#include +#include #include +#include enum BDQFuncType { /** @@ -32,9 +33,21 @@ enum BDQFuncType { BD_Q_FUNC_EPPLEY, }; +enum BDMetricComponent { + BD_METRIC_COMPONENT_RHORHO, + BD_METRIC_COMPONENT_RHOZ, + BD_METRIC_COMPONENT_RHOPHI, + BD_METRIC_COMPONENT_ZRHO, + BD_METRIC_COMPONENT_ZZ, + BD_METRIC_COMPONENT_ZPHI, + BD_METRIC_COMPONENT_PHIRHO, + BD_METRIC_COMPONENT_PHIZ, + BD_METRIC_COMPONENT_PHIPHI, +}; + typedef struct BDContext { /** - * private data + * Solver internals, not to be accessed by the caller */ void *priv; @@ -86,26 +99,27 @@ typedef struct BDContext { int overdet_z; /** - * The difference between the index of the collocation grid used for the ρ - * direction and nb_coeffs_rho. The 'extra' collocation points furthest from - * the origin are discarded. + * The scaling factor used in the basis functions in the ρ direction. * Defaults to 3. */ - unsigned int colloc_grid_offset_rho; + double basis_scale_factor_rho; /** - * Same as colloc_grid_offset_rho, but for the z direction. + * Same as basis_scale_factor_rho, but for the z direction */ - unsigned int colloc_grid_offset_z; + double basis_scale_factor_z; /** - * The scaling factor used in the basis functions in the ρ direction. - * Defaults to 3. + * A callback that will be used to print diagnostic messages. + * + * Defaults to fprintf(stderr, ...) */ - double basis_scale_factor_rho; + void (*log_callback)(const struct BDContext *bd, int level, + const char *fmt, va_list); + /** - * Same as basis_scale_factor_rho, but for the z direction + * Arbitrary user data, e.g. to be used by the log callback. */ - double basis_scale_factor_z; + void *opaque; /********** * output * @@ -121,7 +135,7 @@ typedef struct BDContext { /** * The number of array elements between two rows. */ - ptrdiff_t stride; + unsigned int stride; } BDContext; /** @@ -155,15 +169,17 @@ int bd_solve(BDContext *bd); * itself, diff_order = { 0, 1 } evaluates ∂ψ/∂z etc. * @param psi the array into which the values of ψ will be written. ψ is evaluated on the grid * formed by the outer product of the rho and z vectors. - * I.e. psi[j * nb_coords_rho + i] = ψ(rho[i], z[j]). The length of psi must be - * nb_coords_rho * nb_coords_z. + * I.e. psi[j * psi_stride + i] = ψ(rho[i], z[j]). The length of psi must be + * at least psi_stride * nb_coords_z. + * @param psi_stride the distance (in double-sized elements) in psi between two elements corresponding + * to the same ρ but one step in z. Must be at least nb_coords_rho. * * @return >= 0 on success, a negative error code on failure. */ -int bd_eval_psi(BDContext *bd, const double *rho, int nb_coords_rho, +int bd_eval_psi(const BDContext *bd, const double *rho, int nb_coords_rho, const double *z, int nb_coords_z, const unsigned int diff_order[2], - double *psi); + double *psi, unsigned int psi_stride); /** * Evaluate the 3-metric γ_ij at the specified rectangular grid (in cylindrical @@ -174,21 +190,28 @@ int bd_eval_psi(BDContext *bd, const double *rho, int nb_coords_rho, * @param nb_coords_rho the number of elements in rho. * @param z the array of z coordinates. * @param nb_coords_z the number of elements in z. - * @param comp the component of the metric to evaluate. - * @param diff_order the order of the derivatives of the metric to evaluate. The first element - * specifies the derivative wrt ρ, the second wrt z. I.e. diff_order = { 0, 0 } - * evaluates the metric itself, diff_order = { 0, 1 } evaluates ∂γ/∂z etc. - * @param out the array into which the values of the metric will be written. The metric - * is evaluated on the grid formed by the outer product of the rho - * and z vectors. I.e. out[j * nb_coords_rho + i] = γ_comp[0]comp[1](rho[i], z[j]). - * The length of out must be nb_coords_rho * nb_coords_z. + * @param nb_comp number of needed components of the metric + * @param comp a nb_comp-sized array specifying the components of the metric to evaluate + * @param diff_order a nb_comp-sized array specifying the order of the derivatives of the + * metric to evaluate. The first element specifies the derivative wrt ρ, + * the second wrt z. I.e. diff_order[i] = { 0, 0 } evaluates the metric + * itself, diff_order[i] = { 0, 1 } evaluates ∂γ/∂z etc. + * @param out a nb_comp-sized array of pointers to the arrays into which the values of the + * metric will be written. Each requested component is evaluated on the grid + * formed by the outer product of the rho and z vectors. I.e. + * out[l][j * out_strides[l] + i] = γ_comp[l](rho[i], z[j]). The length of each + * array in out must be nb_coords_rho * nb_coords_z. + * @param out_strides a nb_comp-sized array of distances (in double-sized elements), for each + * array in out, between two elements corresponding to the same ρ but one + * step in z. Each element in out_strides must be at least nb_coords_rho. * * @return >= 0 on success, a negative error code on failure. */ -int bd_eval_metric(BDContext *bd, const double *rho, int nb_coords_rho, +int bd_eval_metric(const BDContext *bd, const double *rho, int nb_coords_rho, const double *z, int nb_coords_z, - const unsigned int comp[2], const unsigned int diff_order[2], - double *out); + int nb_comp, const enum BDMetricComponent *comp, + const unsigned int (*diff_order)[2], + double **out, unsigned int *out_strides); /** * Evaluate the q function at the specified rectangular grid (in cylindrical @@ -204,13 +227,15 @@ int bd_eval_metric(BDContext *bd, const double *rho, int nb_coords_rho, * evaluates the q function itself, diff_order = { 0, 1 } evaluates ∂q/∂z etc. * @param out the array into which the values of the q function will be written. The q function * is evaluated on the grid formed by the outer product of the rho - * and z vectors. I.e. out[j * nb_coords_rho + i] = γ_comp[0]comp[1](rho[i], z[j]). + * and z vectors. I.e. out[j * out_stride + i] = γ_comp[0]comp[1](rho[i], z[j]). * The length of out must be nb_coords_rho * nb_coords_z. + * @param out_stride the distance (in double-sized elements) in out between two elements corresponding + * to the same ρ but one step in z. Must be at least nb_coords_rho. * * @return >= 0 on success, a negative error code on failure. */ -int bd_eval_q(BDContext *bd, const double *rho, int nb_coords_rho, +int bd_eval_q(const BDContext *bd, const double *rho, int nb_coords_rho, const double *z, int nb_coords_z, const unsigned int diff_order[2], - double *out); + double *out, unsigned int out_stride); #endif /* BRILL_DATA_H */ diff --git a/brill_data.py b/brill_data.py index 26c8efb..398a8d8 100644 --- a/brill_data.py +++ b/brill_data.py @@ -20,6 +20,15 @@ import brill_data import ctypes import numpy as np +BD_METRIC_COMPONENT_RHORHO = 0 +BD_METRIC_COMPONENT_RHOZ = 1 +BD_METRIC_COMPONENT_RHOPHI = 2 +BD_METRIC_COMPONENT_ZRHO = 3 +BD_METRIC_COMPONENT_ZZ = 4 +BD_METRIC_COMPONENT_ZPHI = 5 +BD_METRIC_COMPONENT_PHIRHO = 6 +BD_METRIC_COMPONENT_PHIZ = 7 +BD_METRIC_COMPONENT_PHIPHI = 8 class BrillData(object): @@ -37,10 +46,10 @@ class BrillData(object): ("nb_coeffs_z", ctypes.c_uint), ("overdet_rho", ctypes.c_int), ("overdet_z", ctypes.c_int), - ("colloc_grid_offset_rho", ctypes.c_uint), - ("colloc_grid_offset_z", ctypes.c_uint), ("basis_scale_factor_rho", ctypes.c_double), ("basis_scale_factor_z", ctypes.c_double), + ("log_callback", ctypes.c_void_p), + ("opaque", ctypes.c_void_p), ("psi_minus1_coeffs", ctypes.POINTER(ctypes.c_double)), ("stride", ctypes.c_long)] @@ -54,7 +63,7 @@ class BrillData(object): if ret < 0: raise RuntimeError('Error solving the equation') - self.coeffs = np.copy(np.ctypeslib.as_array(self._bdctx.contents.psi_minus1_coeffs, (self._bdctx.contents.nb_coeffs_rho, self._bdctx.contents.nb_coeffs_z))) + self.coeffs = np.copy(np.ctypeslib.as_array(self._bdctx.contents.psi_minus1_coeffs, (self._bdctx.contents.nb_coeffs_z, self._bdctx.contents.nb_coeffs_rho))) def __del__(self): if self._bdctx: @@ -80,7 +89,7 @@ class BrillData(object): c_z = np.ctypeslib.as_ctypes(z) ret = self._libbd.bd_eval_psi(self._bdctx, c_rho, len(c_rho), - c_z, len(c_z), c_diff_order, c_psi) + c_z, len(c_z), c_diff_order, c_psi, len(c_rho)) if ret < 0: raise RuntimeError('Error evaluating psi') @@ -97,18 +106,22 @@ class BrillData(object): c_diff_order[0] = diff_order[0] c_diff_order[1] = diff_order[1] - c_component = (ctypes.c_uint * 2)() - c_component[0] = component[0] - c_component[1] = component[1] + c_component = (ctypes.c_uint * 1)() + c_component[0] = component - metric = np.empty((z.shape[0], rho.shape[0])) - - c_metric = np.ctypeslib.as_ctypes(metric) c_rho = np.ctypeslib.as_ctypes(rho) c_z = np.ctypeslib.as_ctypes(z) + metric = np.empty((z.shape[0], rho.shape[0])) + + c_metric = (ctypes.POINTER(ctypes.c_double) * 1)() + c_metric[0] = ctypes.cast(np.ctypeslib.as_ctypes(metric), type(c_metric[0])) + + c_metric_strides = (ctypes.c_uint * 1)() + c_metric_strides[0] = len(c_rho) + ret = self._libbd.bd_eval_metric(self._bdctx, c_rho, len(c_rho), - c_z, len(c_z), c_component, c_diff_order, c_metric) + c_z, len(c_z), 1, c_component, c_diff_order, c_metric, c_metric_strides) if ret < 0: raise RuntimeError('Error evaluating the metric') @@ -132,7 +145,7 @@ class BrillData(object): c_z = np.ctypeslib.as_ctypes(z) ret = self._libbd.bd_eval_q(self._bdctx, c_rho, len(c_rho), - c_z, len(c_z), c_diff_order, c_q) + c_z, len(c_z), c_diff_order, c_q, len(c_rho)) if ret < 0: raise RuntimeError('Error evaluating q') diff --git a/eval.c b/eval.c new file mode 100644 index 0000000..deb9579 --- /dev/null +++ b/eval.c @@ -0,0 +1,320 @@ +/* + * Copyright 2014-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 . + */ + +/** + * @file + * evaluation of the solution on a grid + */ + +#include +#include +#include +#include + +#include "brill_data.h" +#include "internal.h" + +int bd_eval_psi(const BDContext *bd, + const double *rho, int nb_coords_rho, + const double *z, int nb_coords_z, + const unsigned int diff_order[2], + double *psi, unsigned int psi_stride) +{ + BDPriv *s = bd->priv; + const double sfr = bd->basis_scale_factor_rho, sfz = bd->basis_scale_factor_z; + double (*eval)(double coord, int idx, double sf); + + double *basis_val_rho = NULL, *basis_val_z = NULL; + + double add = 0.0; + + int ret = 0; + + 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; + } + + if (diff_order[0] == 0 && diff_order[1] == 0) + add = 1.0; + + /* precompute the basis values on the grid points */ + basis_val_rho = malloc(sizeof(*basis_val_rho) * bd->nb_coeffs_rho * nb_coords_rho); + basis_val_z = malloc(sizeof(*basis_val_z) * bd->nb_coeffs_z * nb_coords_z); + if (!basis_val_rho || !basis_val_z) { + ret = -ENOMEM; + goto fail; + } + + switch (diff_order[0]) { + case 0: eval = s->basis->eval; break; + case 1: eval = s->basis->eval_diff1; break; + case 2: eval = s->basis->eval_diff2; break; + } + + for (int i = 0; i < nb_coords_rho; i++) { + double rrho = rho[i]; + for (int j = 0; j < bd->nb_coeffs_rho; j++) + basis_val_rho[i * bd->nb_coeffs_rho + j] = eval(rrho, j, sfr); + } + + switch (diff_order[1]) { + case 0: eval = s->basis->eval; break; + case 1: eval = s->basis->eval_diff1; break; + case 2: eval = s->basis->eval_diff2; break; + } + for (int i = 0; i < nb_coords_z; i++) { + double zz = z[i]; + for (int j = 0; j < bd->nb_coeffs_z; j++) + basis_val_z[i * bd->nb_coeffs_z + j] = eval(zz, j, sfz); + } + + for (int i = 0; i < nb_coords_z; i++) { + double *dst = psi + i * psi_stride; + + for (int j = 0; j < nb_coords_rho; j++) { + double ppsi = add; + + for (int m = 0; m < bd->nb_coeffs_z; m++) + for (int n = 0; n < bd->nb_coeffs_rho; n++) + ppsi += s->coeffs[m * bd->nb_coeffs_rho + n] * basis_val_rho[j * bd->nb_coeffs_rho + n] * basis_val_z[i * bd->nb_coeffs_z + m]; + + *dst++ = ppsi; + } + } + +fail: + free(basis_val_rho); + free(basis_val_z); + + return ret; + +} + +static void 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], + const double *psi, + const double *dpsi_rho, const double *dpsi_z, + const double *d2psi_rho, const double *d2psi_z, const double *d2psi_rho_z, + double *out, unsigned int out_stride) +{ + const BDPriv *s = bd->priv; + + 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; + } + +/* the template for the loop over the grid points */ +#define EVAL_LOOP(DO_EVAL) \ +do { \ + 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 double rrho = rho[j]; \ + const double ppsi = psi[i * nb_coords_rho + j]; \ + const double psi2 = SQR(ppsi); \ + const double psi3 = psi2 * ppsi; \ + const double psi4 = SQR(psi2); \ + double val; \ + \ + DO_EVAL; \ + \ + *dst++ = val; \ + } \ + } \ +} while (0) + +#define GRID(arr) (arr[i * nb_coords_rho + j]) + + if (comp == BD_METRIC_COMPONENT_PHIPHI) { + switch (diff_order[0]) { + case 0: + switch (diff_order[1]) { + case 0: EVAL_LOOP(val = SQR(rrho) * psi4); break; + case 1: EVAL_LOOP(val = 4 * SQR(rrho) * psi3 * GRID(dpsi_z)); break; // ∂_z + case 2: EVAL_LOOP(val = 4 * SQR(rrho) * (GRID(d2psi_z) * psi3 + 3 * psi2 * SQR(GRID(dpsi_z)))); break; // ∂_z ∂_z + } + break; + case 1: + switch (diff_order[1]) { + case 0: EVAL_LOOP(val = 2 * rrho * psi4 + 4 * SQR(rrho) * psi3 * GRID(dpsi_rho)); break; // ∂_ρ + case 1: EVAL_LOOP(val = 8 * rrho * psi3 * GRID(dpsi_z) + 12 * SQR(rrho) * psi2 * GRID(dpsi_z) * GRID(dpsi_rho) + 4 * SQR(rrho) * psi3 * GRID(d2psi_rho_z)); break; // ∂_ρ ∂_z + } + break; + case 2: EVAL_LOOP(val = 4 * SQR(rrho) * psi3 * GRID(d2psi_rho) + 12 * SQR(rrho * ppsi * GRID(dpsi_rho)) + 16 * rrho * psi3 * GRID(dpsi_rho) + 2 * psi4); break; // ∂_ρ ∂_ρ + } + } else { + // γ_ρρ / γ_zz + +/* a wrapper around the actual evaluation expression that provides the q function and its + * 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)); \ + 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; \ + \ + DO_EVAL; \ +} while (0) + + switch (diff_order[0]) { + case 0: + switch (diff_order[1]) { + case 0: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = base, 0, 0, 0, 0, 0)); break; + case 1: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = 2 * base * dz, 0, 1, 0, 0, 0)); break; // ∂_z + case 2: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = 4 * SQR(dz) * base + 2 * base * d2z, 0, 1, 0, 1, 0)); break; // ∂_z ∂_z + } + break; + case 1: + switch (diff_order[1]) { + case 0: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = 2 * base * drho, 1, 0, 0, 0, 0)); break; // ∂_ρ + case 1: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = 4 * base * drho * dz + 2 * base * d2rho_z, 1, 1, 0, 0, 1)); break; // ∂_ρ ∂_z + } + break; + case 2: EVAL_LOOP(DO_EVAL_Q_WRAPPER(val = 4 * SQR(drho) * base + 2 * base * d2rho, 1, 0, 1, 0, 0)); break; // ∂_ρ ∂_ρ + } + } +} + +int bd_eval_metric(const BDContext *bd, const double *rho, int nb_coords_rho, + const double *z, int nb_coords_z, + int nb_comp, const enum BDMetricComponent *comp, + const unsigned int (*diff_order)[2], + double **out, unsigned int *out_strides) +{ + const int nb_coords = nb_coords_rho * nb_coords_z; + double *psi = NULL, *dpsi_rho = NULL, *dpsi_z = NULL, *d2psi_rho = NULL, *d2psi_rho_z = NULL, *d2psi_z = NULL; + int ret = 0; + + /* check the parameters for validity */ + for (int i = 0; i < nb_comp; i++) { + if (comp[i] >= 9 || comp[i] < 0) { + bdi_log(bd, 0, "Invalid component %d: %d\n", i, comp[i]); + return -EINVAL; + } + + if (diff_order[i][0] + diff_order[i][1] > 2) { + bdi_log(bd, 0, "At most second order derivatives are supported.\n"); + return -ENOSYS; + } + } + + /* evaluate the conformal factor and its derivatives as necessary */ +#define EVAL_PSI(arr, order) \ +do { \ + if (arr) \ + break; \ + \ + arr = malloc(nb_coords * sizeof(*arr)); \ + if (!arr) { \ + ret = -ENOMEM; \ + goto fail; \ + } \ + \ + ret = bd_eval_psi(bd, rho, nb_coords_rho, z, nb_coords_z, \ + order, arr, nb_coords_rho); \ + if (ret < 0) \ + goto fail; \ +} while (0) + + for (int i = 0; i < nb_comp; i++) { + if (comp[i] != BD_METRIC_COMPONENT_RHORHO && + comp[i] != BD_METRIC_COMPONENT_ZZ && + comp[i] != BD_METRIC_COMPONENT_PHIPHI) + continue; + + EVAL_PSI(psi, ((unsigned int[2]){ 0, 0 })); + if (diff_order[i][0]) + EVAL_PSI(dpsi_rho, ((unsigned int[2]){ 1, 0 })); + if (diff_order[i][0] > 1) + EVAL_PSI(d2psi_rho, ((unsigned int[2]){ 2, 0 })); + if (diff_order[i][1]) + EVAL_PSI(dpsi_z, ((unsigned int[2]){ 0, 1 })); + if (diff_order[i][1] > 1) + EVAL_PSI(d2psi_z, ((unsigned int[2]){ 0, 2 })); + if (diff_order[i][0] && diff_order[i][1]) + EVAL_PSI(d2psi_rho_z, ((unsigned int[2]){ 1, 1 })); + } + + /* evaluate the requested metric components */ + for (int i = 0; i < nb_comp; i++) + eval_metric(bd, rho, nb_coords_rho, z, nb_coords_z, + comp[i], diff_order[i], + psi, dpsi_rho, dpsi_z, d2psi_rho, d2psi_z, d2psi_rho_z, + out[i], out_strides[i]); + +fail: + free(psi); + free(dpsi_rho); + free(dpsi_z); + free(d2psi_rho); + free(d2psi_rho_z); + free(d2psi_z); + + return ret; +} + +int bd_eval_q(const BDContext *bd, const double *rho, int nb_coords_rho, + const double *z, int nb_coords_z, const unsigned int diff_order[2], + double *out, unsigned int out_stride) +{ + BDPriv *s = bd->priv; + double (*eval)(const struct BDContext *bd, double rho, double z); + + 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; + } + + 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]); + } + + return 0; +} diff --git a/init.c b/init.c new file mode 100644 index 0000000..759088f --- /dev/null +++ b/init.c @@ -0,0 +1,133 @@ +/* + * Copyright 2014-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 . + */ + +#include +#include +#include +#include +#include + +#include "brill_data.h" +#include "internal.h" + +static int brill_init_check_options(BDContext *bd) +{ + BDPriv *s = bd->priv; + + 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; + } + + if (!bd->nb_coeffs_z) + bd->nb_coeffs_z = bd->nb_coeffs_rho; + + s->basis = &bdi_sb_even_basis; + + s->nb_coeffs = bd->nb_coeffs_rho * bd->nb_coeffs_z; + + s->nb_colloc_points_rho = bd->nb_coeffs_rho + bd->overdet_rho; + s->nb_colloc_points_z = bd->nb_coeffs_z + bd->overdet_z; + + s->nb_colloc_points = s->nb_colloc_points_rho * s->nb_colloc_points_z; + + s->colloc_grid_order_rho = s->nb_colloc_points_rho; + s->colloc_grid_order_z = s->nb_colloc_points_z; + + return 0; +} + +int bd_solve(BDContext *bd) +{ + BDPriv *s = bd->priv; + int ret; + + if (bd->psi_minus1_coeffs) { + bdi_log(bd, 0, "Solve called multiple times on the same context.\n"); + return -EINVAL; + } + + ret = brill_init_check_options(bd); + if (ret < 0) + return ret; + + ret = bdi_solve(bd); + if (ret < 0) + return ret; + + bd->psi_minus1_coeffs = s->coeffs; + bd->stride = bd->nb_coeffs_rho; + + return 0; +} + +BDContext *bd_context_alloc(void) +{ + BDContext *bd = calloc(1, sizeof(*bd)); + + if (!bd) + return NULL; + + bd->q_func_type = BD_Q_FUNC_GUNDLACH; + bd->amplitude = 1.0; + bd->eppley_n = 5; + + bd->nb_coeffs_rho = 80; + bd->nb_coeffs_z = 0; + + bd->basis_scale_factor_rho = 3.0; + bd->basis_scale_factor_z = 3.0; + + bd->log_callback = bdi_log_default_callback; + + bd->priv = calloc(1, sizeof(BDPriv)); + if (!bd->priv) { + free(bd); + return NULL; + } + + return bd; +} + +void bd_context_free(BDContext **pbd) +{ + BDContext *bd = *pbd; + BDPriv *s; + + if (!bd) + return; + + s = bd->priv; + + free(s->coeffs); + + free(bd->priv); + free(bd); + *pbd = NULL; +} diff --git a/internal.h b/internal.h new file mode 100644 index 0000000..87490e1 --- /dev/null +++ b/internal.h @@ -0,0 +1,83 @@ +/* + * Copyright 2014-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_INTERNAL_H +#define BRILL_DATA_INTERNAL_H + +#include "brill_data.h" + +#define MAX(x, y) ((x) > (y) ? (x) : (y)) +#define SQR(x) ((x) * (x)) + +#define ARRAY_ELEMS(x) (sizeof(x) / sizeof(x[0])) + +/* + * small number to avoid r=0 singularities + */ +#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 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; + const QFunc *q_func; + + int nb_colloc_points; + int nb_colloc_points_rho; + int nb_colloc_points_z; + + int colloc_grid_order_rho; + int colloc_grid_order_z; + int nb_coeffs; + + double *coeffs; +} 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); + +void bdi_log(const BDContext *bd, int level, const char *fmt, ...); +void bdi_log_default_callback(const BDContext *bd, int level, + const char *fmt, va_list vl); + +#endif // BRILL_DATA_INTERNAL_H diff --git a/log.c b/log.c new file mode 100644 index 0000000..9b82f34 --- /dev/null +++ b/log.c @@ -0,0 +1,41 @@ +/* + * Copyright 2014-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 . + */ + +/** + * @file + * logging code + */ + +#include +#include + +#include "brill_data.h" +#include "internal.h" + +void bdi_log_default_callback(const BDContext *bd, int level, + const char *fmt, va_list vl) +{ + vfprintf(stderr, fmt, vl); +} + +void bdi_log(const BDContext *bd, int level, const char *fmt, ...) +{ + va_list vl; + va_start(vl, fmt); + bd->log_callback(bd, level, fmt, vl); + va_end(vl); +} diff --git a/qfunc.c b/qfunc.c new file mode 100644 index 0000000..fea4af9 --- /dev/null +++ b/qfunc.c @@ -0,0 +1,161 @@ +/* + * Copyright 2014-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 . + */ + +/** + * @file + * definitions of the q functions in the exponential in theBrill data + */ + +#include + +#include "brill_data.h" +#include "internal.h" + +static double q_gundlach(const BDContext *bd, double rho, double z) +{ + return bd->amplitude * SQR(rho) * exp(- (SQR(rho) + SQR(z))); +} + +static double dq_rho_gundlach(const BDContext *bd, double rho, double z) +{ + return bd->amplitude * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho)); +} + +static double d2q_rho_gundlach(const BDContext *bd, double rho, double z) +{ + double rho2 = SQR(rho); + return bd->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) +{ + return -bd->amplitude * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho)); +} + +static double dq_z_gundlach(const BDContext *bd, double rho, double z) +{ + return - bd->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z)); +} + +static double d2q_z_gundlach(const BDContext *bd, double rho, double z) +{ + return bd->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 double q_eppley(const BDContext *bd, double rho, double z) +{ + double r2 = SQR(rho) + SQR(z); + return bd->amplitude * SQR(rho) / (1.0 + pow(r2, bd->eppley_n / 2.0)); +} + +static double dq_rho_eppley(const BDContext *bd, double rho, double z) +{ + double A = bd->amplitude; + double n = bd->eppley_n; + double rho2 = SQR(rho); + double z2 = SQR(z); + + double r2 = rho2 + z2; + double r = sqrt(r2); + double rnm2 = pow(r, n - 2); + + 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) +{ + double A = bd->amplitude; + double n = bd->eppley_n; + double rho2 = SQR(rho); + double z2 = SQR(z); + + double r2 = rho2 + z2; + double r = sqrt(r2); + double rnm2 = pow(r, n - 2); + + return - A * n * rho2 * z * rnm2 / SQR(1 + rnm2 * r2); +} + +static double d2q_rho_eppley(const BDContext *bd, double rho, double z) +{ + double A = bd->amplitude; + double n = bd->eppley_n; + double rho2 = SQR(rho); + double z2 = SQR(z); + + double r2 = rho2 + z2; + double r = sqrt(r2); + double rnm4 = pow(r, n - 4); + double rn = rnm4 * SQR(r2); + double rnp1 = 1 + rn; + + 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) +{ + double A = bd->amplitude; + double n = bd->eppley_n; + double rho2 = SQR(rho); + double z2 = SQR(z); + + double r2 = rho2 + z2; + double r = sqrt(r2); + double rnm4 = pow(r, n - 4); + double rn = rnm4 * SQR(r2); + double rnp1 = 1 + rn; + + 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) +{ + double A = bd->amplitude; + double n = bd->eppley_n; + double rho2 = SQR(rho); + double z2 = SQR(z); + + double r2 = rho2 + z2; + double r = sqrt(r2); + double rnm4 = pow(r, n - 4); + double rn = rnm4 * SQR(r2); + double rnp1 = 1 + rn; + + return A * n * rho * z * rnm4 * (rn * (n * rho2 - 2 * z2) - 2 * z2 - n * rho2) / (rnp1 * SQR(rnp1)); +} + +// 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, +}; diff --git a/solve.c b/solve.c new file mode 100644 index 0000000..ab8fc39 --- /dev/null +++ b/solve.c @@ -0,0 +1,230 @@ +/* + * Copyright 2014-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 . + */ + +/** + * @file + * the actual pseudo-spectral solver code + */ + +#include +#include + +#include +#include +#include +#include + +#include + +#include "brill_data.h" +#include "internal.h" + +static int brill_construct_matrix(const BDContext *bd, double *mat, + double *rhs) +{ + BDPriv *s = bd->priv; + const double sfr = bd->basis_scale_factor_rho, sfz = bd->basis_scale_factor_z; + + double *basis_val, *basis_dval, *basis_d2val; + int idx_coeff_rho, idx_coeff_z, idx_grid_rho, idx_grid_z; + + double *basis_r[3]; + double *basis_z[3]; + + /* precompute the basis values we will need */ + for (int i = 0; i < ARRAY_ELEMS(basis_r); i++) + basis_r[i] = malloc(sizeof(*basis_r[i]) * s->nb_colloc_points_rho * bd->nb_coeffs_rho); + for (int i = 0; i < s->nb_colloc_points_rho; i++) { + double coord = s->basis->colloc_point(s->colloc_grid_order_rho, i, sfr); + for (int j = 0; j < bd->nb_coeffs_rho; j++) { + basis_r[0][i * bd->nb_coeffs_rho + j] = s->basis->eval(coord, j, sfr); + basis_r[1][i * bd->nb_coeffs_rho + j] = s->basis->eval_diff1(coord, j, sfr); + basis_r[2][i * bd->nb_coeffs_rho + j] = s->basis->eval_diff2(coord, j, sfr); + } + } + + for (int i = 0; i < ARRAY_ELEMS(basis_z); i++) + basis_z[i] = malloc(sizeof(*basis_z[i]) * s->nb_colloc_points_z * bd->nb_coeffs_z); + for (int i = 0; i < s->nb_colloc_points_z; i++) { + double coord = s->basis->colloc_point(s->colloc_grid_order_z, i, sfr); + for (int j = 0; j < bd->nb_coeffs_z; j++) { + basis_z[0][i * bd->nb_coeffs_z + j] = s->basis->eval(coord, j, sfr); + basis_z[1][i * bd->nb_coeffs_z + j] = s->basis->eval_diff1(coord, j, sfr); + basis_z[2][i * bd->nb_coeffs_z + j] = s->basis->eval_diff2(coord, j, sfr); + } + } + +#define BASIS_RHO (basis_r[0][idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho]) +#define DBASIS_RHO (basis_r[1][idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho]) +#define D2BASIS_RHO (basis_r[2][idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho]) +#define BASIS_Z (basis_z[0][idx_grid_z * bd->nb_coeffs_z + idx_coeff_z]) +#define DBASIS_Z (basis_z[1][idx_grid_z * bd->nb_coeffs_z + idx_coeff_z]) +#define D2BASIS_Z (basis_z[2][idx_grid_z * bd->nb_coeffs_z + idx_coeff_z]) + + for (idx_grid_z = 0; idx_grid_z < s->nb_colloc_points_z; idx_grid_z++) { + double z_physical = s->basis->colloc_point(s->colloc_grid_order_z, idx_grid_z, sfz); + + for (idx_grid_rho = 0; idx_grid_rho < s->nb_colloc_points_rho; idx_grid_rho++) { + double x_physical = s->basis->colloc_point(s->colloc_grid_order_rho, idx_grid_rho, sfr); + double d2q = s->q_func->d2q_rho(bd, x_physical, z_physical) + s->q_func->d2q_z(bd, x_physical, z_physical); + int idx_grid = idx_grid_z * s->nb_colloc_points_rho + idx_grid_rho; + + for (idx_coeff_z = 0; idx_coeff_z < bd->nb_coeffs_z; idx_coeff_z++) + for (idx_coeff_rho = 0; idx_coeff_rho < bd->nb_coeffs_rho; idx_coeff_rho++) { + int idx_coeff = idx_coeff_z * bd->nb_coeffs_rho + idx_coeff_rho; + + double val = D2BASIS_RHO * BASIS_Z + D2BASIS_Z * BASIS_RHO + BASIS_RHO * BASIS_Z * 0.25 * d2q; + if (fabs(x_physical) > EPS) + val += DBASIS_RHO * BASIS_Z / fabs(x_physical); + else + val += D2BASIS_RHO * BASIS_Z; + + mat[idx_grid + s->nb_colloc_points * idx_coeff] = val; + } + rhs[idx_grid] = -0.25 * d2q; + } + } + + for (int i = 0; i < ARRAY_ELEMS(basis_r); i++) + free(basis_r[i]); + for (int i = 0; i < ARRAY_ELEMS(basis_z); i++) + free(basis_z[i]); + + return 0; +} + +static int brill_solve_linear(const BDContext *bd, double *mat, double **px, double **prhs) +{ + const BDPriv *s = bd->priv; + const int stride = s->nb_coeffs; + int *ipiv; + double *mat_f; + double cond, ferr, berr, rpivot; + + double *x = *px; + double *rhs = *prhs; + char equed = 'N'; + int ret = 0; + + ipiv = malloc(stride * sizeof(*ipiv)); + mat_f = malloc(SQR(stride) * sizeof(*mat_f)); + if (!ipiv || !mat_f) { + ret = -ENOMEM; + goto fail; + } + + LAPACKE_dgesvx(LAPACK_COL_MAJOR, 'N', 'N', stride, 1, + mat, stride, mat_f, stride, ipiv, &equed, + NULL, NULL, rhs, stride, x, stride, + &cond, &ferr, &berr, &rpivot); + + bdi_log(bd, 0, "LU factorization solution to a %zdx%zd matrix: " + "condition number %16.16g; forward error %16.16g backward error %16.16g\n", + stride, stride, cond, ferr, berr); +fail: + free(mat_f); + free(ipiv); + + return ret; +} + +static int brill_solve_svd(const BDContext *bd, double *mat, + double **px, double **prhs) +{ + const BDPriv *s = bd->priv; + double *sv; + int rank; + + double *x = *px; + double *rhs = *prhs; + + sv = malloc(sizeof(*sv) * s->nb_coeffs); + if (!sv) + return -ENOMEM; + + LAPACKE_dgelsd(LAPACK_COL_MAJOR, s->nb_colloc_points, s->nb_coeffs, 1, + mat, s->nb_colloc_points, rhs, + MAX(s->nb_coeffs, s->nb_colloc_points), + sv, -1, &rank); + + bdi_log(bd, 0, "Least squares SVD solution to a %zdx%zd matrix: " + "rank %d, condition number %16.16g\n", s->nb_colloc_points, s->nb_coeffs, + rank, sv[s->nb_coeffs - 1] / sv[0]); + + free(sv); + + *px = rhs; + *prhs = x; + + return 0; +} + +/* + * Solve the equation + * Δψ + ¼ ψ (∂²q/∂r² + ∂²q/∂z²) = 0 + * for the coefficients of spectral approximation of ψ: + * ψ(r, z) = 1 + ΣaᵢⱼTᵢ(r)Tⱼ(z) + * where i = { 0, ... , bd->nb_coeffs_rho }; + * j = { 0, ... , bd->nb_coeffs_z }; + * q(r, z) and Tᵢ(x) are defined by bd->q_func, bd->laplace_q_func and + * bd->basis. + * + * The cofficients are exported in bd->priv->coeffs + */ +int bdi_solve(BDContext *bd) +{ + BDPriv *s = bd->priv; + + const int vecsize = MAX(s->nb_coeffs, s->nb_colloc_points); + + double *mat = NULL; + double *rhs = NULL, *coeffs = NULL; + + int ret = 0; + + /* allocate and fill the pseudospectral matrix */ + mat = malloc(sizeof(*mat) * s->nb_coeffs * s->nb_colloc_points); + coeffs = malloc(sizeof(*coeffs) * vecsize); + rhs = malloc(sizeof(*rhs) * vecsize); + if (!mat || !coeffs || !rhs) { + ret = -ENOMEM; + goto fail; + } + + /* fill the matrix */ + ret = brill_construct_matrix(bd, mat, rhs); + if (ret < 0) + goto fail; + + /* solve for the coeffs */ + if (s->nb_colloc_points > s->nb_coeffs) + ret = brill_solve_svd(bd, mat, &coeffs, &rhs); + else + ret = brill_solve_linear(bd, mat, &coeffs, &rhs); + + /* export the result */ + s->coeffs = coeffs; + +fail: + if (ret < 0) + free(coeffs); + + free(rhs); + free(mat); + + return ret; +} -- cgit v1.2.3