Adding shift terms.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/ADMConstraints/trunk@42 b7a48df3-cbbf-4440-997f-b4b717c9f7fc

1 files changed, 73 insertions, 11 deletions

diff --git a/src/ADMConstraints.F b/src/ADMConstraints.F
index 1018d08..1a999dd 100644
--- a/src/ADMConstraints.F
+++ b/src/ADMConstraints.F
@@ -56,8 +56,9 @@ c     Return code from Cactus boundary conditions
       integer ierr
 
       CCTK_REAL :: dx,dy,dz
-
-      CCTK_REAL :: pi,det,uxx,uxy,uxz,uyy,uyz,uzz
+      CCTK_REAL :: det,uxx,uxy,uxz,uyy,uyz,uzz
+      CCTK_REAL :: m_rho,m_sx,m_sy,m_sz
+      CCTK_REAL :: pi,ialp,ialp2
 
 c     Temporaries for the Stress-Energy tensor
       CCTK_REAL :: Ttt,Ttx,Tty,Ttz,Txx,Txy,Txz,Tyy,Tyz,Tzz
@@ -94,10 +95,13 @@ c     Macros from Standard Einstein
 
       Pi =  ACOS(-1D0)
 
-      do k = k1, k2
-         do j = j1, j2
-            do i = i1, i2 
-               
+      do k = k1,k2
+         do j = j1,j2
+            do i = i1,i2
+
+               ialp  = 1.0D0/alp(i,j,k)
+               ialp2 = ialp**2
+
 c              Calculate the stress energy tensor at this point
 c              ------------------------------------------------
 
@@ -120,9 +124,35 @@ c              --------------------------------
 c              Calculate the hamiltonian constraint
 c              ------------------------------------
 
+c              Geometric piece.
+
 #include "CactusEinstein/Einstein/src/macro/HAMADM_guts.h"
 
-               ham(i,j,k) = HAMADM_HAMADM - 16D0*Pi*Ttt/alp(i,j,k)**2
+c              Add matter terms: - 16*pi*rho
+c
+c              with rho defined as:
+c
+c              rho = n_a n_b T^{ab} = n^a n^b T_{ab}
+c                  = (T_00 - 2 beta^i T_{i0} + beta^i beta^j T_{ij})/alpha^2
+
+               m_rho = ialp2*Ttt
+
+               if (shift_state == SHIFT_ACTIVE) then
+
+                  m_rho = m_rho + ialp2
+     &               *(betax(i,j,k)**2*Txx
+     &               + betay(i,j,k)**2*Tyy
+     &               + betaz(i,j,k)**2*Tzz
+     &               +(betax(i,j,k)*betay(i,j,k)*Txy
+     &               + betax(i,j,k)*betaz(i,j,k)*Txz
+     &               + betay(i,j,k)*betaz(i,j,k)*Tyz)*2.0D0
+     &               -(betax(i,j,k)**2*Ttx
+     &               + betay(i,j,k)**2*Tty
+     &               + betaz(i,j,k)**2*Ttz)*2.0D0)
+
+               end if
+
+               ham(i,j,k) = HAMADM_HAMADM - 16.0D0*pi*m_rho
 
 #ifdef THORN_MAHC
                if ((some_hydro .eq. MAHCHYDRO).or.
@@ -136,14 +166,45 @@ c              ------------------------------------
 c              Calculate the Momentum constraints
 c              ----------------------------------
 
+c              Geometric piece.
+
 #include "CactusEinstein/Einstein/src/macro/MOMXADM_guts.h"
 #include "CactusEinstein/Einstein/src/macro/MOMYADM_guts.h"
 #include "CactusEinstein/Einstein/src/macro/MOMZADM_guts.h"
 
-               momx(i,j,k) = MOMXADM_MOMXADM - 8D0*Pi*Ttx/alp(i,j,k)**2
-               momy(i,j,k) = MOMYADM_MOMYADM - 8D0*Pi*Tty/alp(i,j,k)**2
-               momz(i,j,k) = MOMZADM_MOMZADM - 8D0*Pi*Ttz/alp(i,j,k)**2
+c              Add matter terms: - 8*pi*S_i
+c
+c              with S_i defined as:
+c
+c              S_i = - g_{ia} n_b T^{ab} = - g_i^a n^b T_{ab} 
+c                  = - (T_{i0} - beta^j T_{ij})/alpha
+
+               m_sx = - ialp*Ttx
+               m_sy = - ialp*Tty
+               m_sz = - ialp*Ttz
+
+               if (shift_state == SHIFT_ACTIVE) then
+
+                  m_sx = m_sx + ialp
+     &                 *(betax(i,j,k)*Txx
+     &                 + betay(i,j,k)*Txy
+     &                 + betaz(i,j,k)*Txz)
 
+                  m_sy = m_sy + ialp
+     &                 *(betax(i,j,k)*Txy
+     &                 + betay(i,j,k)*Tyy
+     &                 + betaz(i,j,k)*Tyz)
+
+                  m_sz = m_sz + ialp
+     &                 *(betax(i,j,k)*Txz
+     &                 + betay(i,j,k)*Tyz
+     &                 + betaz(i,j,k)*Tzz)
+
+               end if
+
+               momx(i,j,k) = MOMXADM_MOMXADM - 8.0D0*pi*m_sx
+               momy(i,j,k) = MOMYADM_MOMYADM - 8.0D0*pi*m_sy
+               momz(i,j,k) = MOMZADM_MOMZADM - 8.0D0*pi*m_sz
 
 #ifdef THORN_MAHC
                if ((some_hydro .eq. MAHCHYDRO).or.
@@ -168,7 +229,6 @@ c              ----------------------------------
 #include "CactusEinstein/Einstein/src/macro/MOMZADM_undefine.h"
 #include "CactusEinstein/Einstein/src/macro/UPPERMET_undefine.h"
 
-
 c     Synchronize and apply flat boundary conditions
 
       call CartSymGN(ierr,cctkGH,"admconstraints::admconstraints")
@@ -181,5 +241,7 @@ c     Synchronize and apply flat boundary conditions
          call CCTK_SyncGroup(cctkGH,"admconstraints::admconstraints")
       end if
 
+c     End
+
       end subroutine ADMConstraints