From f05be424d4813bf18dfe65b6ecfdb795ddef70e1 Mon Sep 17 00:00:00 2001 From: miguel Date: Fri, 16 Feb 2001 11:10:12 +0000 Subject: Adding shift terms. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/ADMConstraints/trunk@42 b7a48df3-cbbf-4440-997f-b4b717c9f7fc --- src/ADMConstraints.F | 84 +++++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 73 insertions(+), 11 deletions(-) (limited to 'src') 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 -- cgit v1.2.3