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
|
! $Header$
#include "cctk.h"
#include "cctk_Parameters.h"
module tensor
implicit none
DECLARE_CCTK_PARAMETERS
private
public calc_trace
public calc_det
public calc_detderiv
public calc_detdot
public calc_inv
public calc_invderiv
public calc_invdot
contains
subroutine calc_trace (kk, gu, trk)
CCTK_REAL, intent(in) :: kk(3,3)
CCTK_REAL, intent(in) :: gu(3,3)
CCTK_REAL, intent(out) :: trk
integer :: i,j
trk = 0
do i=1,3
do j=1,3
trk = trk + gu(i,j) * kk(i,j)
end do
end do
end subroutine calc_trace
subroutine calc_det (gg, dtg)
CCTK_REAL, intent(in) :: gg(3,3)
CCTK_REAL, intent(out) :: dtg
dtg = gg(1,1) * gg(2,2) * gg(3,3) &
& + 2*gg(1,2) * gg(1,3) * gg(2,3) &
& - gg(1,1) * gg(2,3)**2 &
& - gg(2,2) * gg(1,3)**2 &
& - gg(3,3) * gg(1,2)**2
end subroutine calc_det
subroutine calc_pdet (gg, pgg, pdtg)
CCTK_REAL, intent(in) :: gg(3,3), pgg(3,3)
CCTK_REAL, intent(out) :: pdtg
pdtg = pgg(1,1) * gg(2,2) * gg(3,3) &
& + gg(1,1) * pgg(2,2) * gg(3,3) &
& + gg(1,1) * gg(2,2) * pgg(3,3) &
& + 2*pgg(1,2) * gg(1,3) * gg(2,3) &
& + 2* gg(1,2) * pgg(1,3) * gg(2,3) &
& + 2* gg(1,2) * gg(1,3) * pgg(2,3) &
& - pgg(1,1) * gg(2,3)**2 &
& - gg(1,1) * 2*gg(2,3) * pgg(2,3) &
& - pgg(2,2) * gg(1,3)**2 &
& - gg(2,2) * 2*gg(1,3) * pgg(1,3) &
& - pgg(3,3) * gg(1,2)**2 &
& - gg(3,3) * 2*gg(1,2) * pgg(1,2)
end subroutine calc_pdet
subroutine calc_detderiv (gg, dgg, ddtg)
CCTK_REAL, intent(in) :: gg(3,3), dgg(3,3,3)
CCTK_REAL, intent(out) :: ddtg(3)
integer :: i
do i=1,3
call calc_pdet (gg, dgg(:,:,i), ddtg(i))
end do
end subroutine calc_detderiv
subroutine calc_detdot (gg, gg_dot, dtg_dot)
CCTK_REAL, intent(in) :: gg(3,3), gg_dot(3,3)
CCTK_REAL, intent(out) :: dtg_dot
call calc_pdet (gg, gg_dot, dtg_dot)
end subroutine calc_detdot
subroutine calc_inv (gg, dtg, gu)
CCTK_REAL, intent(in) :: gg(3,3), dtg
CCTK_REAL, intent(out) :: gu(3,3)
gu(1,1) = (gg(2,2) * gg(3,3) - gg(2,3) ** 2 ) / dtg
gu(2,2) = (gg(1,1) * gg(3,3) - gg(1,3) ** 2 ) / dtg
gu(3,3) = (gg(1,1) * gg(2,2) - gg(1,2) ** 2 ) / dtg
gu(1,2) = (gg(1,3) * gg(2,3) - gg(1,2) * gg(3,3)) / dtg
gu(1,3) = (gg(1,2) * gg(2,3) - gg(1,3) * gg(2,2)) / dtg
gu(2,3) = (gg(1,3) * gg(1,2) - gg(2,3) * gg(1,1)) / dtg
gu(2,1) = gu(1,2)
gu(3,1) = gu(1,3)
gu(3,2) = gu(2,3)
end subroutine calc_inv
subroutine calc_pinv (gu, pgg, pgu)
CCTK_REAL, intent(in) :: gu(3,3), pgg(3,3)
CCTK_REAL, intent(out) :: pgu(3,3)
integer :: i,j,k,l
do i=1,3
do j=1,3
pgu(i,j) = 0
do k=1,3
do l=1,3
pgu(i,j) = pgu(i,j) - gu(i,k) * gu(j,l) * pgg(k,l)
end do
end do
end do
end do
end subroutine calc_pinv
subroutine calc_invderiv (gu, dgg, dgu)
CCTK_REAL, intent(in) :: gu(3,3), dgg(3,3,3)
CCTK_REAL, intent(out) :: dgu(3,3,3)
integer :: i
do i=1,3
call calc_pinv (gu, dgg(:,:,i), dgu(:,:,i))
end do
end subroutine calc_invderiv
subroutine calc_invdot (gu, gg_dot, gu_dot)
CCTK_REAL, intent(in) :: gu(3,3), gg_dot(3,3)
CCTK_REAL, intent(out) :: gu_dot(3,3)
call calc_pinv (gu, gg_dot, gu_dot)
end subroutine calc_invdot
end module tensor
|
