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
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
|
#include <inttypes.h>
#include <cl.h>
#include "cctk.h"
#define SQR(x) ((x) * (x))
#define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
#define MAX(x, y) ((x) > (y) ? (x) : (y))
#define MIN(x, y) ((x) > (y) ? (y) : (x))
#define ARRAY_ELEMS(arr) (sizeof(arr) / sizeof(*arr))
/*
* small number to avoid r=0 singularities
*/
#define EPS 1E-08
#define SCALE_FACTOR scale_factor
/* indices (in our code, not cactus structs) of the grid functions which we'll need to
* interpolate on the pseudospectral grid */
enum MetricVars {
GTXX = 0,
GTYY,
GTZZ,
GTXY,
GTXZ,
GTYZ,
PHI,
ATXX,
ATYY,
ATZZ,
ATXY,
ATXZ,
ATYZ,
K,
XTX,
XTY,
XTZ,
BETAX,
BETAY,
BETAZ,
NB_METRIC_VARS,
};
/* indices of the interpolated values of the above grid functions and their derivatives */
enum InterpMetricVars {
I_GTXX = 0,
I_GTXX_DX,
I_GTXX_DY,
I_GTXX_DZ,
I_GTYY,
I_GTYY_DX,
I_GTYY_DY,
I_GTYY_DZ,
I_GTZZ,
I_GTZZ_DX,
I_GTZZ_DY,
I_GTZZ_DZ,
I_GTXY,
I_GTXY_DX,
I_GTXY_DY,
I_GTXY_DZ,
I_GTXZ,
I_GTXZ_DX,
I_GTXZ_DY,
I_GTXZ_DZ,
I_GTYZ,
I_GTYZ_DX,
I_GTYZ_DY,
I_GTYZ_DZ,
I_PHI,
I_PHI_DX,
I_PHI_DY,
I_PHI_DZ,
I_ATXX,
I_ATYY,
I_ATZZ,
I_ATXY,
I_ATXZ,
I_ATYZ,
I_K,
I_K_DX,
I_K_DY,
I_K_DZ,
I_XTX,
I_XTY,
I_XTZ,
I_BETAX,
I_BETAY,
I_BETAZ,
NB_INTERP_VARS,
};
/* a set of basis functions */
typedef struct BasisSet {
/* evaluate the idx-th basis function at the specified point*/
double (*eval) (double coord, int idx);
/* evaluate the first derivative of the idx-th basis function at the specified point*/
double (*eval_diff1)(double coord, int idx);
/* evaluate the second derivative of the idx-th basis function at the specified point*/
double (*eval_diff2)(double coord, int idx);
/**
* 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);
} BasisSet;
extern const BasisSet msa_cheb_basis;
extern const BasisSet msa_full_basis;
extern const BasisSet msa_tb_even_basis;
extern const BasisSet msa_sb_even_basis;
extern double scale_factor;
/* precomputed values for a given refined grid */
typedef struct CoordPatch {
CCTK_REAL origin[3];
CCTK_REAL delta[3];
CCTK_INT size[3];
// basis values on the grid
double *basis_val_r;
double *basis_val_z;
double *transform_z;
double *one;
} CoordPatch;
/* state and scratch storage for the BiCGSTAB solver */
typedef struct BiCGSTABContext {
double *p, *v, *y, *z, *t;
double *res, *res0;
double *k;
cl_mem cl_p, cl_v, cl_y, cl_z, cl_t;
cl_mem cl_res, cl_res0;
cl_mem cl_k, cl_mat;
cl_mem cl_rho, cl_alpha, cl_beta, cl_omega, cl_omega1;
cl_mem cl_tmp, cl_tmp1;
int64_t solve_total;
int64_t iter_total;
int64_t time_total;
} BiCGSTABContext;
typedef struct MaximalSlicingContext {
cGH *gh;
const BasisSet *basis;
BiCGSTABContext bicgstab;
int steps_since_inverse;
int64_t lu_solves_total;
int64_t lu_solves_time;
// the grid of collocation points
CCTK_REAL *grid_x;
CCTK_REAL *grid_z;
// interpolation parameters
int coord_system;
int interp_operator;
int interp_params;
CCTK_REAL *interp_coords[3];
int interp_vars_indices[NB_METRIC_VARS];
CCTK_REAL *interp_values[NB_INTERP_VARS];
CCTK_INT interp_value_codes[NB_INTERP_VARS];
CCTK_REAL *metric_u[6];
CCTK_REAL *kij_kij;
CCTK_REAL *trk;
int nb_coeffs_x;
int nb_coeffs_z;
int nb_coeffs;
int nb_colloc_points_x;
int nb_colloc_points_z;
int nb_colloc_points;
int colloc_grid_order_x;
int colloc_grid_order_z;
double *mat;
double *mat_f;
double *rhs;
double *coeffs;
int *ipiv;
double *basis_x_val;
double *basis_x_dval;
double *basis_x_d2val;
double *basis_z_val;
double *basis_z_dval;
double *basis_z_d2val;
double *basis_val_00;
double *basis_val_20;
double *basis_val_02;
double *basis_val_11;
double *basis_val_10;
double *basis_val_01;
CoordPatch *patches;
int nb_patches;
// OpenCL / CLBLAS stuff
cl_context cl_ctx;
cl_command_queue cl_queue;
cl_mem ocl_coeffs;
} MaximalSlicingContext;
int msa_maximal_solve(MaximalSlicingContext *ms);
#include <sys/time.h>
static inline int64_t gettime(void)
{
struct timeval tv;
gettimeofday(&tv, NULL);
return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
}
|
