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
|
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "../../../CactusEinstein/Einstein/src/Einstein.h"
subroutine finishbrilldata(CCTK_FARGUMENTS)
C Author: Carsten Gundlach.
C
C Set up 3-metric, extrinsic curvature and BM variables
C once the conformal factor has been found.
implicit none
DECLARE_CCTK_FARGUMENTS
DECLARE_CCTK_PARAMETERS
integer i,j,k
integer nx,ny,nz
integer CCTK_Equals
CCTK_REAL x1,y1,z1,rho1,rho2
CCTK_REAL brillq,psi4,e2q,rhofudge
CCTK_REAL phi,phif
CCTK_REAL zero
external brillq
C Numbers.
zero = 0.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
gzz(i,j,k) = psi4*e2q
gxx(i,j,k) = psi4*(e2q + (1.d0 - e2q)*y1**2/rho2)
gyy(i,j,k) = psi4*(e2q + (1.d0 - e2q)*x1**2/rho2)
gxy(i,j,k) = - psi4*(1.d0 - 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.
gzz(i,j,k) = psi4
gxx(i,j,k) = psi4
gyy(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
|
