New function calc_3metricderivdot

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

1 files changed, 189 insertions, 19 deletions

diff --git a/src/adm_metric.F90 b/src/adm_metric.F90
index 9e284cd..a9ef4a4 100644
--- a/src/adm_metric.F90
+++ b/src/adm_metric.F90
@@ -3,6 +3,7 @@
 #include "cctk.h"
 
 module adm_metric
+  use constants
   use tensor
   implicit none
   private
@@ -12,14 +13,19 @@ module adm_metric
   public calc_3metric
   public calc_3metricderivs
   
-  public calc_3metricdot
   public calc_extcurv
+  public calc_3metricdot
+  public calc_3metricderivdot
+!!$  public calc_extcurvdot
+!!$  public calc_3metricdot2
   
   public calc_3metricdot_simple
   public calc_3metricderivdot_simple
   public calc_extcurvdot_simple
   public calc_3metricdot2_simple
   
+!!$  public simplify_t4
+  
 contains
   
   subroutine calc_4metric (gg, alfa, beta, g4)
@@ -235,6 +241,29 @@ contains
   
   
   
+  subroutine calc_extcurv (gg, dgg, gg_dot, alfa, beta, dbeta, kk)
+    CCTK_REAL, intent(in)  :: gg(3,3), dgg(3,3,3), gg_dot(3,3)
+    CCTK_REAL, intent(in)  :: alfa
+    CCTK_REAL, intent(in)  :: beta(3), dbeta(3,3)
+    CCTK_REAL, intent(out) :: kk(3,3)
+    integer :: i,j,k
+    
+    ! d/dt g_ij = -2 alpha K_ij + g_kj beta^k,i + g_ik beta^k,j + beta^k g_ij,k
+    
+    do i=1,3
+       do j=1,3
+          kk(i,j) = - gg_dot(i,j)
+          do k=1,3
+             kk(i,j) = kk(i,j) &
+                  + gg(k,j) * dbeta(k,i) + gg(i,k) * dbeta(k,j) &
+                  + beta(k) * dgg(i,j,k)
+          end do
+          kk(i,j) = kk(i,j) / (2*alfa)
+       end do
+    end do
+    
+  end subroutine calc_extcurv
+  
   subroutine calc_3metricdot (gg, dgg, kk, alfa, beta, dbeta, gg_dot)
     CCTK_REAL, intent(in)  :: gg(3,3), dgg(3,3,3)
     CCTK_REAL, intent(in)  :: kk(3,3)
@@ -257,31 +286,95 @@ contains
     end do
   end subroutine calc_3metricdot
   
-  subroutine calc_extcurv (gg, dgg, gg_dot, alfa, beta, dbeta, kk)
-    CCTK_REAL, intent(in)  :: gg(3,3), dgg(3,3,3), gg_dot(3,3)
-    CCTK_REAL, intent(in)  :: alfa
-    CCTK_REAL, intent(in)  :: beta(3), dbeta(3,3)
-    CCTK_REAL, intent(out) :: kk(3,3)
-    integer :: i,j,k
+  subroutine calc_3metricderivdot (gg, dgg, ddgg, kk, dkk, &
+       alfa, dalfa, beta, dbeta, ddbeta, dgg_dot)
+    CCTK_REAL, intent(in)  :: gg(3,3), dgg(3,3,3), ddgg(3,3,3,3)
+    CCTK_REAL, intent(in)  :: kk(3,3), dkk(3,3,3)
+    CCTK_REAL, intent(in)  :: alfa, dalfa(3)
+    CCTK_REAL, intent(in)  :: beta(3), dbeta(3,3), ddbeta(3,3,3)
+    CCTK_REAL, intent(out) :: dgg_dot(3,3,3)
+    integer :: i,j,k,l
     
     ! d/dt g_ij = -2 alpha K_ij + g_kj beta^k,i + g_ik beta^k,j + beta^k g_ij,k
     
+    ! d/dt g_ij,k = - 2 alpha,k K_ij - 2 alpha K_ij,k
+    !               + g_lj,k beta^l,i + g_lj beta^l,ik
+    !               + g_il,k beta^l,j + g_il beta^l,jk
+    !               + beta^l,k g_ij,l + beta^l g_ij,lk
+    
     do i=1,3
        do j=1,3
-          kk(i,j) = - gg_dot(i,j)
           do k=1,3
-             kk(i,j) = kk(i,j) &
-                  + gg(k,j) * dbeta(k,i) + gg(i,k) * dbeta(k,j) &
-                  + beta(k) * dgg(i,j,k)
+             dgg_dot(i,j,k) = - 2*dalfa(k) * kk(i,j) - 2*alfa * dkk(i,j,k)
+             do l=1,3
+                dgg_dot(i,j,k) = dgg_dot(i,j,k) &
+                     + dgg(l,j,k) * dbeta(l,i) + gg(l,j) * ddbeta(l,i,k) &
+                     + dgg(i,l,k) * dbeta(l,j) + gg(i,l) * ddbeta(l,j,k) &
+                     + dbeta(l,k) * dgg(i,j,l) + beta(l) * ddgg(i,j,l,k)
+             end do
           end do
-          kk(i,j) = kk(i,j) / (2*alfa)
        end do
     end do
-    
-  end subroutine calc_extcurv
+  end subroutine calc_3metricderivdot
+  
+!!$  subroutine calc_extcurvdot (gg,gu,ri, kk, t44, kk_dot)
+!!$    CCTK_REAL, intent(in)  :: gg(3,3), gu(3,3)
+!!$    CCTK_REAL, intent(in)  :: ri(3,3)
+!!$    CCTK_REAL, intent(in)  :: kk(3,3)
+!!$    CCTK_REAL, intent(out) :: kk_dot(3,3)
+!!$    integer :: i,j,l,m
+!!$    
+!!$    ! d/dt K_ij = R_ij - 2 K_il g^lm K_mj + g^lm K_lm K_ij
+!!$    !                  - 8 pi (T_ij - 1/2 g_ij T)
+!!$    
+!!$    do i=1,3
+!!$       do j=1,3
+!!$          kk_dot(i,j) = ri(i,j) - 8*pi * tt(i,j)
+!!$          do l=1,3
+!!$             do m=1,3
+!!$                kk_dot(i,j) = kk_dot(i,j) &
+!!$                     + gu(l,m) * (-2 * kk(i,l) * kk(m,j) + kk(l,m) * kk(i,j)) &
+!!$                     + 4*pi * gg(i,j) * gu(l,m) * tt(l,m)
+!!$             end do
+!!$          end do
+!!$       end do
+!!$    end do
+!!$  end subroutine calc_extcurvdot_simple
+!!$  
+!!$  subroutine calc_3metricdot2 (kk_dot, gg_dot2)
+!!$    CCTK_REAL, intent(in)  :: kk_dot(3,3)
+!!$    CCTK_REAL, intent(out) :: gg_dot2(3,3)
+!!$    integer :: i,j
+!!$    
+!!$    ! d/dt g_ij = -2 alpha K_ij + g_kj beta^k,i + g_ik beta^k,j + beta^k g_ij,k
+!!$    
+!!$    ! d^2/dt^2 g_ij = -2 alpha K_ij + g_kj beta^k,i + g_ik beta^k,j + beta^k g_ij,k
+!!$    
+!!$    do i=1,3
+!!$       do j=1,3
+!!$          gg_dot2(i,j) = -2 * kk_dot(i,j)
+!!$       end do
+!!$    end do
+!!$  end subroutine calc_3metricdot2
   
   
   
+  subroutine calc_extcurv_simple (gg_dot, kk)
+    CCTK_REAL, intent(in)  :: gg_dot(3,3)
+    CCTK_REAL, intent(out) :: kk(3,3)
+    integer :: i,j
+    
+    ! d/dt g_ij = -2 K_ij
+    
+    do i=1,3
+       do j=1,3
+          kk(i,j) = - gg_dot(i,j)
+          kk(i,j) = kk(i,j) / 2
+       end do
+    end do
+    
+  end subroutine calc_extcurv_simple
+  
   subroutine calc_3metricdot_simple (kk, gg_dot)
     CCTK_REAL, intent(in)  :: kk(3,3)
     CCTK_REAL, intent(out) :: gg_dot(3,3)
@@ -312,22 +405,25 @@ contains
     end do
   end subroutine calc_3metricderivdot_simple
   
-  subroutine calc_extcurvdot_simple (gu,ri, kk, kk_dot)
-    CCTK_REAL, intent(in)  :: gu(3,3)
+  subroutine calc_extcurvdot_simple (gg,gu,ri, kk, tt, kk_dot)
+    CCTK_REAL, intent(in)  :: gg(3,3), gu(3,3)
     CCTK_REAL, intent(in)  :: ri(3,3)
     CCTK_REAL, intent(in)  :: kk(3,3)
+    CCTK_REAL, intent(in)  :: tt(3,3)
     CCTK_REAL, intent(out) :: kk_dot(3,3)
     integer :: i,j,l,m
     
     ! d/dt K_ij = R_ij - 2 K_il g^lm K_mj + g^lm K_lm K_ij
+    !                  - 8 pi (T_ij - 1/2 g_ij T)
     
     do i=1,3
        do j=1,3
-          kk_dot(i,j) = ri(i,j)
+          kk_dot(i,j) = ri(i,j) - 8*pi * tt(i,j)
           do l=1,3
              do m=1,3
-                kk_dot(i,j) = kk_dot(i,j) + gu(l,m) &
-                     * (-2 * kk(i,l) * kk(m,j) + kk(l,m) * kk(i,j))
+                kk_dot(i,j) = kk_dot(i,j) &
+                     + gu(l,m) * (-2 * kk(i,l) * kk(m,j) + kk(l,m) * kk(i,j)) &
+                     + 4*pi * gg(i,j) * gu(l,m) * tt(l,m)
              end do
           end do
        end do
@@ -349,4 +445,78 @@ contains
     end do
   end subroutine calc_3metricdot2_simple
   
+  
+  
+!!$  subroutine simplify_t4 (t4, alfa, beta, t4s)
+!!$    CCTK_REAL, intent(in)  :: t4(0:3,0:3)
+!!$    CCTK_REAL, intent(in)  :: alfa
+!!$    CCTK_REAL, intent(in)  :: beta(3)
+!!$    CCTK_REAL, intent(out) :: t4s(0:3,0:3)
+!!$    
+!!$    ! T_ab   g_ij alfa beta^i
+!!$    ! T'_ab   g_ij 1 [0,0,0]
+!!$    
+!!$    ! x^a
+!!$    ! x'^a
+!!$    
+!!$    ! g'_ab = dx^c/dx'^a dx^d/dx'^b g_cd
+!!$    !       = dx^0/dx'^a dx^0/dx'^b g_00
+!!$    !         dx^0/dx'^a dx^l/dx'^b g_0l
+!!$    !         dx^k/dx'^a dx^0/dx'^b g_k0
+!!$    !         dx^k/dx'^a dx^l/dx'^b g_kl
+!!$    
+!!$    ! g'_ab = J^c_a J^d_b g_cd
+!!$    ! T'_ab = J^c_a J^d_b T_cd
+!!$    
+!!$    
+!!$    
+!!$    ! delta^i_b = g'^ia g'_ab = J^c_a J^d_b g_cd g'^ia
+!!$    
+!!$    
+!!$    
+!!$    ! g' = J^T g J
+!!$    ! g' J^(-1) = J^T g
+!!$    ! T' = J^T T J
+!!$    
+!!$    
+!!$    
+!!$    ! -1 = g'_00 = dx^0/dx'^0 dx^0/dx'^0 g_00
+!!$    !              dx^0/dx'^0 dx^l/dx'^0 g_0l
+!!$    !              dx^k/dx'^0 dx^0/dx'^0 g_k0
+!!$    !              dx^k/dx'^0 dx^l/dx'^0 g_kl
+!!$    
+!!$    ! 0 = g'_0j = dx^0/dx'^0 dx^0/dx'^j g_00
+!!$    !             dx^0/dx'^0 dx^l/dx'^j g_0l   !!! <<< !!!
+!!$    !             dx^k/dx'^0 dx^0/dx'^j g_k0
+!!$    !             dx^k/dx'^0 dx^l/dx'^j g_kl
+!!$    
+!!$    ! g_ij = g'_ij = dx^0/dx'^i dx^0/dx'^j g_00
+!!$    !                dx^0/dx'^i dx^l/dx'^j g_0l
+!!$    !                dx^k/dx'^i dx^0/dx'^j g_k0
+!!$    !                dx^k/dx'^i dx^l/dx'^j g_kl
+!!$    
+!!$    ! dx^0/dx'^0 = 1/sqrt(-g_00)
+!!$    ! dx^0/dx'^i = 0
+!!$    ! dx^k/dx'^0 = 0
+!!$    ! dx^k/dx'^i = delta^k_i
+!!$    
+!!$    
+!!$    
+!!$    ! g_00 = -alfa^2 + gama_ij beta^i beta^j
+!!$    ! g_0i = gama_ij beta^j
+!!$    
+!!$    ! alfa^2 = -g_00 + gama_ij beta^i beta^j
+!!$    ! beta^i = gama^ij g_0j
+!!$    
+!!$    ! alfa'^2 = -g'_00 + gama'_ij beta'^i beta'^j
+!!$    ! beta'^i = gama'^ij g'_0j
+!!$    
+!!$    ! -1 = alfa'^2 = -g'_00 + gama'_ij beta'^i beta'^j
+!!$    !  0 = beta'^i = gama'^ij g'_0j
+!!$    
+!!$    ! g'_00 = -1
+!!$    ! g'_0i =  0
+!!$    
+!!$  end subroutine simplify_t4
+  
 end module adm_metric