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
|
/*@@
@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 finishbrilldata(CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
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**2 + y1**2
rho1 = dsqrt(rho2)
phi = phif(x1,y1)
psi4 = brillpsi(i,j,k)**4
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) = psi4*(e2q + (one - e2q)*y1**2/rho2)
gyy(i,j,k) = psi4*(e2q + (one - e2q)*x1**2/rho2)
gzz(i,j,k) = psi4*e2q
gxy(i,j,k) = - psi4*(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) = psi4
gyy(i,j,k) = psi4
gzz(i,j,k) = psi4
gxy(i,j,k) = zero
end if
end do
end do
end do
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
|
