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
|
/*@@
@file setupbrilldata2D.F
@date
@author Carsten Gundlach (Cactus 4, Miguel Alcubierre)
@desc
Set up axisymmetric Brill data for elliptic solve.
@enddesc
@version $Header$
@@*/
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
subroutine setupbrilldata2D(CCTK_ARGUMENTS)
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 = 0
c f z rho
c
c where:
c
c __
c \/ = Flat space Laplacian.
c f
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
integer i,j,k
integer ierr
CCTK_REAL x1,y1,z1,rho1
CCTK_REAL brillq,eps
CCTK_REAL zp,zm,rhop,rhom
CCTK_REAL zero,one
external brillq
c Numbers.
zero = 0.0D0
one = 1.0D0
c Epsilon for finite differencing.
eps = cctk_delta_space(1)
c Set up coefficient arrays for elliptic solve.
c Notice that the Cactus conventions are:
c __
c \/ psi + Mlinear*psi + Nsource = 0
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)
rho1 = dsqrt(x1*x1 + y1*y1)
brillpsi(i,j,k) = one
gxx(i,j,k) = one
gyy(i,j,k) = one
gzz(i,j,k) = one
gxy(i,j,k) = zero
gxz(i,j,k) = zero
gyz(i,j,k) = zero
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.
zp = z1 + eps
zm = z1 - eps
rhop = rho1 + eps
rhom = rho1 - eps
brillMlinear(i,j,k) = 0.25D0
. *(brillq(rho1,zp,zero)
. + brillq(rho1,zm,zero)
. + brillq(rhop,z1,zero)
. + brillq(rhom,z1,zero)
. - 4.0D0*brillq(rho1,z1,zero))/eps**2
brillNsource(i,j,k) = zero
end do
end do
end do
c Synchronization and boundaries.
c call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic')
call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic')
return
end
|
