diff options
Diffstat (limited to 'teukolsky_data.h')
| -rw-r--r-- | teukolsky_data.h | 165 |
1 files changed, 165 insertions, 0 deletions
diff --git a/teukolsky_data.h b/teukolsky_data.h new file mode 100644 index 0000000..78bf824 --- /dev/null +++ b/teukolsky_data.h @@ -0,0 +1,165 @@ +/* + * Copyright 2017 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 TEUKOLSKY_DATA_H +#define TEUKOLSKY_DATA_H + +#include <stdarg.h> +#include <stddef.h> +#include <stdint.h> + +/** + * API usage: + * + * First, allocate the solver context with td_context_alloc(). All interaction + * with the solver is done through this context. + * + * Fill any fields in the context that are described as settable by the caller. + * Call td_solve() to solve the equation determined by the option values. + * + * The metric/extrinsic curvature can be evaluated with td_eval_metric()/td_eval_curv(). + * + * Finally, free the solver context with td_context_free(). + */ + +typedef struct TDContext { + /** + * Solver internals, not to be accessed by the caller + */ + void *priv; + + /******************************** + * options, set by the caller * + ********************************/ + + /** + * A callback that will be used to print diagnostic messages. + * + * Defaults to fprintf(stderr, ...) + */ + void (*log_callback)(const struct TDContext *td, int level, + const char *fmt, va_list); + + /** + * Arbitrary user data, e.g. to be used by the log callback. + */ + void *opaque; + + /******************************** + * initial data parameters * + ********************************/ + + /** + * The amplitude in the I function. + * Defaults to 1. + */ + double amplitude; + + /******************************** + * solver options * + ********************************/ + + /** + * The number of basis functions in each direction. + * [0] - radial, [1] - angular + */ + unsigned int nb_coeffs[2]; + + /** + * The scaling factor used in the basis functions. + * [0] - radial, [1] - angular + */ + double basis_scale_factor[2]; + + /** + * Maximum number of Newton iterations. + */ + unsigned int max_iter; + + /** + * Absolute tolerance. The solver is deemed to have converged + * after maximum difference between iterations is below this bound. + */ + double atol; + + /** + * Number of threads to use for parallelization. Set to 0 to automatically + * detect according to the number of CPU cores. + */ + unsigned int nb_threads; + + double *coeffs[3]; +} TDContext; + +/** + * Allocate and initialize the solver. + */ +TDContext *td_context_alloc(void); + +/** + * Free the solver and everything associated with it. + */ +void td_context_free(TDContext **td); + +/** + * Solve the equation for the conformal factor ψ and export the expansion coefficients in the + * context. + * + * @return >= 0 on success, a negative error code on failure + */ +int td_solve(TDContext *td, double *coeffs_init[3]); + +/** + * Evaluate the 3-metric γ_ij at the specified rectangular grid (in cylindrical + * coordinates { ρ, z, φ }). + * + * @param td the solver context + * @param rho the array of ρ coordinates. + * @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 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 td_eval_psi(const TDContext *td, + size_t nb_coords, const double *r, const double *theta, + const unsigned int diff_order[2], + double *out); +int td_eval_krr(const TDContext *td, + size_t nb_coords, const double *r, const double *theta, + const unsigned int diff_order[2], + double *out); +int td_eval_kpp(const TDContext *td, + size_t nb_coords, const double *r, const double *theta, + const unsigned int diff_order[2], + double *out); + +#endif /* TEUKOLSKY_DATA_H */ |