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
|
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"
subroutine qlm_killing_transport (CCTK_ARGUMENTS, hn)
use cctk
use constants
use lapack
use qlm_boundary
use qlm_killing_transportation
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
integer :: hn
interface
integer function TAT_isnan (x)
implicit none
CCTK_REAL x
end function TAT_isnan
end interface
integer, parameter :: lik = lapack_integer_kind
CCTK_REAL :: xi(2,3), chi(3)
CCTK_REAL :: vec(3,3)
CCTK_REAL :: wr(3), wi(3), vl(3,3), vr(3,3)
integer :: i0, j0
integer :: n
integer(lik) :: info
character :: msg*1000
integer(lik), parameter :: lwork = 100
CCTK_REAL :: work(lwork)
if (veryverbose/=0) then
call CCTK_INFO ("Transporting Killing vector field")
end if
! latitude of "equator"
i0 = (qlm_ntheta(hn)+1)/2
! longitude of zero meridian
j0 = 1+qlm_nghostsphi(hn)
vec = delta3
if (veryverbose/=0) then
write (msg, '("Initial vectors (xi^t xi^p chi):")')
call CCTK_INFO (msg)
do n=1,3
write (msg, '(i2,3g18.8)') n, vec(:,n)
call CCTK_INFO (msg)
end do
end if
xi(1,:) = vec(1,:)
xi(2,:) = vec(2,:)
chi(:) = vec(3,:)
do n=1,3
call transport_along_equator (CCTK_PASS_FTOF, hn, i0, xi(:,n), chi(n))
end do
vec(1,:) = xi(1,:)
vec(2,:) = xi(2,:)
vec(3,:) = chi(:)
if (veryverbose/=0) then
write (msg, '("Final vectors (xi^t xi^p chi):")')
call CCTK_INFO (msg)
do n=1,3
write (msg, '(i2,3g18.8)') n, vec(:,n)
call CCTK_INFO (msg)
end do
end if
if (TAT_isnan(sum(vec)) /= 0) then
! qlm_calc_error(hn) = 1
qlm_have_killing_vector(hn) = 0
call CCTK_WARN (3, "There are nans in the final vectors")
goto 9999
end if
call geev ('n', 'v', 3_lik, vec, 3_lik, wr, wi, vl, 3_lik, vr, 3_lik, &
work, lwork, info)
if (info/=0) then
! qlm_calc_error(hn) = 1
qlm_have_killing_vector(hn) = 0
write (msg, '("Error in call to GEEV, info=",i2)') info
call CCTK_WARN (3, msg)
goto 9999
end if
if (veryverbose/=0) then
write (msg, '("Eigenvectors and eigenvalues:")')
call CCTK_INFO (msg)
do n=1,3
write (msg, '(i2,3g18.8,(" (",g18.8,",",g18.8,")"))') &
n, vr(:,n), wr(n), wi(n)
call CCTK_INFO (msg)
end do
end if
xi(1,:) = vr(1,:)
xi(2,:) = vr(2,:)
chi(:) = vr(3,:)
! TODO:
! This transport scheme is not ideal.
! It leads to large fluctuations in the phi direction.
n=3
qlm_killing_eigenvalue_re(hn) = wr(n)
qlm_killing_eigenvalue_im(hn) = wi(n)
if (abs(cmplx(wr(n),wi(n),kind(wr)) - (1,0)) > 1.0d-4) then
call CCTK_WARN (3, "Did not manage to find an eigenvector with the eigenvalue 1")
end if
call transport_along_equator (CCTK_PASS_FTOF, hn, i0, xi(:,n), chi(n))
call transport_along_meridians (CCTK_PASS_FTOF, hn, i0)
call set_boundary (CCTK_PASS_FTOF, hn, qlm_xi_t(:,:,hn), -1)
call set_boundary (CCTK_PASS_FTOF, hn, qlm_xi_p(:,:,hn), -1)
call set_boundary (CCTK_PASS_FTOF, hn, qlm_chi (:,:,hn), +1)
9999 continue
end subroutine qlm_killing_transport
|
