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
|
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"
SUBROUTINE qlm_calc_3determinant (CCTK_ARGUMENTS, hn)
USE adm_metric
USE cctk
USE qlm_boundary
USE qlm_derivs
USE qlm_variables
USE tensor
IMPLICIT NONE
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
INTEGER :: hn
CCTK_REAL :: gg(3,3)
CCTK_REAL :: alfa
CCTK_REAL :: beta(3)
CCTK_REAL :: g4(0:3,0:3)
CCTK_COMPLEX :: mm(0:3)
CCTK_REAL :: rr(3), rr_p(3), rr_p_p(3)
CCTK_REAL :: Ttilde(0:3)
CCTK_REAL :: a
CCTK_COMPLEX :: d
CCTK_REAL :: t0, t1, t2
logical :: ce0, ce1, ce2
CCTK_REAL :: delta_space(2)
INTEGER :: i, j
CCTK_REAL :: theta, phi
IF (veryverbose/=0) THEN
CALL CCTK_INFO ("Computing 3-Volume Element of Surface")
END IF
t0 = qlm_time(hn)
t1 = qlm_time_p(hn)
t2 = qlm_time_p_p(hn)
ce0 = qlm_have_valid_data(hn) == 0
ce1 = qlm_have_valid_data_p(hn) == 0
ce2 = qlm_have_valid_data_p_p(hn) == 0
delta_space(:) = (/ qlm_delta_theta(hn), qlm_delta_phi(hn) /)
! Calculate the coordinates
DO j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn)
DO i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn)
theta = qlm_origin_theta(hn) + (i-1)*qlm_delta_theta(hn)
phi = qlm_origin_phi(hn) + (j-1)*qlm_delta_phi(hn)
! Get the stuff from the arrays
gg(1,1) = qlm_gxx(i,j)
gg(1,2) = qlm_gxy(i,j)
gg(1,3) = qlm_gxz(i,j)
gg(2,2) = qlm_gyy(i,j)
gg(2,3) = qlm_gyz(i,j)
gg(3,3) = qlm_gzz(i,j)
gg(2,1) = gg(1,2)
gg(3,1) = gg(1,3)
gg(3,2) = gg(2,3)
alfa = qlm_alpha(i,j)
beta(1) = qlm_betax(i,j)
beta(2) = qlm_betay(i,j)
beta(3) = qlm_betaz(i,j)
mm(0) = qlm_m0(i,j,hn)
mm(1) = qlm_m1(i,j,hn)
mm(2) = qlm_m2(i,j,hn)
mm(3) = qlm_m3(i,j,hn)
! Build the four-metric
CALL calc_4metric (gg,alfa,beta, g4)
! Find a 3rd vector of triad
rr(1) = qlm_x(i,j,hn)
rr(2) = qlm_y(i,j,hn)
rr(3) = qlm_z(i,j,hn)
rr_p(1) = qlm_x_p(i,j,hn)
rr_p(2) = qlm_y_p(i,j,hn)
rr_p(3) = qlm_z_p(i,j,hn)
rr_p_p(1) = qlm_x_p_p(i,j,hn)
rr_p_p(2) = qlm_y_p_p(i,j,hn)
rr_p_p(3) = qlm_z_p_p(i,j,hn)
Ttilde(0) = 1
Ttilde(1:3) = timederiv (rr, rr_p, rr_p_p, t0,t1,t2, ce0,ce1,ce2)
! Compute some scalar products
a=SUM(MATMUL(g4, Ttilde)*Ttilde)
d=SUM(MATMUL(g4, mm)*Ttilde)
! This is the determinant
qlm_3det(i,j,hn)=a-2*CONJG(d)*d
END DO
END DO
CALL set_boundary (CCTK_PASS_FTOF, hn, qlm_3det(:,:,hn), +1)
END SUBROUTINE qlm_calc_3determinant
|
