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
|
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* Using Cactus infrastructure */
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
static const char *rcsid = "$Header$";
CCTK_FILEVERSION(CactusWave_WaveBinarySource_WaveBinary_c)
int IndexCeilC(cGH *GH, CCTK_REAL coord_value, int d);
int IndexFloorC(cGH *GH, CCTK_REAL coord_value, int d);
void WaveBinaryC(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
int i,j,k,d,index;
/* the start/end points in local index space of the binary charge */
int lowerloc[3];
int upperloc[3];
/* lower/upper physical coord. of one binary charge */
CCTK_REAL lbin[3],ubin[3];
CCTK_REAL minc[3],maxc[3];
/* lower/upper grid points, that we own (no ghostzones) */
int mingp[3], maxgp[3];
/* some flags */
int onthisproc, xproc, sign;
static int firstcall;
CCTK_REAL xs,ys,zs,rad;
CCTK_REAL zmin,zmax,charge_factor;
/* Because binary_charge and binary_size are a steerable parameters now,
charge_factor needs to be recomputed at every iteration. */
charge_factor = 3.0*binary_charge/
(4.0*3.1415*binary_size*binary_size*binary_size);
/* Initialize the range arrays */
for (d=0;d<3;d++)
{
mingp[d] = 0;
maxgp[d] = cctk_lsh[d]-1;
}
/* The physical minimum/maximum coordinates of our grid chunk;
will use those to check if the charges are on our proc. */
index = CCTK_GFINDEX3D(cctkGH,0,0,0);
minc[0] = x[index];
minc[1] = y[index];
minc[2] = z[index];
index = CCTK_GFINDEX3D(cctkGH,cctk_lsh[0]-1,cctk_lsh[1]-1,cctk_lsh[2]-1);
maxc[0] = x[index];
maxc[1] = y[index];
maxc[2] = z[index];
/* we have two charges opposite, origin is the center
sign mulitplies the charge xy-positions by +1 and -1 */
for (sign=1; sign>=-1; sign-=2)
{
/* default flag: the charge is on this proc. */
onthisproc = 0;
/* calculate the center of a single binary source */
CCTK_CoordRange(cctkGH, &zmin, &zmax, -1, "z", "cart3d");
xs = sign * binary_radius * cos(binary_omega*(cctk_time-cctk_delta_time));
ys = sign * binary_radius * sin(binary_omega*(cctk_time-cctk_delta_time));
zs = (zmax-zmin)/2 + zmin;
/* shortcuts for the physical boundaries of
a single charge */
lbin[0] = xs - binary_size;
lbin[1] = ys - binary_size;
lbin[2] = zs - binary_size;
ubin[0] = xs + binary_size;
ubin[1] = ys + binary_size;
ubin[2] = zs + binary_size;
/* loop over all three dimensions */
for (d=0;d<3;d++)
{
/* check where our binaries are:
- the can be on our grid chunk completly
- we have a lower part of the source on our grid chunk
- we have a upper part
- the source covers a grid chunk completetly
Get the floor-index of the lower bound and
ceiling-index of the upper bound.
Increase onthisproc : if equal to three, we have the source
*/
xproc = 0;
lowerloc[d]=mingp[d];
upperloc[d]=maxgp[d];
if ((lbin[d]<=minc[d]) && (ubin[d]>=maxc[d]))
{
xproc = 1;
}
else
{
if ((minc[d]<=lbin[d]) && (lbin[d]<=maxc[d]))
{
lowerloc[d]=IndexFloorC(cctkGH, lbin[d],d);
xproc = 1;
}
if ((minc[d]<=ubin[d]) && (ubin[d]<=maxc[d]))
{
upperloc[d]=IndexCeilC (cctkGH, ubin[d],d);
xproc = 1;
}
}
onthisproc+=xproc;
}
/* Debugging debugging */
if (((firstcall==1)&&(CCTK_EQUALS(binary_verbose,"yes"))) ||
(CCTK_EQUALS(binary_verbose,"debug")))
{
printf("Charge: %f \n",charge_factor);
printf("On this proc: %d \n",onthisproc);
printf("cctk_lsh: %d %d %d\n",cctk_lsh[0],cctk_lsh[0],cctk_lsh[0]);
printf("Charge center: %f %f %f \n", xs,ys,zs);
printf("Charge extension: \n");
printf(" x-extension: %f -> %f \n",lbin[0],ubin[0]);
printf(" y-extension: %f -> %f \n",lbin[1],ubin[1]);
printf(" z-extension: %f -> %f \n",lbin[2],ubin[2]);
printf("Charge local index range: \n");
printf(" x-index %d %d \n", lowerloc[0],upperloc[0]);
printf(" y-index %d %d \n", lowerloc[1],upperloc[1]);
printf(" z-index %d %d \n", lowerloc[2],upperloc[2]);
index = CCTK_GFINDEX3D(cctkGH, lowerloc[0],lowerloc[1],lowerloc[2]);
printf("Charge corner LOW: [%f %f %f] \n",x[index],y[index],z[index]);
index = CCTK_GFINDEX3D(cctkGH, upperloc[0],upperloc[1],upperloc[2]);
printf("Charge corner UP : [%f %f %f] \n",x[index],y[index],z[index]);
printf(" MINC/MAXC [%f %f %f] [%f %f %f] \n\n\n",
minc[0],minc[1],minc[2],
maxc[0],maxc[1],maxc[2]);
}
/* we have the source, no loop over the index boudary */
if (onthisproc==3)
{
for (k=lowerloc[2];k<=upperloc[2];k++)
{
for (j=lowerloc[1];j<=upperloc[1];j++)
{
for (i=lowerloc[0];i<=upperloc[0];i++)
{
index = CCTK_GFINDEX3D(cctkGH,i,j,k);
rad =
((x[index]-xs)*(x[index]-xs)) +
((y[index]-ys)*(y[index]-ys)) +
((z[index]-zs)*(z[index]-zs));
if (rad<binary_size*binary_size)
{
phi[index] = phi[index] + charge_factor;
}
}
}
}
}
/* end of the sign loop */
}
/* reset firstcall to zero */
firstcall = 0;
/* Note, that we do not need to sync anything, since each grid
patch has filled out its ghostzones. */
return;
}
|
