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
|
! $Header$
#include "cctk.h"
module ricci
implicit none
private
public calc_connections
public calc_connectionderivs
public calc_ricci
contains
subroutine calc_connections (gu, dg, gamma)
CCTK_REAL, intent(in) :: gu(3,3), dg(3,3,3)
CCTK_REAL, intent(out) :: gamma(3,3,3)
integer :: i,j,k,l
! Gamma^i_jk = 1/2 g^il (g_lj,k + g_lk,j - g_jk,l)
do i=1,3
do j=1,3
do k=1,3
gamma(i,j,k) = 0
do l=1,3
gamma(i,j,k) = gamma(i,j,k) + 0.5d0 * gu(i,l) &
* (dg(l,j,k) + dg(l,k,j) - dg(j,k,l))
end do
end do
end do
end do
end subroutine calc_connections
subroutine calc_connectionderivs (gu, dgg, dgu, ddgg, dgamma)
CCTK_REAL, intent(in) :: gu(3,3), dgg(3,3,3), dgu(3,3,3), ddgg(3,3,3,3)
CCTK_REAL, intent(out) :: dgamma(3,3,3,3)
integer :: i,j,k,l,m
! Gamma^i_jk,l = 1/2 g^im,l (g_mj,k + g_mk,j - g_jk,m)
! + 1/2 g^im (g_mj,kl + g_mk,jl - g_jk,ml)
do i=1,3
do j=1,3
do k=1,3
do l=1,3
dgamma(i,j,k,l) = 0
do m=1,3
dgamma(i,j,k,l) = dgamma(i,j,k,l) &
+ 0.5d0 * dgu(i,m,l) * (dgg(m,j,k) + dgg(m,k,j) - dgg(j,k,m)) &
+ 0.5d0 * gu(i,m) * (ddgg(m,j,k,l) + ddgg(m,k,j,l) - ddgg(j,k,m,l))
end do
end do
end do
end do
end do
end subroutine calc_connectionderivs
subroutine calc_ricci (gamma, dgamma, ri)
CCTK_REAL, intent(in) :: gamma(3,3,3), dgamma(3,3,3,3)
CCTK_REAL, intent(out) :: ri(3,3)
CCTK_REAL :: sum, cnt
integer :: i,j,k,l
! R_ij = Gamma^k_ij,k - Gamma^k_ik,j
! + Gamma^k_lk Gamma^l_ij - Gamma^k_lj Gamma^l_ki
do i=1,3
do j=1,3
ri(i,j) = 0
do k=1,3
ri(i,j) = ri(i,j) + dgamma(k,i,j,k) - dgamma(k,i,k,j)
do l=1,3
ri(i,j) = ri(i,j) + gamma(k,l,k) * gamma(l,i,j) &
& - gamma(k,l,j) * gamma(l,k,i)
end do
end do
end do
end do
! check symmetries
sum = 0
cnt = 0
do i=1,3
do j=1,3
sum = sum + (ri(i,j) - ri(j,i))**2
cnt = cnt + ri(i,j)**2
end do
end do
if (sum > 1.0e-12 * cnt) call CCTK_WARN (0, "Ricci tensor is not symmetric")
end subroutine calc_ricci
end module ricci
|
