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
|
c Using Cactus infrastructure
#include "cctk.h"
c Using Cactus parameters
#include "cctk_Parameters.h"
c Using Cactus arguments lists
#include "cctk_Arguments.h"
/*@@
@routine WaveBinarySlow
@date Fri Jan 7 09:10:28 2000
@author Gerd Lanfermann
@desc
Provides the sources for rotating binary charges.
It does this in a straight forward, slow approach,
by looping over all grid points.
@enddesc
@calls CCTK_CoordRange
@calledby
@@*/
subroutine WaveBinarySlow(CCTK_FARGUMENTS)
implicit none
c Declare variables in argument list
DECLARE_CCTK_FARGUMENTS
c Declare parameters
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_FUNCTIONS
INTEGER i,j,k,ierr
INTEGER istart, jstart, kstart, iend, jend, kend
CCTK_REAL dx,dy,dz,dt
CCTK_REAL xs,ys,zs,rad
CCTK_REAL zmin,zmax,charge_factor
c Set up shorthands
c -----------------
dx = CCTK_DELTA_SPACE(1)
dy = CCTK_DELTA_SPACE(2)
dz = CCTK_DELTA_SPACE(3)
dt = CCTK_DELTA_TIME
charge_factor = 3.0d0*binary_charge/
$ (4.0d0*3.1415*binary_size*binary_size*binary_size)
call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,"z");
xs = binary_radius * cos(binary_omega*(cctk_time-dt))
ys = binary_radius * sin(binary_omega*(cctk_time-dt))
zs = (zmax-zmin)/2 +zmin
istart = 2
jstart = 2
kstart = 2
iend = cctk_lsh(1)-1
jend = cctk_lsh(2)-1
kend = cctk_lsh(3)-1
do k = kstart, kend
do j = jstart, jend
do i = istart, iend
rad =((x(i,j,k)-xs)**2)+
$ ((y(i,j,k)-ys)**2)+
$ ((z(i,j,k)-zs)**2)
if (rad.le.(binary_size**2)) then
phi(i,j,k) =phi(i,j,k)+ charge_factor
end if
rad =((x(i,j,k)+xs)**2)+
$ ((y(i,j,k)+ys)**2)+
$ ((z(i,j,k)-zs)**2)
if (rad.le.(binary_size**2)) then
phi(i,j,k) =phi(i,j,k)+ charge_factor
end if
end do
end do
end do
return
end
/*@@
@routine WaveBinary
@date Fri Jan 7 09:12:02 2000
@author Gerd Lanfermann
@desc
Provides the sources for rotating binary charges.
It does this in a more intelligent approach
by looping over the only gridpoints that are covered by
charge.
@enddesc
@calls CCTK_CoordRange IndexFloor IndexCeil IndexFloor IndexCeil IndexFloor IndexCeil
@calledby
@history
@@*/
subroutine WaveBinary(CCTK_FARGUMENTS)
implicit none
c Declare variables in argument list
DECLARE_CCTK_FARGUMENTS
c Declare parameters
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_FUNCTIONS
INTEGER i,j,k,ierr,d,f
c the start/end points in local index space of the binary charge
INTEGER lowerloc(3)
INTEGER upperloc(3)
c lower/upper grid points, that we own (no ghostzones)
INTEGER mingp(3),maxgp(3)
c some flags
INTEGER nothere,sign
CCTK_REAL dx,dy,dz,dt
CCTK_REAL xs,ys,zs,rad
CCTK_REAL zmin,zmax,charge_factor
c localgridstart/end: global physical coordinates where the
c patch of grid starts/ends
CCTK_REAL locgridstart(3),locgridend(3)
c Things to do on first iteration
INTEGER firstcall
DATA firstcall /1/
SAVE firstcall,charge_factor
if (firstcall.eq.1) then
charge_factor = 3.0d0*binary_charge/
$ (4.0d0*3.1415*binary_size*binary_size*binary_size)
end if
c Initialize the range arrays
do d=1,3
if (cctk_bbox(2*d-1).eq.1) then
mingp(d) = 1
else
mingp(d) = cctk_nghostzones(d)
end if
if (cctk_bbox(2*d).eq.1) then
maxgp(d) = cctk_lsh(d)
else
maxgp(d) = cctk_lsh(d)-cctk_nghostzones(d)
end if
end do
c we have two charges opposite, origin is the center
c sign mulitplies the charge xy-positions by +1 and -1
do sign=1,-1,-2
c default flag: the charge is not here
nothere = 1
c if we need to update phiold on first call
c calculate the position of the binary sources
call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,"z");
xs = sign*binary_radius * cos(binary_omega*cctk_time)
ys = sign*binary_radius * sin(binary_omega*cctk_time)
zs = (zmax-zmin)/2 +zmin
c get the local indices for the extension of the binary source
call IndexFloor(cctkGH, xs-binary_size, lowerloc(1), 1)
call IndexCeil (cctkGH, xs+binary_size, upperloc(1), 1)
call IndexFloor(cctkGH, ys-binary_size, lowerloc(2), 2)
call IndexCeil (cctkGH, ys+binary_size, upperloc(2), 2)
call IndexFloor(cctkGH, zs-binary_size, lowerloc(3), 3)
call IndexCeil (cctkGH, zs+binary_size, upperloc(3), 3)
do d=1,3
c Find out if the charge border can be found on our grid patch
c (return value in lowerloc >= 0)
if ((lowerloc(d).ge.0).or.(upperloc(d).ge.0)) then
nothere = 0
if ((lowerloc(d).ge.0).and.(upperloc(d).lt.0)) then
upperloc(d) = cctk_lsh(d)
end if
if ((upperloc(d).ge.0).and.(lowerloc(d).lt.0)) then
lowerloc(d) = 1
end if
else
c else check if the region of charge is spread across
c more than one processor: we have to make sure that the processor in
c between knows that it has to loop as well! While lower/upperloc
c have the values (-1), they should rather be 1 and cctk_lsh.
c We need to find out if the size of the binary is covering the
c physical size of the grid. To do that we look up the xyz value with
c the min/max index we actually own (no ghostzones).
if (d.eq.1) then
if ((xs-binary_size.lt.x(mingp(1),mingp(2),mingp(3))).and.
$ (xs+binary_size.gt.x(maxgp(1),maxgp(2),maxgp(3)))) then
lowerloc(d)=1
upperloc(d)=cctk_lsh(d)
nothere = 0
end if
else if (d.eq.2) then
if ((ys-binary_size.lt.y(mingp(1),mingp(2),mingp(3))).and.
$ (ys+binary_size.gt.y(maxgp(1),maxgp(2),maxgp(3)))) then
lowerloc(d)=1
upperloc(d)=cctk_lsh(d)
nothere = 0
end if
else if (d.eq.3) then
if ((zs-binary_size.lt.z(mingp(1),mingp(2),mingp(3))).and.
$ (zs+binary_size.gt.z(maxgp(1),maxgp(2),maxgp(3)))) then
lowerloc(d)=1
upperloc(d)=cctk_lsh(d)
nothere = 0
end if
end if
end if
end do
if (((firstcall.eq.1).and.(CCTK_EQUALS(binary_verbose,"yes")))
$ .or.(CCTK_EQUALS(binary_verbose,"debug"))) then
write (*,*)
write (*,*) "Charge: ",charge_factor
write (*,*) "Charge center: xs,ys,zs ", xs,ys,zs
write (*,*) "Charge extension: "
write (*,*) " x-extension: ",xs-binary_size, xs+binary_size
write (*,*) " y-extension: ",ys-binary_size, ys+binary_size
write (*,*) " z-extension: ",zs-binary_size, zs+binary_size
write (*,*) "Charge local index range"
write (*,*) " x-index", lowerloc(1),upperloc(1)
write (*,*) " y-index", lowerloc(2),upperloc(2)
write (*,*) " z-index", lowerloc(3),upperloc(3)
write (*,*)
end if
c Now loop over the grid points, that are covered by the charge
if (nothere.eq.0) then
do i=lowerloc(1),upperloc(1)
do j=lowerloc(2),upperloc(2)
do k=lowerloc(3),upperloc(3)
rad =((x(i,j,k)-xs)**2)+
$ ((y(i,j,k)-ys)**2)+
$ ((z(i,j,k)-zs)**2)
if (rad.le.(binary_size**2)) then
phi(i,j,k) = phi(i,j,k)+charge_factor
end if
end do
end do
end do
end if
c end of the sign-loop
end do
c we reset firstcall to zero
firstcall=0
c Note, that we do not need to sync anything, since each grid
c patch has filled out its ghostzones.
return
end
|
