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
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
|
c This subroutine calculates the 4-metric and its inverse at an event,
c taking into account an optional Lorentz boost and an optional rotation.
c The model is first rotated and then boosted, such that the boost is
c applied to the rotated model.
c $Header$
c
c The coordinates are
c Cx(a) = Cactus $x^a$
c Mx(a) = Model $X^a$
c The 4-metrics are
c Cgdd(a,b) = Cactus $g_{ab}$ Cguu(a,b) = Cactus $g^{ab}$
c Mgdd(a,b) = Model $g_{ab}$ Mguu(a,b) = Model $g^{ab}$
c
c For a definition of the Euler angles in the conventions used below, see
c http://mathworld.wolfram.com/EulerAngles.html
c Another useful resource may be
c http://en.wikipedia.org/wiki/Euler_angles
c although this uses (on 2006-11-29) different conventions.
c
c This file is copyright (c) 2003 by Jonathan Thornburg <jthorn@aei.mpg.de>.
c This file is covered by the GNU GPL license; see the files ../README
c and ../COPYING for details.
c
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "param_defs.inc"
subroutine Exact__metric(
$ decoded_exact_model,
$ x, y, z, t,
$ gdtt, gdtx, gdty, gdtz,
$ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
$ gutt, gutx, guty, gutz,
$ guxx, guyy, guzz, guxy, guyz, guxz,
$ psi)
implicit none
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
c input arguments
CCTK_INT decoded_exact_model
CCTK_REAL x, y, z, t
c output arguments
CCTK_REAL gdtt, gdtx, gdty, gdtz,
$ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
$ gutt, gutx, guty, gutz,
$ guxx, guyy, guzz, guxy, guyz, guxz,
$ psi
c intrinsic functions called
CCTK_REAL sqrt
c static local variables describing Lorentz transformation
logical firstcall
data firstcall /.true./
CCTK_REAL gamma
CCTK_REAL vv(3), nn(3)
CCTK_REAL parallel(3,3), perp(3,3)
CCTK_REAL Cx_par(3), Cx_perp(3)
CCTK_REAL partial_Mx_wrt_Cx(0:3,0:3)
CCTK_REAL partial_Cx_wrt_Mx(0:3,0:3)
CCTK_REAL R(0:3,0:3)
save firstcall
save gamma
save vv, nn
save parallel, perp
save Cx_par, Cx_perp
save partial_Mx_wrt_Cx
save partial_Cx_wrt_Mx
save R
c coordinates and 4-metric
CCTK_REAL Cx(0:3)
CCTK_REAL Cgdd(0:3,0:3), Cguu(0:3,0:3)
CCTK_REAL Mx(0:3)
CCTK_REAL Mgdd(0:3,0:3), Mguu(0:3,0:3)
CCTK_REAL Nx(0:3)
CCTK_REAL Ngdd(0:3,0:3), Nguu(0:3,0:3)
c misc temps
CCTK_REAL vnorm, vnormsq
CCTK_REAL delta_ij
CCTK_REAL Cx_par_i, Cx_perp_i
CCTK_REAL vdotCx
CCTK_REAL Cgdd_ab, Cguu_ab
CCTK_REAL cos_phi, sin_phi
CCTK_REAL cos_theta, sin_theta
CCTK_REAL cos_psi, sin_psi
CCTK_REAL R_phi(0:3,0:3), R_theta(0:3,0:3), R_psi(0:3,0:3)
character*1000 warn_buffer
c flags, array indices, etc
logical Tmunu_flag
integer i, j, k, l
integer Ca, Cb, MA, MB
c constants
integer n
parameter (n = 3)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c optimized fast-path if no Lorentz boost and no rotation
c
if ( (boost_vx .eq. 0.0)
$ .and. (boost_vy .eq. 0.0)
$ .and. (boost_vz .eq. 0.0)
$ .and. (rotation_euler_phi .eq. 0.0)
$ .and. (rotation_euler_theta .eq. 0.0)
$ .and. (rotation_euler_psi .eq. 0.0)) then
call Exact__metric_for_model(
$ decoded_exact_model,
$ x, y, z, t,
$ gdtt, gdtx, gdty, gdtz,
$ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
$ gutt, gutx, guty, gutz,
$ guxx, guyy, guzz, guxy, guyz, guxz,
$ psi,
$ Tmunu_flag)
return
end if
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c the rest of this function is the Lorentz-boost case:
c - Lorentz-transform Cactus coordinates --> Model coordinates
c - compute Model 4-metric and inverse at Model coordinates
c - tensor-transform 4-metric and inverse from Model coordinates
c --> Cactus coordinates
c
c All the equations used are given in ../doc/documentation.tex
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c compute Lorentz transformation information on first call
c
if (firstcall) then
firstcall = .false.
c boost velocity
vv(1) = boost_vx
vv(2) = boost_vy
vv(3) = boost_vz
c Lorentz gamma factor, unit vector in direction of boost velocity
vnormsq = 0.0
do i = 1,n
vnormsq = vnormsq + vv(i)*vv(i)
end do
gamma = 1.0 / sqrt(1.0 - vnormsq)
vnorm = sqrt(vnormsq)
if ( (boost_vx .eq. 0.0)
$ .and. (boost_vy .eq. 0.0)
$ .and. (boost_vz .eq. 0.0)) then
nn(1) = 1
nn(2) = 0
nn(3) = 0
else
do i = 1,n
nn(i) = vv(i) / vnorm
end do
end if
c projection operators parallel(*,*) and perp(*,*)
do j = 1,n
do i = 1,n
parallel(i,j) = nn(i) * nn(j)
if (i .eq. j) then
delta_ij = 1.0
else
delta_ij = 0.0
end if
perp(i,j) = delta_ij - parallel(i,j)
end do
end do
c partial derivatives of Model coordinates with respect to Cactus coordinates
partial_Mx_wrt_Cx(0,0) = gamma
do i = 1,n
partial_Mx_wrt_Cx(0,i) = -gamma*vv(i)
end do
do i = 1,n
partial_Mx_wrt_Cx(i,0) = -gamma*vv(i)
do j=1,n
partial_Mx_wrt_Cx(i,j) = gamma*parallel(i,j) + perp(i,j)
end do
end do
c partial derivatives of Cactus coordinates with respect to Model coordinates
partial_Cx_wrt_Mx(0,0) = gamma
do i = 1,n
partial_Cx_wrt_Mx(0,i) = + gamma*vv(i)
end do
do i = 1,n
partial_Cx_wrt_Mx(i,0) = + gamma*vv(i)
do j=1,n
partial_Cx_wrt_Mx(i,j) = gamma*parallel(i,j) + perp(i,j)
end do
end do
c Sines and cosines of rotation angles
cos_phi = cos (rotation_euler_phi)
sin_phi = sin (rotation_euler_phi)
cos_theta = cos (rotation_euler_theta)
sin_theta = sin (rotation_euler_theta)
cos_psi = cos (rotation_euler_psi)
sin_psi = sin (rotation_euler_psi)
c Set up individual rotation matrices
R_phi(0,0) = 1
R_phi(0,1) = 0
R_phi(0,2) = 0
R_phi(0,3) = 0
R_phi(1,0) = 0
R_phi(1,1) = + cos_phi
R_phi(1,2) = + sin_phi
R_phi(1,3) = 0
R_phi(2,0) = 0
R_phi(2,1) = - sin_phi
R_phi(2,2) = + cos_phi
R_phi(2,3) = 0
R_phi(3,0) = 0
R_phi(3,1) = 0
R_phi(3,2) = 0
R_phi(3,3) = 1
R_theta(0,0) = 1
R_theta(0,1) = 0
R_theta(0,2) = 0
R_theta(0,3) = 0
R_theta(1,0) = 0
R_theta(1,1) = 1
R_theta(1,2) = 0
R_theta(1,3) = 0
R_theta(2,0) = 0
R_theta(2,1) = 0
R_theta(2,2) = + cos_theta
R_theta(2,3) = + sin_theta
R_theta(3,0) = 0
R_theta(3,1) = 0
R_theta(3,2) = - sin_theta
R_theta(3,3) = + cos_theta
R_psi(0,0) = 1
R_psi(0,1) = 0
R_psi(0,2) = 0
R_psi(0,3) = 0
R_psi(1,0) = 0
R_psi(1,1) = + cos_psi
R_psi(1,2) = + sin_psi
R_psi(1,3) = 0
R_psi(2,0) = 0
R_psi(2,1) = - sin_psi
R_psi(2,2) = + cos_psi
R_psi(2,3) = 0
R_psi(3,0) = 0
R_psi(3,1) = 0
R_psi(3,2) = 0
R_psi(3,3) = 1
c Combine individual rotation matrices
do i = 0,n
do j = 0,n
R(i,j) = 0
do k = 0,n
do l = 0,n
R(i,j) = R(i,j) + R_psi(i,k) * R_theta(k,l) * R_phi(l,j)
end do
end do
end do
end do
c Notes that help me (Erik Schnetter) think:
c This considers a rotation with phi=0, theta=pi/2, psi=0.
c Nx(i) = Nx(i) + R(j,i) * Mx(j)
c Nx(1) = - Mx(3)
c Nx(3) = Mx(1)
c Mgxx(1,1) = Ngxx(3,3)
c Mgxx(1,3) = - Ngxx(1,3)
c Mgxx(3,3) = Ngxx(1,1)
c Mgxx(i,j) = R(i,k) R(j,l) Ngxx(k,l) [correct]
c Mgxx(i,j) = R(k,i) R(l,j) Ngxx(k,l) [correct]
c Mbetax(1) = - Nbetax(3)
c Mbetax(3) = Nbetax(1)
c Mbetax(i) + R(i,j) Nbetax(j) [wrong]
c Mbetax(i) + R(j,i) Nbetax(j) [correct]
end if
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c compute flat-space components of Cx(*) parallel and perpendicular to vv(*)
c
Cx(0) = t
Cx(1) = x
Cx(2) = y
Cx(3) = z
do i=1,n
Cx_par_i = 0.0
Cx_perp_i = 0.0
do j=1,n
Cx_par_i = Cx_par_i + parallel(i,j)*Cx(j)
Cx_perp_i = Cx_perp_i + perp (i,j)*Cx(j)
end do
Cx_par (i) = Cx_par_i
Cx_perp(i) = Cx_perp_i
end do
c
c Lorentz-transform and rotate the Cactus coordinate
c to get the Model coordinates
c
c Boost
vdotCx = 0.0
do i = 1,n
vdotCx = vdotCx + vv(i)*Cx(i)
end do
Mx(0) = gamma * (Cx(0) - vdotCx)
do i=1,n
Mx(i) = gamma * (Cx_par(i) - vv(i)*Cx(0)) + Cx_perp(i)
end do
c Rotation
do i=0,n
Nx(i) = 0
do j = 0,n
Nx(i) = Nx(i) + R(i,j) * Mx(j)
end do
end do
c
c compute the Model 4-metric and inverse 4-metric at the Model coordinates
c
call Exact__metric_for_model(
$ decoded_exact_model,
$ Nx(1), Nx(2), Nx(3), Nx(0),
$ Ngdd(0,0), Ngdd(0,1), Ngdd(0,2), Ngdd(0,3),
$ Ngdd(1,1), Ngdd(2,2), Ngdd(3,3),
$ Ngdd(1,2), Ngdd(2,3), Ngdd(1,3),
$ Nguu(0,0), Nguu(0,1), Nguu(0,2), Nguu(0,3),
$ Nguu(1,1), Nguu(2,2), Nguu(3,3),
$ Nguu(1,2), Nguu(2,3), Nguu(1,3),
$ psi, Tmunu_flag)
if (Tmunu_flag) then
write (warn_buffer, '(a,i8,a,a)')
$ 'exact_model = ', decoded_exact_model,
$ 'sets the stress-energy tensor',
$ ' ==> we cannot Lorentz-boost or rotate it!'
call CCTK_WARN(0, warn_buffer)
end if
c
c symmetrize the Model 4-metric and inverse 4-metric arrays
c (the Exact__metric_for_model() call only set the upper triangles)
c
Ngdd(1,0) = Ngdd(0,1)
Ngdd(2,0) = Ngdd(0,2)
Ngdd(2,1) = Ngdd(1,2)
Ngdd(3,0) = Ngdd(0,3)
Ngdd(3,1) = Ngdd(1,3)
Ngdd(3,2) = Ngdd(2,3)
Nguu(1,0) = Nguu(0,1)
Nguu(2,0) = Nguu(0,2)
Nguu(2,1) = Nguu(1,2)
Nguu(3,0) = Nguu(0,3)
Nguu(3,1) = Nguu(1,3)
Nguu(3,2) = Nguu(2,3)
c
c tensor-transorm (the upper triangle of) the 4-metric and inverse 4-metric
c
c Rotations
do i = 0,n
do j = 0,n
Mgdd(i,j) = 0
Mguu(i,j) = 0
do k = 0,n
do l = 0,n
Mgdd(i,j) = Mgdd(i,j) + R(k,i) * R(l,j) * Ngdd(k,l)
c The inverse of R is also its transpose. That means that the
c transpose of the inverse, which you would use for g^ij, is just R
c again.
Mguu(i,j) = Mguu(i,j) + R(k,i) * R(l,j) * Nguu(k,l)
end do
end do
end do
end do
c Boost
do Ca = 0,n
do Cb = Ca,n
Cgdd_ab = 0.0
do Ma = 0,n
do Mb = 0,n
Cgdd_ab = Cgdd_ab
$ + Mgdd(Ma,Mb)
$ * partial_Mx_wrt_Cx(Ma,Ca)
$ * partial_Mx_wrt_Cx(Mb,Cb)
end do
end do
Cgdd(Ca,Cb) = Cgdd_ab
end do
end do
do Ca = 0,n
do Cb = Ca,n
Cguu_ab = 0.0
do Ma = 0,n
do Mb = 0,n
Cguu_ab = Cguu_ab
$ + Mguu(Ma,Mb)
$ * partial_Cx_wrt_Mx(Ca,Ma)
$ * partial_Cx_wrt_Mx(Cb,Mb)
end do
end do
Cguu(Ca,Cb) = Cguu_ab
end do
end do
c
c unpack the Cactus-coordinates 4-metric and inverse 4-metric
c into the corresponding output arguments
c
gdtt = Cgdd(0,0)
gdtx = Cgdd(0,1)
gdty = Cgdd(0,2)
gdtz = Cgdd(0,3)
gdxx = Cgdd(1,1)
gdxy = Cgdd(1,2)
gdxz = Cgdd(1,3)
gdyy = Cgdd(2,2)
gdyz = Cgdd(2,3)
gdzz = Cgdd(3,3)
gutt = Cguu(0,0)
gutx = Cguu(0,1)
guty = Cguu(0,2)
gutz = Cguu(0,3)
guxx = Cguu(1,1)
guxy = Cguu(1,2)
guxz = Cguu(1,3)
guyy = Cguu(2,2)
guyz = Cguu(2,3)
guzz = Cguu(3,3)
end
|