/*@@ @file adjADMConstr.F @date May 2002 @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 j_i @enddesc @history @hdate May 2002 @hauthor Hisaaki Shinkai @hdesc copied from ADMConstraints.F in EinsteinBase @endhistory @version $Header$ @@*/ c----------------------------------------------------------------------- c this routine is for adjusted ADM version, to evaluate constraints c during the evolution. c c coded by Hisaaki Shinkai c c----------------------------------------------------------------------- #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_DefineThorn.h" #ifdef BETATHORNS_CARTOON2D #include "BetaThorns/Cartoon2D/src/Cartoon2D_tensors.h" #endif subroutine adjADMConstraints(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. #include "EinsteinBase/ADMMacros/src/macro/ADM_Spacing_declare.h" 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 "EinsteinBase/ADMMacros/src/macro/HAMADM_declare.h" #include "EinsteinBase/ADMMacros/src/macro/MOMXADM_declare.h" #include "EinsteinBase/ADMMacros/src/macro/MOMYADM_declare.h" #include "EinsteinBase/ADMMacros/src/macro/MOMZADM_declare.h" #include "EinsteinBase/ADMMacros/src/macro/DETG_declare.h" #include "EinsteinBase/ADMMacros/src/macro/UPPERMET_declare.h" c -------------------------------------------------------------- c Grid parameters. #include "EinsteinBase/ADMMacros/src/macro/ADM_Spacing.h" nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) c Fill with zeros. aham = 0.0D0 amomx = 0.0D0 amomy = 0.0D0 amomz = 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 "EinsteinBase/ADMMacros/src/macro/DETG_guts.h" det = DETG_DETCG #include "EinsteinBase/ADMMacros/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 "EinsteinBase/ADMMacros/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 .eq. 1) 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 aham(i,j,k) = HAMADM_HAMADM - 16.0D0*pi*m_rho c Calculate the Momentum constraints c ---------------------------------- c Geometric piece. #include "EinsteinBase/ADMMacros/src/macro/MOMXADM_guts.h" #include "EinsteinBase/ADMMacros/src/macro/MOMYADM_guts.h" #include "EinsteinBase/ADMMacros/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 .eq. 1) 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 amomx(i,j,k) = MOMXADM_MOMXADM - 8.0D0*pi*m_sx amomy(i,j,k) = MOMYADM_MOMYADM - 8.0D0*pi*m_sy amomz(i,j,k) = MOMZADM_MOMZADM - 8.0D0*pi*m_sz end do end do end do #include "EinsteinBase/ADMMacros/src/macro/DETG_undefine.h" #include "EinsteinBase/ADMMacros/src/macro/UPPERMET_undefine.h" #include "EinsteinBase/ADMMacros/src/macro/HAMADM_undefine.h" #include "EinsteinBase/ADMMacros/src/macro/MOMXADM_undefine.h" #include "EinsteinBase/ADMMacros/src/macro/MOMYADM_undefine.h" #include "EinsteinBase/ADMMacros/src/macro/MOMZADM_undefine.h" c Apply symmetry boundary conditions. c call CartSymGN(ierr,cctkGH,"admconstraints::hamiltonian") c call CartSymGN(ierr,cctkGH,"admconstraints::momentum") call CartSymGN(ierr,cctkGH,"adm::adjADM_constraints") c Apply flat boundary conditions at outer boundaries. if (CCTK_Equals(bound,"flat") == 1) then c call BndFlatGN(ierr,cctkGH,sw,"admconstraints::hamiltonian") c call BndFlatGN(ierr,cctkGH,sw,"admconstraints::momentum") call BndFlatGN(ierr,cctkGH,sw,"adm::adjADM_constraints") end if call CCTK_SyncGroup(ierr,cctkGH,"adm::adjADM_constraints") c Synchronize. c check later c if (adjconstraint_communication.eq.1) then c call CCTK_SyncGroup(ierr,cctkGH,"admconstraints::hamiltonian") c call CCTK_SyncGroup(ierr,cctkGH,"admconstraints::momentum") c 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 adjADMConstraints