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
|
/*@@
@file finishbrilldata.F
@date
@author Carsten Gundlach (Cactus 4, Miguel Alcubierre)
@desc
Reconstruct metric from conformal factor.
@enddesc
@version $Header$
@@*/
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
subroutine IDBrillData_Finish(CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_FUNCTIONS
integer i,j,k
integer nx,ny,nz
CCTK_REAL x1,y1,z1,rho1,rho2
CCTK_REAL phi,psi4,e2q,rhofudge
CCTK_REAL zero,one
CCTK_REAL brillq,phif
c Numbers.
zero = 0.0D0
one = 1.0D0
c Set up grid size.
nx = cctk_lsh(1)
ny = cctk_lsh(2)
nz = cctk_lsh(3)
c Parameters.
rhofudge = brill_rhofudge
c Replace flat metric left over from elliptic solve by
c Brill metric calculated from q and Psi.
do k=1,nz
do j=1,ny
do i=1,nx
x1 = x(i,j,k)
y1 = y(i,j,k)
z1 = z(i,j,k)
rho2 = x1*x1 + y1*y1
rho1 = dsqrt(rho2)
phi = phif(x1,y1)
e2q = dexp(2.d0*brillq(rho1,z1,phi))
c Fudge division by rho^2 on axis. (Physically, y^/rho^2,
c x^2/rho^2 and xy/rho^2 are of course regular.)
c Transform Brills form of the physical metric to Cartesian
c coordinates via
c
c e^2q (drho^2 + dz^2) + rho^2 dphi^2 =
c e^2q (dx^2 + dy^2 + dz^2) + (1-e^2q) (xdy-ydx)^2/rho^2
c
c The individual coefficients can be read off as
if (rho1.gt.rhofudge) then
gxx(i,j,k) = (e2q + (one - e2q)*y1*y1/rho2)
gyy(i,j,k) = (e2q + (one - e2q)*x1*x1/rho2)
gzz(i,j,k) = e2q
gxy(i,j,k) = - (one - e2q)*x1*y1/rho2
else
c This fudge assumes that q = O(rho^2) near the axis. Which
c it should be, or the data will be singular.
gxx(i,j,k) = 0.0d0
gyy(i,j,k) = 0.0d0
gzz(i,j,k) = 0.0d0
gxy(i,j,k) = 0.0d0
end if
end do
end do
end do
if (CCTK_EQUALS(metric_type,"static conformal")) then
conformal_state = 1
do k=1,nz
do j=1,ny
do i=1,nx
psi(i,j,k) = brillpsi(i,j,k)
end do
end do
end do
else
conformal_state = 0
do k=1,nz
do j=1,ny
do i=1,nx
psi4 = brillpsi(i,j,k)**4
gxx(i,j,k) = psi4*gxx(i,j,k)
gyy(i,j,k) = psi4*gyy(i,j,k)
gzz(i,j,k) = psi4*gzz(i,j,k)
gxy(i,j,k) = psi4*gxy(i,j,k)
end do
end do
end do
end if
c In any case,
gxz = zero
gyz = zero
c Vanishing extrinsic curvature completes the Cauchy data.
kxx = zero
kyy = zero
kzz = zero
kxy = zero
kxz = zero
kyz = zero
return
end
|
