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
|
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "../../../CactusEinstein/Einstein/src/Einstein.h"
subroutine setupbrilldata2D(CCTK_FARGUMENTS)
C Author: Carsten Gundlach.
C
C Set up axisymmetric Brill data for elliptic solve. The elliptic
C equation we need to solve is:
C
C __ 2 2
C \/ psi = psi ( d q + d q ) / 4
C f z rho
C
C where:
C
C __
C \/ = Flat space Laplacian.
C f
implicit none
DECLARE_CCTK_FARGUMENTS
DECLARE_CCTK_PARAMETERS
integer i,j,k
integer nx,ny,nz
CCTK_REAL x1,y1,z1,rho1
CCTK_REAL brillq,eps
CCTK_REAL zero,one
external brillq
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.
eps = brill_eps
C Initialize psi.
brillpsi = one
C Set up flat metric for the elliptic solve.
C [Delta(g) + Mlinear] psi = Nsource.
psi = one
psix = zero
psiy = zero
psiz = zero
gxx = one
gyy = one
gzz = one
gxy = zero
gxz = zero
gyz = zero
C Set up coefficient Mlinear = - (1/4) (q,zz + q,rhorho).
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)
rho1 = dsqrt(x1**2 + y1**2)
C Centered derivatives. Note that here we may be calling brillq
C with a small negative rho, but that should be ok as long as
C brillq is even in rho - physically it must be, or the data
C will not be regular on the axis.
Mlinear(i,j,k) = - 0.25d0
$ *(brillq(rho1,z1+eps,zero)
$ + brillq(rho1,z1-eps,zero)
$ + brillq(rho1+eps,z1,zero)
$ + brillq(rho1-eps,z1,zero)
$ - 4.d0*brillq(rho1,z1,zero))/eps**2
end do
end do
end do
C Set coefficient Nsource = 0.
Nsource = zero
return
end
|