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
|
/*@@
@file AHFinder_mask.F
@date November 1998
@author Miguel Alcubierre
@desc
Set up horizon mask for black hole excision.
@enddesc
@@*/
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
subroutine AHFinder_mask(CCTK_ARGUMENTS,rhor)
! Set up the horizon mask. All this routine does is
! set up the value of the grid function "ahmask" to zero
! for points inside the horizon.
!
! It also finds out the radius of the largest sphere that
! fits inside the horizon and it centered in (xc,yc,zc).
! This is useful for simple excision.
use AHFinder_dat
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_FUNCTIONS
integer i,j,k
integer reduce_handle
CCTK_REAL zero,one
CCTK_REAL rhor,rhortemp
CCTK_REAL buffer,dd,aux
CCTK_REAL xa,ya,za
! *******************
! *** NUMBERS ***
! *******************
zero = 0.0D0
one = 1.0D0
! *********************************
! *** FIND HORIZON FUNCTION ***
! *********************************
call AHFinder_fun(CCTK_ARGUMENTS)
! ***********************************
! *** FIND BUFFER ZONE FACTOR ***
! ***********************************
! The mask will always be set inside the horizon minus
! a safety buffer zone. To find out how big is the buffer
! zone I check the value of the parameter "ahf_maskshrink".
! This parameter says by which factor of the monopole term
! to shrink the mask with respect to the real horizon.
! Here I use the fact that the function "ahfgrid" is zero
! at the horizon, and is linear on the monopole term
! (which measures the overall scale). What I then do is
! set up the value of "buffer" to a negative number that
! roughly measures the desired coordinate distance from
! the mask to the horizon. The reason why it is negative
! is because the horizon function is negative inside the
! horizon, and I just want to compare its value with the
! value of "buffer" to decide when a given point is masked.
buffer = - (one-ahf_maskshrink)*c0(0)
! Now I make sure that there are at least "ahf_maskbuffer"
! grid points in the buffer zone. The final buffer zone will
! be whatever is largest between the number of grid points
! specified by "ahf_maskbuffer" and the shrink factor
! specified by "ahf_maskshrink".
dd = max(dx,dy,dz)
aux = ahf_maskbuffer*dd
if (abs(buffer).lt.aux) buffer = - aux
! *****************************************
! *** FIND RADIUS OF LARGEST SPHERE ***
! *****************************************
! Find the radius of the largest sphere that
! fits inside the horizon.
rhor = 1.0D10
do k=1,nz
do j=1,ny
do i=1,nx
if (ahfgrid(i,j,k).ge.buffer) then
aux = sqrt((x(i,j,k)-xc)**2 + (y(i,j,k)-yc)**2
. + (z(i,j,k)-zc)**2)
if (rhor.gt.aux) rhor = aux
end if
end do
end do
end do
! Now find the minimum of rhor across processors.
reduce_handle = -1
call CCTK_ReductionArrayHandle(reduce_handle,"minimum")
if (reduce_handle.lt.0) then
call CCTK_WARN(1,"Cannot get handle for minimum reduction ! Forgot to activate an implementation providing reduction operators ??")
end if
call CCTK_ReduceLocalScalar(ierr ,cctkGH,-1,reduce_handle,
. rhor,rhortemp,CCTK_VARIABLE_REAL)
rhor = rhortemp
if (ierr.ne.0) then
call CCTK_WARN(1,"AHFinder_mask:Reduction of rhor failed!")
end if
! ***********************
! *** LEGO SPHERE ***
! ***********************
! If we want a lego shere mask life is simple.
! Just set the mask inside of horizon (minus a
! buffer zone) equal to 0.
if (CCTK_EQUALS(ahf_masktype,'lego')) then
do k=1,nz
do j=1,ny
do i=1,nx
if (ahfgrid(i,j,k).lt.buffer) then
ahmask(i,j,k) = zero
end if
end do
end do
end do
! ************************
! *** CUBICAL MASK ***
! ************************
! If we want a cubical mask, set the mask based on the value
! of rhor. Notice that if we want the corners of the cube to
! touch the sphere of radius rhor, we must make the side of
! the cube equal to rhor/sqrt(3), which is about 0.57*rhor.
else if (CCTK_EQUALS(ahf_masktype,'cube')) then
do k=1,nz
do j=1,ny
do i=1,nx
aux = 0.57D0*rhor
xa = abs(x(i,j,k)-xc)
ya = abs(y(i,j,k)-yc)
za = abs(z(i,j,k)-zc)
if ((xa.lt.aux).and.(ya.lt.aux).and.(za.lt.aux)) then
ahmask(i,j,k) = zero
end if
end do
end do
end do
! **************************
! *** POLYHEDRA MASK ***
! **************************
! Here we want to excise a polyhedra. But not just any
! polyhedra, a polyhedra defined by taking a cube and
! cutting away the corners all the way up to the middle
! of the edges. This results in a figure made up of
! 6 squares and 8 equilateral triangles.
! Notice that for such a figure, if we want the corners
! to touch a sphere of radius rhor, then the sides have
! to be rhor/sqrt(2) which is roughly 0.7*rhor.
else if (CCTK_EQUALS(ahf_masktype,'poly')) then
do k=1,nz
do j=1,ny
do i=1,nx
aux = 0.7D0*rhor
xa = abs(x(i,j,k)-xc)
ya = abs(y(i,j,k)-yc)
za = abs(z(i,j,k)-zc)
if ((xa.lt.aux).and.(ya.lt.aux).and.(za.lt.aux).and.
. ((xa+ya+za).lt.2.0D0*aux)) then
ahmask(i,j,k) = zero
end if
end do
end do
end do
end if
! *****************************************************
! *** APPLY SYMMETRY BOUNDARIES AND SYNCHRONIZE ***
! *****************************************************
call CCTK_SyncGroup(ierr,cctkGH,"ahfinder::ahfmask")
call CartSymGN(ierr,cctkGH,"ahfinder::ahfmask")
! If desired, copy ahf_mask to spacemask mask (emask).
if (use_mask.eq.1) then
emask = ahmask
end if
! ***************
! *** END ***
! ***************
end subroutine AHFinder_mask
|
