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
|
#include <cassert>
#include <cstdio>
#include <vector>
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
#include <bin_bh.h>
using namespace std;
extern "C"
void ID_Bin_BH_initialise (CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
CCTK_INFO ("Setting up LORENE Bin_BH initial data");
CCTK_INFO ("Setting up coordinates");
int const npoints = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
vector<double> xx(npoints), yy(npoints), zz(npoints);
#pragma omp parallel for
for (int i=0; i<npoints; ++i) {
xx[i] = x[i];
yy[i] = y[i];
zz[i] = z[i];
}
CCTK_VInfo (CCTK_THORNSTRING, "Reading from file \"%s\"", filename);
Bin_BH bin_bh (npoints, &xx[0], &yy[0], &zz[0], 1, filename);
CCTK_VInfo (CCTK_THORNSTRING, "Omega [1/a]: %g", bin_bh.omega);
CCTK_VInfo (CCTK_THORNSTRING, "dist [a]: %g", bin_bh.dist);
CCTK_VInfo (CCTK_THORNSTRING, "radius1 [a]: %g", 1.0);
CCTK_VInfo (CCTK_THORNSTRING, "radius2 [a]: %g", bin_bh.radius2);
assert (bin_bh.np == npoints);
CCTK_INFO ("Filling in Cactus grid points");
#pragma omp parallel for
for (int i=0; i<npoints; ++i) {
alp[i] = bin_bh.nnn[i];
betax[i] = bin_bh.beta_x[i];
betay[i] = bin_bh.beta_y[i];
betaz[i] = bin_bh.beta_z[i];
CCTK_REAL g[3][3];
g[0][0] = bin_bh.g_xx[i];
g[0][1] = bin_bh.g_xy[i];
g[0][2] = bin_bh.g_xz[i];
g[1][1] = bin_bh.g_yy[i];
g[1][2] = bin_bh.g_yz[i];
g[2][2] = bin_bh.g_zz[i];
g[1][0] = g[0][1];
g[2][0] = g[0][2];
g[2][1] = g[1][2];
CCTK_REAL ku[3][3];
ku[0][0] = bin_bh.k_xx[i];
ku[0][1] = bin_bh.k_xy[i];
ku[0][2] = bin_bh.k_xz[i];
ku[1][1] = bin_bh.k_yy[i];
ku[1][2] = bin_bh.k_yz[i];
ku[2][2] = bin_bh.k_zz[i];
ku[1][0] = ku[0][1];
ku[2][0] = ku[0][2];
ku[2][1] = ku[1][2];
CCTK_REAL k[3][3];
for (int a=0; a<3; ++a) {
for (int b=0; b<3; ++b) {
k[a][b] = 0.0;
for (int c=0; c<3; ++c) {
for (int d=0; d<3; ++d) {
k[a][b] += g[a][c] * g[b][d] * ku[c][d];
}
}
}
}
gxx[i] = g[0][0];
gxy[i] = g[0][1];
gxz[i] = g[0][2];
gyy[i] = g[1][1];
gyz[i] = g[1][2];
gzz[i] = g[2][2];
kxx[i] = k[0][0];
kxy[i] = k[0][1];
kxz[i] = k[0][2];
kyy[i] = k[1][1];
kyz[i] = k[1][2];
kzz[i] = k[2][2];
} // for i
CCTK_INFO ("Done.");
}
|
