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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
|
C This routine sets the lapse and/or shift by calling a routine
C that does it pointwise. Note that it could be easily modified
C to set the Bona-Masso variables B_xx etc.
C $Header$
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
subroutine Exact__gauge(CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_FUNCTIONS
integer i,j,k
integer nx,ny,nz
logical set_lapse, set_dtlapse, set_shift, set_dtshift
CCTK_REAL tt, xx, yy, zz
CCTK_REAL gxxtmp, gyytmp, gzztmp,
$ gxytmp, gyztmp, gxztmp,
$ hxxtmp, hyytmp, hzztmp,
$ hxytmp, hyztmp, hxztmp,
$ dxgxxtmp, dxgyytmp, dxgzztmp,
$ dxgxytmp, dxgyztmp, dxgxztmp,
$ dygxxtmp, dygyytmp, dygzztmp,
$ dygxytmp, dygyztmp, dygxztmp,
$ dzgxxtmp, dzgyytmp, dzgzztmp,
$ dzgxytmp, dzgyztmp, dzgxztmp,
$ alptmp, dtalptmp, axtmp, aytmp, aztmp,
$ betaxtmp, betaytmp, betaztmp,
$ dtbetaxtmp, dtbetaytmp, dtbetaztmp,
$ bxxtmp, bxytmp, bxztmp,
$ byxtmp, byytmp, byztmp,
$ bzxtmp, bzytmp, bzztmp
CCTK_REAL
$ exact_psi,
$ exact_psix, exact_psiy, exact_psiz,
$ exact_psixx, exact_psiyy, exact_psizz,
$ exact_psixy, exact_psiyz, exact_psixz
LOGICAL is_initial_slice, is_later_slice
C are we on the initial slice or some later slice?
C n.b. the logical expressions later in this function involving
C these flags below would be *so* much nicer if Fortran
C grokked C-style conditional expressions... :) :)
is_initial_slice = cctk_iteration .eq. 0
is_later_slice = cctk_iteration .ne. 0
C Grid parameters.
nx = cctk_lsh(1)
ny = cctk_lsh(2)
nz = cctk_lsh(3)
C This code used to set t = time + dt/2 to get 2nd order accuracy,
C but this leads to the initial data being set at the wrong time. :(
C In the context of MoL, we want to set variables at the standard Cactus
C time (cctk_time), because MoL takes care of calling us at each MoL
C iteration, and updating the field variables appropriately.
C
C Alas, setting at cctk_time probably gives O(dt) errors for non-MoL
C evoutions where Exact is used to set stuff at each time step.
C Fixing this [unless we just declare all non-MoL stuff obselete :) ]
C probably requires cleaning up our (++messy) schedule.ccl , which
C is why this remains a bug for now... :( :(
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Set lapse and/or shift?
C
if ( is_initial_slice ) then
set_lapse = CCTK_Equals(initial_lapse, "exact").ne.0
set_shift = CCTK_Equals(initial_shift, "exact").ne.0
set_dtlapse = CCTK_Equals(initial_dtlapse, "exact").ne.0
set_dtshift = CCTK_Equals(initial_dtshift, "exact").ne.0
end if
if ( is_later_slice ) then
set_lapse = CCTK_Equals(lapse_evolution_method, "exact").ne.0
set_shift = CCTK_Equals(shift_evolution_method, "exact").ne.0
set_dtlapse = CCTK_Equals(dtlapse_evolution_method, "exact").ne.0
set_dtshift = CCTK_Equals(dtshift_evolution_method, "exact").ne.0
end if
if ( set_lapse .or. set_shift .or. set_dtlapse .or. set_dtshift) then
C$omp parallel do private (
C$omp$ i, j, k,
C$omp$ tt, xx, yy, zz,
C$omp$ alptmp, dtalptmp, axtmp, aytmp, aztmp,
C$omp$ betaxtmp, betaytmp, betaztmp,
C$omp$ dtbetaxtmp, dtbetaytmp, dtbetaztmp,
C$omp$ bxxtmp, bxytmp, bxztmp,
C$omp$ byxtmp, byytmp, byztmp,
C$omp$ bzxtmp, bzytmp, bzztmp,
C$omp$ dxgxxtmp, dxgyytmp, dxgzztmp,
C$omp$ dxgxytmp, dxgyztmp, dxgxztmp,
C$omp$ dygxxtmp, dygyytmp, dygzztmp,
C$omp$ dygxytmp, dygyztmp, dygxztmp,
C$omp$ dzgxxtmp, dzgyytmp, dzgzztmp,
C$omp$ dzgxytmp, dzgyztmp, dzgxztmp,
C$omp$ exact_psi,
C$omp$ exact_psix, exact_psiy, exact_psiz,
C$omp$ exact_psixx, exact_psiyy, exact_psizz,
C$omp$ exact_psixy, exact_psiyz, exact_psixz)
do k=1,nz
do j=1,ny
do i=1,nx
tt = cctk_time
xx = x(i,j,k) - cctk_time * shift_add_x
yy = y(i,j,k) - cctk_time * shift_add_y
zz = z(i,j,k) - cctk_time * shift_add_z
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
call Exact__Bona_Masso_data(
$ decoded_exact_model,
$ xx, yy, zz, tt,
$ gxxtmp, gyytmp, gzztmp,
$ gxytmp, gyztmp, gxztmp,
$ hxxtmp, hyytmp, hzztmp,
$ hxytmp, hyztmp, hxztmp,
$ exact_psi,
$ exact_psix, exact_psiy, exact_psiz,
$ exact_psixx, exact_psiyy, exact_psizz,
$ exact_psixy, exact_psiyz, exact_psixz,
$ dxgxxtmp, dxgyytmp, dxgzztmp,
$ dxgxytmp, dxgyztmp, dxgxztmp,
$ dygxxtmp, dygyytmp, dygzztmp,
$ dygxytmp, dygyztmp, dygxztmp,
$ dzgxxtmp, dzgyytmp, dzgzztmp,
$ dzgxytmp, dzgyztmp, dzgxztmp,
$ alptmp, dtalptmp, axtmp, aytmp, aztmp,
$ betaxtmp, betaytmp, betaztmp,
$ dtbetaxtmp, dtbetaytmp, dtbetaztmp,
$ bxxtmp, bxytmp, bxztmp,
$ byxtmp, byytmp, byztmp,
$ bzxtmp, bzytmp, bzztmp)
if ( set_lapse ) then
alp(i,j,k) = alptmp
end if
if ( set_shift ) then
betax(i,j,k) = betaxtmp + shift_add_x
betay(i,j,k) = betaytmp + shift_add_y
betaz(i,j,k) = betaztmp + shift_add_z
end if
if ( set_dtlapse ) then
dtalp(i,j,k) = dtalptmp
end if
if ( set_dtshift ) then
dtbetax(i,j,k) = dtbetaxtmp
dtbetay(i,j,k) = dtbetaytmp
dtbetaz(i,j,k) = dtbetaztmp
end if
end do
end do
end do
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
else
call CCTK_WARN(1,'Exact__gauge has been called without doing anything')
end if
return
end
|
