/*@@ @file ADMConstraints.F77 @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_Arguments.h" #include "cctk_Functions.h" #include "cctk_Parameters.h" subroutine ADMConstraints(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS integer i,j,k integer nx,ny,nz c Stencil width used for calculating constraints c (for outer boundary condition) integer sw(3) logical docalc c Return code from Cactus sync routine and boundary conditions. integer ierr c Various real variables. CCTK_REAL m_rho,m_sx,m_sy,m_sz CCTK_REAL pi,ialp,ialp2 CCTK_REAL det,uxx,uyy,uzz,uxy,uxz,uyz #include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing_declare.h" 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/ADMMacros/src/macro/HAMADM_declare.h" #include "CactusEinstein/ADMMacros/src/macro/MOMXADM_declare.h" #include "CactusEinstein/ADMMacros/src/macro/MOMYADM_declare.h" #include "CactusEinstein/ADMMacros/src/macro/MOMZADM_declare.h" #include "CactusEinstein/ADMMacros/src/macro/DETG_declare.h" #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h" c -------------------------------------------------------------- sw(1) = 1 sw(2) = 1 sw(3) = 1 c Grid parameters. #include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing.h" nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) c Fill with zeros. do k=1,nz do j=1,ny do i=1,nx ham(i,j,k) = 0.0D0 hamnormalized(i,j,k) = 0.0D0 momx(i,j,k) = 0.0D0 momy(i,j,k) = 0.0D0 momz(i,j,k) = 0.0D0 end do end do end do c Calculate constraints. pi = acos(-1.0D0) do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 if ( (i<3).or.(i>cctk_lsh(1)-2).or. . (j<3).or.(j>cctk_lsh(2)-2).or. . (k<3).or.(k>cctk_lsh(3)-2) ) then local_spatial_order = 2 else local_spatial_order = spatial_order end if docalc = .TRUE. if (use_mask .eq. 1) then if (abs(emask(i,j,k)-1) > 0.001) then docalc = .FALSE. end if end if if (docalc) then 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/ADMMacros/src/macro/DETG_guts.h" det = DETG_DETCG #include "CactusEinstein/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 "CactusEinstein/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 ham(i,j,k) = HAMADM_HAMADM - 16.0D0*pi*m_rho if ((HAMADM_HAMADMABS + abs(16.0D0*pi*m_rho))==0) then hamnormalized(i,j,k) = abs(ham(i,j,k)) else hamnormalized(i,j,k) = abs(ham(i,j,k))/ & (HAMADM_HAMADMABS + abs(16.0D0*pi*m_rho)) end if c Calculate the Momentum constraints c ---------------------------------- c Geometric piece. #include "CactusEinstein/ADMMacros/src/macro/MOMXADM_guts.h" #include "CactusEinstein/ADMMacros/src/macro/MOMYADM_guts.h" #include "CactusEinstein/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 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 else ham(i,j,k) = excised_value momx(i,j,k) = excised_value momy(i,j,k) = excised_value momz(i,j,k) = excised_value end if end do end do end do #include "CactusEinstein/ADMMacros/src/macro/DETG_undefine.h" #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_undefine.h" #include "CactusEinstein/ADMMacros/src/macro/HAMADM_undefine.h" #include "CactusEinstein/ADMMacros/src/macro/MOMXADM_undefine.h" #include "CactusEinstein/ADMMacros/src/macro/MOMYADM_undefine.h" #include "CactusEinstein/ADMMacros/src/macro/MOMZADM_undefine.h" c Apply symmetry boundary conditions. call CartSymGN(ierr,cctkGH,"admconstraints::hamiltonian") call CartSymGN(ierr,cctkGH,"admconstraints::normalized_hamiltonian") call CartSymGN(ierr,cctkGH,"admconstraints::momentum") c Apply flat boundary conditions at outer boundaries. if (CCTK_Equals(bound,"flat") == 1) then ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, $ "admconstraints::hamiltonian", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, $ "admconstraints::normalized_hamiltonian", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, $ "admconstraints::momentum", "Flat") else ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, $ "admconstraints::hamiltonian", "None") ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, $ "admconstraints::normalized_hamiltonian", "None") ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, $ "admconstraints::momentum", "None") end if c Synchronize. if (constraint_communication.eq.1) then call CCTK_SyncGroup(ierr,cctkGH,"admconstraints::hamiltonian") call CCTK_SyncGroup(ierr,cctkGH,"admconstraints::normalized_hamiltonian") call CCTK_SyncGroup(ierr,cctkGH,"admconstraints::momentum") end if c Cartoon. if (cartoon==1) then if (CCTK_IsThornActive("Cartoon2D")==1) then ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 0, -1, $ "admconstraints::ham", "None") ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 0, -1, $ "admconstraints::hamnormalized", "None") ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 0, -1, $ "admconstraints::momentum", "None") else call CCTK_WARN(0,"You have not activated Cartoon2D") end if end if return end