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
|
c The metric given here is not a solution of
c Einsteins equations! It is nevertheless
c useful for tests since it has a particularly
c nice geometry. It is a static, spherically
c symmetric metric with no shift.
c
c In spherical coordinates, the metric has the
c form:
c
c 2 2 2 2
c ds = dr + R(r) d Omega
c
c Clearly, r measures radial proper distance, and R(r)
c is the areal (Schwarzschild) radius.
c
c I choose a form of R(r) such that:
c
c R --> r r<<1, r>>1
c
c So close in, and far away we have a flat metric.
c In the middle region, I take R to be smaller than
c r, but still larger than zero. This deficit in
c areal radius produces the geometry of a "bag of gold".
c
c The size of the deviation from a flat geometry
c is controled by the parameter "bowl__strength".
c If this parameter is 0, we are in flat space.
c The width of the curved region is controled by
c the paramter "bowl__sigma", and the place where the
c curvature becomes significant (the center of the
c deformation) is controled by "bowl__center".
c
c The specific form of the function R(r) is:
c
c R(r) = ( r - a f(r) )
c
c and the form of thr function f(r) depends on the
c parameter bowl__shape:
c 2 2 2
c bowl__shape = "Gaussian" f(r) = exp(-(r-c) / s )
c
c bowl__shape = "Fermi" f(r) = 1 / ( 1 + exp(-s(r-c)) )
c
c There are three extra paramters
c (bowl__x_scale,bowl__y_scale,bowl__z_scale) that set the
c scales for the (x,y,z) axis respectively. Their default
c value are all 1. These parameters are useful to hide the
c spherical symmetry of the metric.
C
C Author: unknown
C Copyright/License: unknown
c
c $Header$
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Functions.h"
subroutine Exact__bowl(
$ x, y, z, t,
$ gdtt, gdtx, gdty, gdtz,
$ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
$ gutt, gutx, guty, gutz,
$ guxx, guyy, guzz, guxy, guyz, guzx,
$ psi, Tmunu_flag)
implicit none
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_FUNCTIONS
c input arguments
CCTK_REAL x, y, z, t
c output arguments
CCTK_REAL gdtt, gdtx, gdty, gdtz,
$ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
$ gutt, gutx, guty, gutz,
$ guxx, guyy, guzz, guxy, guyz, guzx
CCTK_REAL psi
LOGICAL Tmunu_flag
c local static variables
logical firstcall,evolve
CCTK_REAL a,c,s
CCTK_REAL dx,dy,dz
CCTK_REAL t0,st
data firstcall /.true./
save firstcall,evolve,type,a,c,s,dx,dy,dz,t0,st
c local variables
integer type
character*100 warn_buffer
CCTK_REAL r,r2,rr2
CCTK_REAL xr,yr,zr,xr2,yr2,zr2
CCTK_REAL fac,det
CCTK_REAL tfac
c constants
CCTK_REAL zero,one,two
parameter (zero=0.0d0, one=1.0d0, two=2.0d0)
C This is a vacuum spacetime with no cosmological constant
Tmunu_flag = .false.
c Get parameters of the metric.
if (firstcall) then
a = bowl__strength
c = bowl__center
s = bowl__sigma
dx = bowl__x_scale
dy = bowl__y_scale
dz = bowl__z_scale
if (CCTK_Equals(bowl__shape,"Gaussian").ne.0) then
type = 1
else if (CCTK_Equals(bowl__shape,"Fermi").ne.0) then
type = 2
else
write (warn_buffer, '(a,a,a)')
$ 'Unknown bowl__shape = "', bowl__shape, '"'
call CCTK_WARN(0, warn_buffer)
end if
if (bowl__evolve.eq.1) then
evolve = .true.
t0 = bowl__t0
st = bowl__sigma_t
else
evolve = .false.
end if
firstcall = .false.
end if
c Multiply the bowl strength "a" with a time evolution factor.
c The time evolution factor is taken to be a Fermi step centered
c in "t0" and with a width "st". The size of this step is always
c 1 so that far in the past we will always have flat space, and
c far in the future we will have the static bowl.
if (evolve) then
tfac = one/(one + exp(-st*(t-t0)))
else
tfac = one
end if
a = a*tfac
c Find {r2,r}.
r2 = (x/dx)**2 + (y/dy)**2 + (z/dz)**2
r = sqrt(r2)
c Find the form function rr2
c
c 2 2 2
c rr2 = (r - a f) / r = (1 - a f / r)
if (type.eq.1) then
c Gaussian bowl:
c 2 2 2
c rr2 = [ 1 - a exp(-(r-c) /s ) / r ]
c
c Notice that this really does not go to 1 at the
c origin. To fix this, I multiply the gaussian
c with the factor:
c
c fac = 1 - sech(4r)
c
c This goes smoothly to 0 at the origin, and climbs
c fast to a limiting value of 1 (at r=1 it is already
c equal to 0.96).
fac = one - two/(exp(4.0d0*r) + exp(-4.0d0*r))
rr2 = (one - a*fac*exp(-((r-c)/s)**2)/r)**2
else if (type.eq.2) then
c Fermi bowl:
c 2
c rr2 = [ 1 - 1 / ( 1 + exp(-s(r-c)) ) / r ]
c
c Again, this doesnt really go to 1 at the origin, so
c I use the same trick as above.
fac = one - two/(exp(4.0d0*r) + exp(-4.0d0*r))
rr2 = (one - a*fac/(one + exp(-s*(r-c)))/r)**2
else
write (warn_buffer, '(a,i8)')
$ 'Unknown type = ', type
call CCTK_WARN(0, warn_buffer)
end if
c Give metric components.
gdtt = - one
gdtx = zero
gdty = zero
gdtz = zero
if (r.ne.0) then
xr = (x/dx)/r
yr = (y/dy)/r
zr = (z/dz)/r
xr2 = xr**2
yr2 = yr**2
zr2 = zr**2
gdxx = (xr2 + rr2*(yr2 + zr2))/dx**2
gdyy = (yr2 + rr2*(xr2 + zr2))/dy**2
gdzz = (zr2 + rr2*(xr2 + yr2))/dz**2
gdxy = xr*yr*(one - rr2)/(dx*dy)
gdyz = yr*zr*(one - rr2)/(dy*dz)
gdzx = xr*zr*(one - rr2)/(dx*dz)
else
gdxx = one
gdyy = one
gdzz = one
gdxy = zero
gdyz = zero
gdzx = zero
end if
c Find inverse metric.
gutt = - one
gutx = zero
guty = zero
gutz = zero
det = gdxx*gdyy*gdzz + two*gdxy*gdzx*gdyz
. - gdxx*gdyz**2 - gdyy*gdzx**2 - gdzz*gdxy**2
guxx = (gdyy*gdzz - gdyz**2)/det
guyy = (gdxx*gdzz - gdzx**2)/det
guzz = (gdxx*gdyy - gdxy**2)/det
guxy = (gdzx*gdyz - gdzz*gdxy)/det
guyz = (gdxy*gdzx - gdxx*gdyz)/det
guzx = (gdxy*gdyz - gdyy*gdzx)/det
return
end
|
