Add routines for derivatives of the Ricci tensor

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinUtils/TGRtensor/trunk@23 b716e942-a2de-43ad-8f52-f3dfe468e4e7

1 files changed, 69 insertions, 3 deletions

diff --git a/src/ricci.F90 b/src/ricci.F90
index ebe01fc..f14eb49 100644
--- a/src/ricci.F90
+++ b/src/ricci.F90
@@ -7,12 +7,14 @@ module ricci
   private
   public calc_connections
   public calc_connectionderivs
+  public calc_connectionderivs2
   public calc_ricci
+  public calc_ricciderivs
   
 contains
   
-  subroutine calc_connections (gu, dg, gamma)
-    CCTK_REAL, intent(in)  :: gu(3,3), dg(3,3,3)
+  subroutine calc_connections (gu, dgg, gamma)
+    CCTK_REAL, intent(in)  :: gu(3,3), dgg(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)
@@ -22,7 +24,7 @@ contains
              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))
+                     * (dgg(l,j,k) + dgg(l,k,j) - dgg(j,k,l))
              end do
           end do
        end do
@@ -51,6 +53,32 @@ contains
     end do
   end subroutine calc_connectionderivs
   
+  subroutine calc_connectionderivs2 (gu, dgg, dgu, ddgg, ddgu, dddgg, ddgamma)
+    CCTK_REAL, intent(in)  :: gu(3,3), dgg(3,3,3), dgu(3,3,3), ddgg(3,3,3,3), ddgu(3,3,3,3), dddgg(3,3,3,3,3)
+    CCTK_REAL, intent(out) :: ddgamma(3,3,3,3,3)
+    integer   :: i,j,k,l,m,n
+    ! 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
+                do m=1,3
+                   ddgamma(i,j,k,l,m) = 0
+                   do n=1,3
+                      ddgamma(i,j,k,l,m) = ddgamma(i,j,k,l,m) &
+                           + 0.5d0 * ddgu(i,n,l,m) * (dgg(n,j,k) + dgg(n,k,j) - dgg(j,k,n)) &
+                           + 0.5d0 * dgu(i,n,l) * (ddgg(n,j,k,m) + ddgg(n,k,j,m) - ddgg(j,k,n,m)) &
+                           + 0.5d0 * dgu(i,n,m) * (ddgg(n,j,k,l) + ddgg(n,k,j,l) - ddgg(j,k,n,l)) &
+                           + 0.5d0 * gu(i,n) * (dddgg(n,j,k,l,m) + dddgg(n,k,j,l,m) - dddgg(j,k,n,l,m))
+                   end do
+                end do
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_connectionderivs2
+  
   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)
@@ -82,4 +110,42 @@ contains
     if (sum > 1.0e-12 * cnt) call CCTK_WARN (0, "Ricci tensor is not symmetric")
   end subroutine calc_ricci
   
+  subroutine calc_ricciderivs (gamma, dgamma, ddgamma, dri)
+    CCTK_REAL, intent(in)  :: gamma(3,3,3), dgamma(3,3,3,3), ddgamma(3,3,3,3,3)
+    CCTK_REAL, intent(out) :: dri(3,3,3)
+    CCTK_REAL :: sum, cnt
+    integer :: i,j,k,l,m
+    ! 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
+          do k=1,3
+             dri(i,j,k) = 0
+             do l=1,3
+                dri(i,j,k) = dri(i,j,k) + ddgamma(l,i,j,l,k) &
+                     &                  - ddgamma(l,i,l,j,k)
+                do m=1,3
+                   dri(i,j,k) = dri(i,j,k) + dgamma(l,m,l,k) * gamma(m,i,j) &
+                        &                  + gamma(l,m,l) * dgamma(m,i,j,k) &
+                        &                  - dgamma(l,m,j,k) * gamma(m,l,i) &
+                        &                  - gamma(l,m,j) * dgamma(m,l,i,k)
+                end do
+             end do
+          end do
+       end do
+    end do
+    ! check symmetries
+    sum = 0
+    cnt = 0
+    do i=1,3
+       do j=1,3
+          do k=1,3
+             sum = sum + (dri(i,j,k) - dri(j,i,k))**2
+             cnt = cnt + dri(i,j,k)**2
+          end do
+       end do
+    end do
+    if (sum > 1.0e-12 * cnt) call CCTK_WARN (0, "Ricci tensor is not symmetric")
+  end subroutine calc_ricciderivs
+  
 end module ricci