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
|
/*@@
@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
CCTK_REAL x1,y1,z1,rho1,rho2
CCTK_REAL phi,psi4,e2q
CCTK_REAL zero,one
CCTK_REAL brillq,phif
c Numbers.
zero = 0.0D0
one = 1.0D0
c Replace flat metric left over from elliptic solve by
c Brill metric calculated from q and Psi.
do k=1,cctk_lsh(3)
do j=1,cctk_lsh(2)
do i=1,cctk_lsh(1)
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) = zero
gyy(i,j,k) = zero
gzz(i,j,k) = zero
gxy(i,j,k) = zero
end if
gxz(i,j,k) = zero
gyz(i,j,k) = zero
kxx(i,j,k) = zero
kyy(i,j,k) = zero
kzz(i,j,k) = zero
kxy(i,j,k) = zero
kxz(i,j,k) = zero
kyz(i,j,k) = zero
end do
end do
end do
if (CCTK_EQUALS(metric_type,"static conformal")) then
conformal_state = 1
do k=1,cctk_lsh(3)
do j=1,cctk_lsh(2)
do i=1,cctk_lsh(1)
psi(i,j,k) = brillpsi(i,j,k)
end do
end do
end do
else
conformal_state = 0
do k=1,cctk_lsh(3)
do j=1,cctk_lsh(2)
do i=1,cctk_lsh(1)
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
return
end
|
