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
|
/*
* Copyright 2016 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 template for the loop over the grid points
*/
#define FUNC3(a, b, c) a ## _ ## b ## _ ## c
#define FUNC2(a, b, c) FUNC3(a, b, c)
#define FUNC(x) FUNC2(x, COMP, DIFF)
#define DIFF_RHO (DIFF % 3)
#define DIFF_Z (DIFF / 3)
static void FUNC(eval_metric)(const BDContext *bd, double *out, ptrdiff_t out_stride,
int nb_coords_rho, const double *rho,
int nb_coords_z, const double *z,
double *psi[9], double *q[9])
{
BDPriv *s = bd->priv;
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 int idx = i * nb_coords_rho + j;
const double rrho = rho[j];
const double ppsi = psi[0][idx];
const double psi2 = SQR(ppsi);
const double psi3 = psi2 * ppsi;
const double psi4 = SQR(psi2);
const double base = psi4 * exp(2 * q[0][idx]);
const double dqrho = DIFF_RHO > 0 ? q[1][idx] + 2 * psi[1][idx] / ppsi : 0;
const double dqz = DIFF_Z > 0 ? q[3][idx] + 2 * psi[3][idx] / ppsi : 0;
const double dq2rho = DIFF_RHO > 1 ? q[2][idx] + 2 * psi[2][idx] / ppsi - 2 * SQR(psi[1][idx] / ppsi) : 0;
const double dq2z = DIFF_Z > 1 ? q[6][idx] + 2 * psi[6][idx] / ppsi - 2 * SQR(psi[3][idx] / ppsi) : 0;
const double dq2rho_z = DIFF_RHO && DIFF_Z ? q[4][idx] + 2 * psi[4][idx] / ppsi - 2 * psi[1][idx] * psi[3][idx] / psi2 : 0;
double val;
DO_EVAL;
*dst++ = val;
}
}
}
#undef FUNC3
#undef FUNC2
#undef FUNC
|
