/*@@ @file ADMConstraints.F @date August 98 @desc Calculate the ADM Constraints for output: Hamiltonian Constraint is: H = R - K^i_j K^j_i + trK^2 - 16 Pi rho Momentum Constraints are: M_i = Del_j K_i^j - Del_i trK - 8 Pi S_i @enddesc @version $Header$ @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_DefineThorn.h" #include "CactusEinstein/Einstein/src/Einstein.h" #ifdef BETATHORNS_CARTOON2D #include "BetaThorns/Cartoon2D/src/Cartoon2D_tensors.h" #endif subroutine ADMConstraints(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS integer :: i,j,k integer :: nx,ny,nz #ifdef EXCISION_LEGOEXCISION CCTK_REAL, allocatable, dimension (:,:,:) :: dirx,diry,dirz,aux #endif c Stencil width used for calculating constraints c (for outer boundary condition) integer, dimension(3),parameter :: sw = 1 c Return code from Cactus sync routine and boundary conditions. integer ierr c Various real variables. CCTK_REAL :: dx,dy,dz CCTK_REAL :: m_rho,m_sx,m_sy,m_sz CCTK_REAL :: pi,ialp,ialp2 CCTK_REAL :: det,uxx,uyy,uzz,uxy,uxz,uyz c Temporaries for the Stress-Energy tensor. CCTK_REAL :: Ttt,Ttx,Tty,Ttz,Txx,Txy,Txz,Tyy,Tyz,Tzz c Matter declarations. #include "CalcTmunu_temps.inc" c Macros from Standard Einstein. #include "CactusEinstein/Einstein/src/macro/HAMADM_declare.h" #include "CactusEinstein/Einstein/src/macro/MOMXADM_declare.h" #include "CactusEinstein/Einstein/src/macro/MOMYADM_declare.h" #include "CactusEinstein/Einstein/src/macro/MOMZADM_declare.h" #include "CactusEinstein/Einstein/src/macro/DETG_declare.h" #include "CactusEinstein/Einstein/src/macro/UPPERMET_declare.h" c -------------------------------------------------------------- c Grid parameters. dx = cctk_delta_space(1) dy = cctk_delta_space(2) dz = cctk_delta_space(3) nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) c Fill with zeros. ham = 0.0D0 momx = 0.0D0 momy = 0.0D0 momz = 0.0D0 c Calculate constraints. pi = acos(-1.0D0) do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 ialp = 1.0D0/alp(i,j,k) ialp2 = ialp**2 c Calculate the stress energy tensor at this point c ------------------------------------------------ c This may be needed for CalcTmunu #include "CactusEinstein/Einstein/src/macro/DETG_guts.h" det = DETG_DETCG #include "CactusEinstein/Einstein/src/macro/UPPERMET_guts.h" uxx = UPPERMET_UXX; uxy = UPPERMET_UXY; uxz = UPPERMET_UXZ uyy = UPPERMET_UYY; uyz = UPPERMET_UYZ; uzz = UPPERMET_UZZ c Initialize stress-energy tensor components. Ttt = 0.0D0 Ttx = 0.0D0; Tty = 0.0D0; Ttz = 0.0D0 Txx = 0.0D0; Tyy = 0.0D0; Tzz = 0.0D0 Txy = 0.0D0; Txz = 0.0D0; Tyz = 0.0D0 c Include macro for stress energy tensor. #include "CalcTmunu.inc" c Calculate the hamiltonian constraint c ------------------------------------ c Geometric piece. #include "CactusEinstein/Einstein/src/macro/HAMADM_guts.h" 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)*Ttx & + betay(i,j,k)*Tty & + betaz(i,j,k)*Ttz)*2.0D0) end if ham(i,j,k) = HAMADM_HAMADM - 16.0D0*pi*m_rho 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" 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 end do end do end do #include "CactusEinstein/Einstein/src/macro/DETG_undefine.h" #include "CactusEinstein/Einstein/src/macro/UPPERMET_undefine.h" #include "CactusEinstein/Einstein/src/macro/HAMADM_undefine.h" #include "CactusEinstein/Einstein/src/macro/MOMXADM_undefine.h" #include "CactusEinstein/Einstein/src/macro/MOMYADM_undefine.h" #include "CactusEinstein/Einstein/src/macro/MOMZADM_undefine.h" c LegoExcision (must be done before symmetries are applied). if (excise==1) then #ifdef EXCISION_LEGOEXCISION allocate(dirx(nx,ny,nz),diry(nx,ny,nz),dirz(nx,ny,nz)) allocate(aux(nx,ny,nz)) aux = 0.0D0 call excision_findboundary(ierr,emask,nx,ny,nz) call CCTK_SyncGroup(ierr,cctkGH,"einstein::mask") call CartSymGN(ierr,cctkGH,"einstein::mask") call excision_findnormals(ierr,emask,dirx,diry,dirz,nx,ny,nz) call excision_extrapolate(ierr,ham ,aux,emask, . dirx,diry,dirz,nx,ny,nz,0.0D0) call excision_extrapolate(ierr,momx,aux,emask, . dirx,diry,dirz,nx,ny,nz,0.0D0) call excision_extrapolate(ierr,momy,aux,emask, . dirx,diry,dirz,nx,ny,nz,0.0D0) call excision_extrapolate(ierr,momz,aux,emask, . dirx,diry,dirz,nx,ny,nz,0.0D0) deallocate(dirx,diry,dirz) deallocate(aux) #else call CCTK_WARN(0,"You have not compiled with LegoExcision") #endif end if c Apply symmetry boundary conditions. call CartSymGN(ierr,cctkGH,"admconstraints::hamiltonian") call CartSymGN(ierr,cctkGH,"admconstraints::momentum") c Apply flat boundary conditions at outer boundaries. if (CCTK_Equals(bound,"flat") == 1) then call BndFlatGN(ierr,cctkGH,sw,"admconstraints::hamiltonian") call BndFlatGN(ierr,cctkGH,sw,"admconstraints::momentum") end if c Synchronize. if (constraint_communication.eq.1) then call CCTK_SyncGroup(ierr,cctkGH,"admconstraints::hamiltonian") call CCTK_SyncGroup(ierr,cctkGH,"admconstraints::momentum") end if c Cartoon. if (cartoon==1) then #ifdef BETATHORNS_CARTOON2D call BndCartoon2DGN(ierr,cctkGH,TENSORTYPE_SCALAR,"admconstraints::hamiltonian") call BndCartoon2DGN(ierr,cctkGH,TENSORTYPE_U,"admconstraints::momentum") #else call CCTK_WARN(0,"You have not compiled with Cartoon2D") #endif end if c End end subroutine ADMConstraints