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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
C $Header$
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
subroutine Exact__boundary(CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
integer i,j,k
integer nx,ny,nz
CCTK_REAL tplusone
CCTK_REAL
$ dxgxxjunk, dxgyyjunk, dxgzzjunk,
$ dxgxyjunk, dxgyzjunk, dxgxzjunk,
$ dygxxjunk, dygyyjunk, dygzzjunk,
$ dygxyjunk, dygyzjunk, dygxzjunk,
$ dzgxxjunk, dzgyyjunk, dzgzzjunk,
$ dzgxyjunk, dzgyzjunk, dzgxzjunk,
$ axjunk, ayjunk, azjunk,
$ bxxjunk, bxyjunk, bxzjunk,
$ byxjunk, byyjunk, byzjunk,
$ bzxjunk, bzyjunk, bzzjunk
CCTK_REAL
$ exact_psi,
$ exact_psix, exact_psiy, exact_psiz,
$ exact_psixx, exact_psiyy, exact_psizz,
$ exact_psixy, exact_psiyz, exact_psixz
C Grid parameters.
nx = cctk_lsh(1)
ny = cctk_lsh(2)
nz = cctk_lsh(3)
C Initialize the psi of exact
C (also to tell the models about the conformal_state)
if (conformal_state .ne. 0) then
exact_psi = 1.0D0
else
exact_psi = 0.0D0
end if
exact_psix = 0.0D0
exact_psiy = 0.0D0
exact_psiz = 0.0D0
exact_psixx = 0.0D0
exact_psiyy = 0.0D0
exact_psizz = 0.0D0
exact_psixy = 0.0D0
exact_psiyz = 0.0D0
exact_psixz = 0.0D0
C Set all initial data including dijk and vi on all points which
C are on the boundary of the domain if it really is the boundary
C of the complete grid. Treat all six sides of the grid cube this way.
c Set t = time + dt. This is necessary here because by the time
c we reach this point the geometry has been evolved one time step
c but the variable `time' still hasn't been updated.
tplusone = cctk_time + cctk_delta_time
C Note we also always set the lapse and shift at the boundaries at
C time t+1. This is to provide boundary conditions for testing
C elliptic gauge conditions. If they are not used, they will be
C overwritten by Exact__gauge.
#define EXACTDATAPOINT \
call Exact__Bona_Masso_data( \
decoded_exact_model, \
x(i,j,k), y(i,j,k), z(i,j,k), tplusone, \
gxx(i,j,k), gyy(i,j,k), gzz(i,j,k), \
gxy(i,j,k), gyz(i,j,k), gxz(i,j,k), \
kxx(i,j,k), kyy(i,j,k), kzz(i,j,k), \
kxy(i,j,k), kyz(i,j,k), kxz(i,j,k), \
exact_psi, \
exact_psix, exact_psiy, exact_psiz, \
exact_psixx, exact_psiyy, exact_psizz, \
exact_psixy, exact_psiyz, exact_psixz, \
dxgxxjunk, dxgyyjunk, dxgzzjunk, \
dxgxyjunk, dxgyzjunk, dxgxzjunk, \
dygxxjunk, dygyyjunk, dygzzjunk, \
dygxyjunk, dygyzjunk, dygxzjunk, \
dzgxxjunk, dzgyyjunk, dzgzzjunk, \
dzgxyjunk, dzgyzjunk, dzgxzjunk, \
alp(i,j,k), dtalp(i,j,k), \
axjunk, ayjunk, azjunk, \
betax(i,j,k), betay(i,j,k), betaz(i,j,k), \
dtbetax(i,j,k), dtbetay(i,j,k), dtbetaz(i,j,k), \
bxxjunk, bxyjunk, bxzjunk, \
byxjunk, byyjunk, byzjunk, \
bzxjunk, bzyjunk, bzzjunk)
if (cctk_bbox(1) .eq. 1) then
i=1
do j=1,ny
do k=1,nz
EXACTDATAPOINT
end do
end do
end if
if (cctk_bbox(2) .eq. 1) then
i=nx
do j=1,ny
do k=1,nz
EXACTDATAPOINT
end do
end do
end if
if (cctk_bbox(3) .eq. 1) then
j=1
do i=1,nx
do k=1,nz
EXACTDATAPOINT
end do
end do
end if
if (cctk_bbox(4) .eq. 1) then
j=ny
do i=1,nx
do k=1,nz
EXACTDATAPOINT
end do
end do
end if
if (cctk_bbox(5) .eq. 1) then
k=1
do j=1,ny
do i=1,nx
EXACTDATAPOINT
end do
end do
end if
if (cctk_bbox(6) .eq. 1) then
k=nz
do j=1,ny
do i=1,nx
EXACTDATAPOINT
end do
end do
end if
return
end
|
