1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
|
/*
* 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 <stdint.h>
#include <stddef.h>
enum BDQFuncType {
/**
* q(ρ, z) = A ρ^2 exp(-ρ^2 - z^2)
*/
BD_Q_FUNC_GUNDLACH,
/**
* q(ρ, z) = A ρ^2 / (1 + r^n); where r^2 = ρ^2 + z^2
*/
BD_Q_FUNC_EPPLEY,
};
typedef struct BDContext {
/**
* private data
*/
void *priv;
/*******************************
* options, set by the caller *
*******************************/
/**
* The choice of the q function in the exponent.
* Defaults to BD_Q_FUNC_GUNDLACH.
*/
enum BDQFuncType q_func_type;
/**
* The amplitude in the q function.
* Defaults to 1.
*/
double amplitude;
/**
* For BD_Q_FUNC_EPPLEY, the power in the denominator.
* Must be >= 4.
* Defaults to 5.
*/
unsigned int eppley_n;
/* the solver parameters */
/**
* The number of basis functions in the ρ direction.
* Defaults to 80.
*/
unsigned int nb_coeffs_rho;
/**
* The number of basis functions in the z direction. 0 means the same as
* nb_coeffs_rho. Defaults to 0.
* Using values other than 0 or nb_coeffs_rho is unsupported for now.
*/
unsigned int nb_coeffs_z;
/**
* The difference between the number of collocation points and the number
* of basis functions in the rho direction. I.e. the number of collocation
* points used will be nb_coeffs_rho + overdet_rho.
* Defaults to 0;
*/
int overdet_rho;
/**
* Same as overdet_rho, but for the z direction.
*/
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.
* Defaults to 3.
*/
unsigned int colloc_grid_offset_rho;
/**
* Same as colloc_grid_offset_rho, but for the z direction.
*/
unsigned int colloc_grid_offset_z;
/**
* The scaling factor used in the basis functions in the ρ direction.
* Defaults to 3.
*/
double basis_scale_factor_rho;
/**
* Same as basis_scale_factor_rho, but for the z direction
*/
double basis_scale_factor_z;
/**********
* output *
**********/
/**
* The coefficients of the solution expanded in the basis.
* The ρ index increases along rows, z along columns.
*
* The data is owned and managed by this library and is read-only for
* the caller.
*/
double *psi_minus1_coeffs;
/**
* The number of array elements between two rows.
*/
ptrdiff_t stride;
} BDContext;
/**
* Allocate and initialize the solver.
*/
BDContext *bd_context_alloc(void);
/**
* Free the solver and everything associated with it.
*/
void bd_context_free(BDContext **bd);
/**
* 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 bd_solve(BDContext *bd);
/**
* Evaluate the conformal factor ψ or its derivatives at the specified rectangular grid.
*
* @param bd 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 diff_order the order of the derivatives of ψ to evaluate. The first element specifies the
* derivative wrt ρ, the second wrt z. I.e. diff_order = { 0, 0 } evaluates ψ
* 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.
*
* @return >= 0 on success, a negative error code on failure.
*/
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);
/**
* Evaluate the 3-metric γ_ij at the specified rectangular grid (in cylindrical
* coordinates { ρ, z, φ }).
*
* @param bd 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 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.
*
* @return >= 0 on success, a negative error code on failure.
*/
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);
|
