aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAnton Khirnov <anton@khirnov.net>2015-04-14 16:00:24 +0200
committerAnton Khirnov <anton@khirnov.net>2015-04-14 16:02:58 +0200
commit41f29cbf3076db51b96f240d27d432ef31b8aa71 (patch)
treeaedc8f140385f4c94054b561e09358a04ffbc296
parent307d35ed2bc4e073ae3a6ceac3dcb7a74f25e98a (diff)
A major rewrite.
Split the code into multiple files, drop the GSL dependency, introduce configurable logging, random cleanups.
-rw-r--r--Makefile13
-rw-r--r--README6
-rw-r--r--basis.c74
-rw-r--r--brill_data.c843
-rw-r--r--brill_data.h89
-rw-r--r--brill_data.py37
-rw-r--r--eval.c320
-rw-r--r--init.c133
-rw-r--r--internal.h83
-rw-r--r--log.c41
-rw-r--r--qfunc.c161
-rw-r--r--solve.c230
12 files changed, 1136 insertions, 894 deletions
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 <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file
+ * definitions of the basis functions for the spectral method
+ */
+
+#include <math.h>
+
+#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 <anton@khirnov.net>
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- */
-
-#include <errno.h>
-#include <limits.h>
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include <gsl/gsl_blas.h>
-#include <gsl/gsl_linalg.h>
-
-#include <lapacke.h>
-
-#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 <stdint.h>
+#include <stdarg.h>
#include <stddef.h>
+#include <stdint.h>
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 <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file
+ * evaluation of the solution on a grid
+ */
+
+#include <errno.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#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 <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <errno.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+
+#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 <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#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 <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file
+ * logging code
+ */
+
+#include <stdarg.h>
+#include <stdio.h>
+
+#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 <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file
+ * definitions of the q functions in the exponential in theBrill data
+ */
+
+#include <math.h>
+
+#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 <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file
+ * the actual pseudo-spectral solver code
+ */
+
+#include <cblas.h>
+#include <lapacke.h>
+
+#include <errno.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <stdio.h>
+
+#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;
+}