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
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
|
#include <math.h>
#include <cctk.h>
#ifdef CCTK_CXX_RESTRICT
# undef restrict
# define restrict CCTK_CXX_RESTRICT
#endif
#define KRANC_C
extern "C" {
#include <GenericFD.h>
}
// Adapted from BSSN_MoL's files NewRad.F and newrad.h
static
void newrad_kernel (cGH const* restrict const cctkGH,
int const* restrict const bmin,
int const* restrict const bmax,
int const* restrict const dir,
CCTK_REAL const* restrict const var,
CCTK_REAL * restrict const rhs,
CCTK_REAL const* restrict const x,
CCTK_REAL const* restrict const y,
CCTK_REAL const* restrict const z,
CCTK_REAL const* restrict const r,
CCTK_REAL const& var0,
CCTK_REAL const& v0,
CCTK_REAL const& radpower)
{
int const ni = cctkGH->cctk_lsh[0];
int const nj = cctkGH->cctk_lsh[1];
int const nk = cctkGH->cctk_lsh[2];
int const si = -dir[0];
int const sj = -dir[1];
int const sk = -dir[2];
int const di = 1;
int const dj = ni;
int const dk = ni*nj;
CCTK_REAL const dx = cctkGH->cctk_delta_space[0] / cctkGH->cctk_levfac[0];
CCTK_REAL const dy = cctkGH->cctk_delta_space[1] / cctkGH->cctk_levfac[1];
CCTK_REAL const dz = cctkGH->cctk_delta_space[2] / cctkGH->cctk_levfac[2];
CCTK_REAL const idx = 1.0/dx;
CCTK_REAL const idy = 1.0/dy;
CCTK_REAL const idz = 1.0/dz;
for (int k=bmin[2]; k<bmax[2]; ++k) {
for (int j=bmin[1]; j<bmax[1]; ++j) {
for (int i=bmin[0]; i<bmax[0]; ++i) {
int const ind = CCTK_GFINDEX3D(cctkGH, i,j,k);
{
// The main part of the boundary condition assumes that we
// have an outgoing radial wave with some speed v0:
//
// var = var0 + u(r-v0*t)/r
//
// This implies the following differential equation:
//
// d_t var = - v^i d_i var - v0 (var - var0) / r
//
// where vi = v0 xi/r
// Find local wave speeds
CCTK_REAL const rp = r[ind];
CCTK_REAL const rpi = 1.0/rp;
CCTK_REAL const vx = v0*x[ind]*rpi;
CCTK_REAL const vy = v0*y[ind]*rpi;
CCTK_REAL const vz = v0*z[ind]*rpi;
int const svi = i==0 ? +1 : i==ni-1 ? -1 : vx>0 ? +1 : -1;
int const svj = j==0 ? +1 : j==nj-1 ? -1 : vy>0 ? +1 : -1;
int const svk = k==0 ? +1 : k==nk-1 ? -1 : vz>0 ? +1 : -1;
// Find x derivative
CCTK_REAL derivx;
if (i>0 and i<ni-1) {
derivx = 0.5*(var[ind+di]-var[ind-di])*idx;
} else {
derivx = svi*0.5*(3*var[ind] - 4*var[ind-svi*di]
+ var[ind-2*svi*di])*idx;
}
// Find y derivative
CCTK_REAL derivy;
if (j>0 and j<nj-1) {
derivy = 0.5*(var[ind+dj]-var[ind-dj])*idy;
} else {
derivy = svj*0.5*(3*var[ind] - 4*var[ind-svj*dj]
+ var[ind-2*svj*dj])*idy;
}
// Find z derivative
CCTK_REAL derivz;
if (k>0 and k<nk-1) {
derivz = 0.5*(var[ind+dk]-var[ind-dk])*idz;
} else {
derivz = svk*0.5*(3*var[ind] - 4*var[ind-svk*dk]
+ var[ind-2*svk*dk])*idz;
}
// Calculate source term
rhs[ind] =
- vx*derivx - vy*derivy - vz*derivz - v0*(var[ind] - var0)*rpi;
}
if (radpower >= 0) {
// *****************************************
// *** EXTRAPOLATION OF MISSING PART ***
// *****************************************
//
// Here we try to extrapolate for the part of the boundary
// that does not behave as a pure wave (i.e. Coulomb type
// terms caused by infall of the coordinate lines).
//
// This we do by comparing the source term one grid point
// away from the boundary (which we already have), to what
// we would have obtained if we had used the boundary
// condition there. The difference gives us an idea of the
// missing part and we extrapolate that to the boundary
// assuming a power-law decay.
int const ip = i-si;
int const jp = j-sj;
int const kp = k-sk;
// Find local wave speeds
int const indp = CCTK_GFINDEX3D(cctkGH, ip,jp,kp);
CCTK_REAL const rp = r[indp];
CCTK_REAL const rpi = 1.0/rp;
CCTK_REAL const vx = v0*x[indp]*rpi;
CCTK_REAL const vy = v0*y[indp]*rpi;
CCTK_REAL const vz = v0*z[indp]*rpi;
int const svi = i==0 ? +1 : i==ni-1 ? -1 : vx>0 ? +1 : -1;
int const svj = j==0 ? +1 : j==nj-1 ? -1 : vy>0 ? +1 : -1;
int const svk = k==0 ? +1 : k==nk-1 ? -1 : vz>0 ? +1 : -1;
// Find x derivative
CCTK_REAL derivx;
if (ip>0 and ip<ni-1) {
derivx = 0.5*(var[indp+di]-var[indp-di])*idx;
} else {
derivx = svi*0.5*(3*var[indp] - 4*var[indp-svi*di]
+ var[indp-2*svi*di])*idx;
}
// Find y derivative
CCTK_REAL derivy;
if (jp>0 and jp<nj-1) {
derivy = 0.5*(var[indp+dj]-var[indp-dj])*idy;
} else {
derivy = svj*0.5*(3*var[indp] - 4*var[indp-svj*dj]
+ var[indp-2*svj*dj])*idy;
}
// Find z derivative
CCTK_REAL derivz;
if (kp>0 and kp<nk-1) {
derivz = 0.5*(var[indp+dk]-var[indp-dk])*idz;
} else {
derivz = svk*0.5*(3*var[indp] - 4*var[indp-svk*dk]
+ var[indp-2*svk*dk])*idz;
}
// Find difference in sources
CCTK_REAL const aux =
rhs[indp] +
vx*derivx + vy*derivy + vz*derivz + v0*(var[indp] - var0)*rpi;
// Extrapolate difference and add it to source in boundary
rhs[ind] += aux*pow(rp/r[ind],radpower);
} // if radpower>=0
} // for i j k
}
}
}
// Adapted from Kranc's KrancNumericalTools/GenericFD's file
// GenericFD.c
static
void newrad_loop (cGH const* restrict const cctkGH,
CCTK_REAL const* restrict const var,
CCTK_REAL * restrict const rhs,
CCTK_REAL const* restrict const x,
CCTK_REAL const* restrict const y,
CCTK_REAL const* restrict const z,
CCTK_REAL const* restrict const r,
CCTK_REAL const& var0,
CCTK_REAL const& v0,
CCTK_REAL const& radpower,
int const width)
{
int imin[3], imax[3], is_symbnd[6], is_physbnd[6], is_ipbnd[6];
GenericFD_GetBoundaryInfo
(cctkGH, cctkGH->cctk_lsh, cctkGH->cctk_lssh, cctkGH->cctk_bbox,
cctkGH->cctk_nghostzones,
imin, imax, is_symbnd, is_physbnd, is_ipbnd);
// Loop over all faces
for (int dir2=-1; dir2<=+1; ++dir2) {
for (int dir1=-1; dir1<=+1; ++dir1) {
for (int dir0=-1; dir0<=+1; ++dir0) {
int const dir[3] = { dir0, dir1, dir2 };
// one of tahe faces is a boundary
bool have_bnd = false;
// all boundary faces are physical boundaries
bool all_physbnd = true;
int bmin[3], bmax[3];
for (int d=0; d<3; ++d) {
switch (dir[d]) {
case -1:
bmin[d] = 0;
bmax[d] = imin[d];
have_bnd = true;
all_physbnd = all_physbnd and is_physbnd[2*d+0];
break;
case 0:
bmin[d] = imin[d];
bmax[d] = imax[d];
break;
case +1:
bmin[d] = imax[d];
bmax[d] = cctkGH->cctk_lssh[CCTK_LSSH_IDX(0,d)];
have_bnd = true;
all_physbnd = all_physbnd and is_physbnd[2*d+1];
break;
}
}
if (have_bnd and all_physbnd) {
newrad_kernel (cctkGH, bmin, bmax, dir,
var, rhs, x,y,z,r, var0, v0, radpower);
}
} // for dir0 dir1 dir2
}
}
}
extern "C"
CCTK_INT NewRad_Apply1 (CCTK_POINTER_TO_CONST const cctkGH_,
CCTK_REAL const* restrict const var,
CCTK_REAL * restrict const rhs,
CCTK_REAL const var0,
CCTK_REAL const v0,
CCTK_REAL const radpower,
CCTK_INT const width)
{
cGH const* restrict const cctkGH = static_cast<cGH const*> (cctkGH_);
if (not cctkGH) {
CCTK_WARN (CCTK_WARN_ABORT,
"cctkGH is NULL");
}
#if 0
CCTK_REAL const* restrict const var = CCTK_VarDataPtr (cctkGH, 0, varname);
if (not var) {
CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
"Cannot access variable \"%s\"", varname);
}
CCTK_REAL * restrict const rhs = CCTK_VarDataPtr (cctkGH, 0, rhsname);
if (not rhs) {
CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
"Cannot access RHS variable \"%s\"", rhsname);
}
#endif
if (not var) {
CCTK_WARN (CCTK_WARN_ABORT,
"Pointer to variable is NULL");
}
if (not rhs) {
CCTK_WARN (CCTK_WARN_ABORT,
"Pointer to RHS is NULL");
}
CCTK_REAL const* restrict const x =
static_cast<CCTK_REAL const*> (CCTK_VarDataPtr (cctkGH, 0, "grid::x"));
CCTK_REAL const* restrict const y =
static_cast<CCTK_REAL const*> (CCTK_VarDataPtr (cctkGH, 0, "grid::y"));
CCTK_REAL const* restrict const z =
static_cast<CCTK_REAL const*> (CCTK_VarDataPtr (cctkGH, 0, "grid::z"));
CCTK_REAL const* restrict const r =
static_cast<CCTK_REAL const*> (CCTK_VarDataPtr (cctkGH, 0, "grid::r"));
if (not x or not y or not z or not z) {
CCTK_WARN (CCTK_WARN_ABORT,
"Cannot access coordinate variables x, y, z, and r");
}
newrad_loop (cctkGH, var, rhs, x,y,z,r, var0, v0, radpower, width);
return 0;
}
|
