/*@@ @file ADM_Constraints.F @date Converted to cactus4 in March 2000. @author Mark Miller @desc Calculates the constraints in 3+1 general relativity. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "CactusEinstein/Einstein/src/Einstein.h" /*@@ @routine ADM_Constraints @date Converted to cactus4 in March 2000. @author Mark Miller @desc Calculates the constraints in 3+1 general relativity. @enddesc @calls @calledby @history @endhistory @@*/ subroutine ADM_Constraints(CCTK_FARGUMENTS) implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS INTEGER CCTK_Equals c Local Vars integer i,j,k,nx,ny,nz,ierr integer, dimension(3),parameter :: sw = 1 logical wilson_style_mom CCTK_REAL dx,dy,dz CCTK_REAL det,uxx,uxy,uxz,uyy,uyz,uzz CCTK_REAL scalar_curv,two,pi,half #include "EUAstro_SpaceTime/ADM_Constraints/src/deriv_temps_con.inc" #include "EUAstro_SpaceTime/ADM_Constraints/src/ricci_temps.inc" CCTK_REAL my_trk(cctk_lsh(1),cctk_lsh(2),3) CCTK_REAL kkxx(cctk_lsh(1),cctk_lsh(2),3) CCTK_REAL kkxy(cctk_lsh(1),cctk_lsh(2),3) CCTK_REAL kkxz(cctk_lsh(1),cctk_lsh(2),3) CCTK_REAL kkyx(cctk_lsh(1),cctk_lsh(2),3) CCTK_REAL kkyy(cctk_lsh(1),cctk_lsh(2),3) CCTK_REAL kkyz(cctk_lsh(1),cctk_lsh(2),3) CCTK_REAL kkzx(cctk_lsh(1),cctk_lsh(2),3) CCTK_REAL kkzy(cctk_lsh(1),cctk_lsh(2),3) CCTK_REAL kkzz(cctk_lsh(1),cctk_lsh(2),3) CCTK_REAL dxxx,dxxy,dxyx,dxxz,dxzx,dxyy,dxyz,dxzy,dxzz CCTK_REAL dyxx,dyxy,dyyx,dyxz,dyzx,dyyy,dyyz,dyzy,dyzz CCTK_REAL dzxx,dzxy,dzyx,dzxz,dzzx,dzyy,dzyz,dzzy,dzzz CCTK_REAL Gxxx,Gxxy,Gxxz,Gxyy,Gxyz,Gxzz CCTK_REAL Gyxx,Gyxy,Gyxz,Gyyy,Gyyz,Gyzz CCTK_REAL Gzxx,Gzxy,Gzxz,Gzyy,Gzyz,Gzzz CCTK_REAL Ttt,Ttx,Tty,Ttz,Txx,Txy,Txz,Tyy,Tyz,Tzz 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) wilson_style_mom = (admcons_wilson_momentum .eq. 1) ham = 0.0d0 momx = 0.0d0 momy = 0.0d0 momz = 0.0d0 two = 2.0d0 pi = 3.14159265358979d0 half = 0.5d0 c calculate k=1 and k=2 of non-local variables do k=1,2 do j=1,ny do i=1,nx det = gxx(i,j,k)*gyy(i,j,k)*gzz(i,j,k) + . 2.*gxy(i,j,k)*gxz(i,j,k)*gyz(i,j,k) - . gxx(i,j,k) * (gyz(i,j,k)**2) - . gyy(i,j,k) * (gxz(i,j,k)**2) - . gzz(i,j,k) * (gxy(i,j,k)**2) uxx=(-gyz(i,j,k)**2 + gyy(i,j,k)*gzz(i,j,k))/det uxy=(gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k))/det uyy=(-gxz(i,j,k)**2 + gxx(i,j,k)*gzz(i,j,k))/det uxz=(-gxz(i,j,k)*gyy(i,j,k) + gxy(i,j,k)*gyz(i,j,k))/det uyz=(gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k))/det uzz=(-gxy(i,j,k)**2 + gxx(i,j,k)*gyy(i,j,k))/det my_trk(i,j,k+1) = kxx(i,j,k) * uxx + kyy(i,j,k) * uyy + . kzz(i,j,k) * uzz + two * kxy(i,j,k) * uxy + . two * kxz(i,j,k) * uxz + two * kyz(i,j,k) * uyz c calculate kkab = K^{a}_{b} kkxx(i,j,k+1) = uxx*kxx(i,j,k) + uxy*kxy(i,j,k) + uxz*kxz(i,j,k) kkxy(i,j,k+1) = uxx*kxy(i,j,k) + uxy*kyy(i,j,k) + uxz*kyz(i,j,k) kkxz(i,j,k+1) = uxx*kxz(i,j,k) + uxy*kyz(i,j,k) + uxz*kzz(i,j,k) kkyx(i,j,k+1) = uxy*kxx(i,j,k) + uyy*kxy(i,j,k) + uyz*kxz(i,j,k) kkyy(i,j,k+1) = uxy*kxy(i,j,k) + uyy*kyy(i,j,k) + uyz*kyz(i,j,k) kkyz(i,j,k+1) = uxy*kxz(i,j,k) + uyy*kyz(i,j,k) + uyz*kzz(i,j,k) kkzx(i,j,k+1) = uxz*kxx(i,j,k) + uyz*kxy(i,j,k) + uzz*kxz(i,j,k) kkzy(i,j,k+1) = uxz*kxy(i,j,k) + uyz*kyy(i,j,k) + uzz*kyz(i,j,k) kkzz(i,j,k+1) = uxz*kxz(i,j,k) + uyz*kyz(i,j,k) + uzz*kzz(i,j,k) enddo enddo enddo do k=2,nz-1 c first, move non-locals and calculate new k=3 row do j=1,ny do i=1,nx my_trk(i,j,1) = my_trk(i,j,2) kkxx(i,j,1) = kkxx(i,j,2) kkxy(i,j,1) = kkxy(i,j,2) kkxz(i,j,1) = kkxz(i,j,2) kkyx(i,j,1) = kkyx(i,j,2) kkyy(i,j,1) = kkyy(i,j,2) kkyz(i,j,1) = kkyz(i,j,2) kkzx(i,j,1) = kkzx(i,j,2) kkzy(i,j,1) = kkzy(i,j,2) kkzz(i,j,1) = kkzz(i,j,2) my_trk(i,j,2) = my_trk(i,j,3) kkxx(i,j,2) = kkxx(i,j,3) kkxy(i,j,2) = kkxy(i,j,3) kkxz(i,j,2) = kkxz(i,j,3) kkyx(i,j,2) = kkyx(i,j,3) kkyy(i,j,2) = kkyy(i,j,3) kkyz(i,j,2) = kkyz(i,j,3) kkzx(i,j,2) = kkzx(i,j,3) kkzy(i,j,2) = kkzy(i,j,3) kkzz(i,j,2) = kkzz(i,j,3) c calculate new row det = gxx(i,j,k+1)*gyy(i,j,k+1)*gzz(i,j,k+1) + . 2.0d0*gxy(i,j,k+1)*gxz(i,j,k+1)*gyz(i,j,k+1) - . gxx(i,j,k+1) * (gyz(i,j,k+1)**2) - . gyy(i,j,k+1) * (gxz(i,j,k+1)**2) - . gzz(i,j,k+1) * (gxy(i,j,k+1)**2) uxx=(-gyz(i,j,k+1)**2 + gyy(i,j,k+1)*gzz(i,j,k+1))/det uxy=(gxz(i,j,k+1)*gyz(i,j,k+1) - gxy(i,j,k+1)*gzz(i,j,k+1))/det uyy=(-gxz(i,j,k+1)**2 + gxx(i,j,k+1)*gzz(i,j,k+1))/det uxz=(-gxz(i,j,k+1)*gyy(i,j,k+1) + gxy(i,j,k+1)*gyz(i,j,k+1))/det uyz=(gxy(i,j,k+1)*gxz(i,j,k+1) - gxx(i,j,k+1)*gyz(i,j,k+1))/det uzz=(-gxy(i,j,k+1)**2 + gxx(i,j,k+1)*gyy(i,j,k+1))/det my_trk(i,j,3) = kxx(i,j,k+1) * uxx + kyy(i,j,k+1) * uyy + . kzz(i,j,k+1) * uzz + two * kxy(i,j,k+1) * uxy + . two * kxz(i,j,k+1) * uxz + two * kyz(i,j,k+1) * uyz c calculate kkab = K^{a}_{b} kkxx(i,j,3) = uxx*kxx(i,j,k+1) + uxy*kxy(i,j,k+1) + uxz*kxz(i,j,k+1) kkxy(i,j,3) = uxx*kxy(i,j,k+1) + uxy*kyy(i,j,k+1) + uxz*kyz(i,j,k+1) kkxz(i,j,3) = uxx*kxz(i,j,k+1) + uxy*kyz(i,j,k+1) + uxz*kzz(i,j,k+1) kkyx(i,j,3) = uxy*kxx(i,j,k+1) + uyy*kxy(i,j,k+1) + uyz*kxz(i,j,k+1) kkyy(i,j,3) = uxy*kxy(i,j,k+1) + uyy*kyy(i,j,k+1) + uyz*kyz(i,j,k+1) kkyz(i,j,3) = uxy*kxz(i,j,k+1) + uyy*kyz(i,j,k+1) + uyz*kzz(i,j,k+1) kkzx(i,j,3) = uxz*kxx(i,j,k+1) + uyz*kxy(i,j,k+1) + uzz*kxz(i,j,k+1) kkzy(i,j,3) = uxz*kxy(i,j,k+1) + uyz*kyy(i,j,k+1) + uzz*kyz(i,j,k+1) kkzz(i,j,3) = uxz*kxz(i,j,k+1) + uyz*kyz(i,j,k+1) + uzz*kzz(i,j,k+1) enddo enddo c calculate constraints do j=2,ny-1 do i=2,nx-1 det= -(gxz(i,j,k)**2*gyy(i,j,k)) + . 2.0d0*gxy(i,j,k)*gxz(i,j,k)*gyz(i,j,k) - . gxx(i,j,k)*gyz(i,j,k)**2 - . gxy(i,j,k)**2*gzz(i,j,k) + . gxx(i,j,k)*gyy(i,j,k)*gzz(i,j,k) uxx=(-gyz(i,j,k)**2 + gyy(i,j,k)*gzz(i,j,k))/det uxy=(gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k))/det uyy=(-gxz(i,j,k)**2 + gxx(i,j,k)*gzz(i,j,k))/det uxz=(-gxz(i,j,k)*gyy(i,j,k) + gxy(i,j,k)*gyz(i,j,k))/det uyz=(gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k))/det uzz=(-gxy(i,j,k)**2 + gxx(i,j,k)*gyy(i,j,k))/det #include "EUAstro_SpaceTime/ADM_Constraints/src/setup_derivs_con.inc" #include "EUAstro_SpaceTime/ADM_Constraints/src/ricci_calc.inc" c calculate scalar curvature scalar_curv = riccixx*uxx + ricciyy*uyy + riccizz*uzz + . two*riccixy*uxy + two*riccixz*uxz + two*ricciyz*uyz c Hamiltonian constraint ham(i,j,k) = scalar_curv + (my_trk(i,j,2))**2 - . (kkxx(i,j,2)*kkxx(i,j,2) + kkxy(i,j,2)*kkyx(i,j,2) + . kkxz(i,j,2)*kkzx(i,j,2) + kkyx(i,j,2)*kkxy(i,j,2) + . kkyy(i,j,2)*kkyy(i,j,2) + kkyz(i,j,2)*kkzy(i,j,2) + . kkzx(i,j,2)*kkxz(i,j,2) + kkzy(i,j,2)*kkyz(i,j,2) + . kkzz(i,j,2)*kkzz(i,j,2)) c now, calculate momentum constraints c these d variables correspond to BM D_{i}_{j}^{k} variables c (thus, the factor of 2) dxxx = (dexgxx*uxx + dexgxy*uxy + dexgxz*uxz)/2.0d0 dxxy = (dexgxx*uxy + dexgxy*uyy + dexgxz*uyz)/2.0d0 dxxz = (dexgxx*uxz + dexgxy*uyz + dexgxz*uzz)/2.0d0 dxyx = (dexgxy*uxx + dexgyy*uxy + dexgyz*uxz)/2.0d0 dxyy = (dexgxy*uxy + dexgyy*uyy + dexgyz*uyz)/2.0d0 dxyz = (dexgxy*uxz + dexgyy*uyz + dexgyz*uzz)/2.0d0 dxzx = (dexgxz*uxx + dexgyz*uxy + dexgzz*uxz)/2.0d0 dxzy = (dexgxz*uxy + dexgyz*uyy + dexgzz*uyz)/2.0d0 dxzz = (dexgxz*uxz + dexgyz*uyz + dexgzz*uzz)/2.0d0 dyxx = (deygxx*uxx + deygxy*uxy + deygxz*uxz)/2.0d0 dyxy = (deygxx*uxy + deygxy*uyy + deygxz*uyz)/2.0d0 dyxz = (deygxx*uxz + deygxy*uyz + deygxz*uzz)/2.0d0 dyyx = (deygxy*uxx + deygyy*uxy + deygyz*uxz)/2.0d0 dyyy = (deygxy*uxy + deygyy*uyy + deygyz*uyz)/2.0d0 dyyz = (deygxy*uxz + deygyy*uyz + deygyz*uzz)/2.0d0 dyzx = (deygxz*uxx + deygyz*uxy + deygzz*uxz)/2.0d0 dyzy = (deygxz*uxy + deygyz*uyy + deygzz*uyz)/2.0d0 dyzz = (deygxz*uxz + deygyz*uyz + deygzz*uzz)/2.0d0 dzxx = (dezgxx*uxx + dezgxy*uxy + dezgxz*uxz)/2.0d0 dzxy = (dezgxx*uxy + dezgxy*uyy + dezgxz*uyz)/2.0d0 dzxz = (dezgxx*uxz + dezgxy*uyz + dezgxz*uzz)/2.0d0 dzyx = (dezgxy*uxx + dezgyy*uxy + dezgyz*uxz)/2.0d0 dzyy = (dezgxy*uxy + dezgyy*uyy + dezgyz*uyz)/2.0d0 dzyz = (dezgxy*uxz + dezgyy*uyz + dezgyz*uzz)/2.0d0 dzzx = (dezgxz*uxx + dezgyz*uxy + dezgzz*uxz)/2.0d0 dzzy = (dezgxz*uxy + dezgyz*uyy + dezgzz*uyz)/2.0d0 dzzz = (dezgxz*uxz + dezgyz*uyz + dezgzz*uzz)/2.0d0 Gxxx = two*dxxx - dexgxx*half*uxx - deygxx*half*uxy - . dezgxx*half*uxz Gxxy = dxyx + dyxx - dexgxy*half*uxx - deygxy*half*uxy - . dezgxy*half*uxz Gxxz = dxzx + dzxx - dexgxz*half*uxx - deygxz*half*uxy - . dezgxz*half*uxz Gxyy = two*dyyx - dexgyy*half*uxx - deygyy*half*uxy - . dezgyy*half*uxz Gxyz = dyzx + dzyx - dexgyz*half*uxx - deygyz*half*uxy - . dezgyz*half*uxz Gxzz = two*dzzx - dexgzz*half*uxx - deygzz*half*uxy - . dezgzz*half*uxz Gyxx = two*dxxy - dexgxx*half*uxy - deygxx*half*uyy - . dezgxx*half*uyz Gyxy = dxyy + dyxy - dexgxy*half*uxy - deygxy*half*uyy - . dezgxy*half*uyz Gyxz = dxzy + dzxy - dexgxz*half*uxy - deygxz*half*uyy - . dezgxz*half*uyz Gyyy = two*dyyy - dexgyy*half*uxy - deygyy*half*uyy - . dezgyy*half*uyz Gyyz = dyzy + dzyy - dexgyz*half*uxy - deygyz*half*uyy - . dezgyz*half*uyz Gyzz = two*dzzy - dexgzz*half*uxy - deygzz*half*uyy - . dezgzz*half*uyz Gzxx = two*dxxz - dexgxx*half*uxz - deygxx*half*uyz - . dezgxx*half*uzz Gzxy = dxyz + dyxz - dexgxy*half*uxz - deygxy*half*uyz - . dezgxy*half*uzz Gzxz = dxzz + dzxz - dexgxz*half*uxz - deygxz*half*uyz - . dezgxz*half*uzz Gzyy = two*dyyz - dexgyy*half*uxz - deygyy*half*uyz - . dezgyy*half*uzz Gzyz = dyzz + dzyz - dexgyz*half*uxz - deygyz*half*uyz - . dezgyz*half*uzz Gzzz = two*dzzz - dexgzz*half*uxz - deygzz*half*uyz - . dezgzz*half*uzz momx(i,j,k) = (kkxx(i+1,j,2) - kkxx(i-1,j,2))/(two*dx) + . (kkyx(i,j+1,2) - kkyx(i,j-1,2))/(two*dy) + . (kkzx(i,j,3) - kkzx(i,j,1))/(two*dz) + . Gxxx*kkxx(i,j,2) + Gyxy*kkxx(i,j,2) + Gzxz*kkxx(i,j,2) + . Gxxy*kkyx(i,j,2) + Gyyy*kkyx(i,j,2) + Gzyz*kkyx(i,j,2) + . Gxxz*kkzx(i,j,2) + Gyyz*kkzx(i,j,2) + Gzzz*kkzx(i,j,2) - . Gxxx*kkxx(i,j,2) - Gxxy*kkyx(i,j,2) - Gxxz*kkzx(i,j,2) - . Gyxx*kkxy(i,j,2) - Gyxy*kkyy(i,j,2) - Gyxz*kkzy(i,j,2) - . Gzxx*kkxz(i,j,2) - Gzxy*kkyz(i,j,2) - Gzxz*kkzz(i,j,2) - . (my_trk(i+1,j,2) - my_trk(i-1,j,2))/(two*dx) momy(i,j,k) = (kkxy(i+1,j,2) - kkxy(i-1,j,2))/(two*dx) + . (kkyy(i,j+1,2) - kkyy(i,j-1,2))/(two*dy) + . (kkzy(i,j,3) - kkzy(i,j,1))/(two*dz) + . Gxxx*kkxy(i,j,2) + Gyxy*kkxy(i,j,2) + Gzxz*kkxy(i,j,2) + . Gxxy*kkyy(i,j,2) + Gyyy*kkyy(i,j,2) + Gzyz*kkyy(i,j,2) + . Gxxz*kkzy(i,j,2) + Gyyz*kkzy(i,j,2) + Gzzz*kkzy(i,j,2) - . Gxxy*kkxx(i,j,2) - Gxyy*kkyx(i,j,2) - Gxyz*kkzx(i,j,2) - . Gyxy*kkxy(i,j,2) - Gyyy*kkyy(i,j,2) - Gyyz*kkzy(i,j,2) - . Gzxy*kkxz(i,j,2) - Gzyy*kkyz(i,j,2) - Gzyz*kkzz(i,j,2) - . (my_trk(i,j+1,2) - my_trk(i,j-1,2))/(two*dy) momz(i,j,k) = (kkxz(i+1,j,2) - kkxz(i-1,j,2))/(two*dx) + . (kkyz(i,j+1,2) - kkyz(i,j-1,2))/(two*dy) + . (kkzz(i,j,3) - kkzz(i,j,1))/(two*dz) + . Gxxx*kkxz(i,j,2) + Gyxy*kkxz(i,j,2) + Gzxz*kkxz(i,j,2) + . Gxxy*kkyz(i,j,2) + Gyyy*kkyz(i,j,2) + Gzyz*kkyz(i,j,2) + . Gxxz*kkzz(i,j,2) + Gyyz*kkzz(i,j,2) + Gzzz*kkzz(i,j,2) - . Gxxz*kkxx(i,j,2) - Gxyz*kkyx(i,j,2) - Gxzz*kkzx(i,j,2) - . Gyxz*kkxy(i,j,2) - Gyyz*kkyy(i,j,2) - Gyzz*kkzy(i,j,2) - . Gzxz*kkxz(i,j,2) - Gzyz*kkyz(i,j,2) - Gzzz*kkzz(i,j,2) - . (my_trk(i,j,3) - my_trk(i,j,1))/(two*dz) if(mahc_evolution.eq.1) then ham(i,j,k) = ham(i,j,k) - 16.0d0 * pi * . ( (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k)) * . w_lorentz(i,j,k)**2 - press(i,j,k)) momx(i,j,k) = momx(i,j,k) - . 8.0d0 * pi * sx(i,j,k) / sqrt(det) momy(i,j,k) = momy(i,j,k) - . 8.0d0 * pi * sy(i,j,k) / sqrt(det) momz(i,j,k) = momz(i,j,k) - . 8.0d0 * pi * sz(i,j,k) / sqrt(det) if(wilson_style_mom) then if(shift_state .eq. SHIFT_ACTIVE) then momx(i,j,k) = momx(i,j,k) + 8.0d0 * pi * . w_lorentz(i,j,k)**2 * (rho(i,j,k)*(1.0d0 + eps(i,j,k)) + . press(i,j,k)) * (1.0d0/alp(i,j,k)) * . (gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + . gxz(i,j,k)*betaz(i,j,k)) momy(i,j,k) = momy(i,j,k) + 8.0d0 * pi * . w_lorentz(i,j,k)**2 * (rho(i,j,k)*(1.0d0 + eps(i,j,k)) + . press(i,j,k)) * (1.0d0/alp(i,j,k)) * . (gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + . gyz(i,j,k)*betaz(i,j,k)) momz(i,j,k) = momz(i,j,k) + 8.0d0 * pi * . w_lorentz(i,j,k)**2 * (rho(i,j,k)*(1.0d0 + eps(i,j,k)) + . press(i,j,k)) * (1.0d0/alp(i,j,k)) * . (gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + . gzz(i,j,k)*betaz(i,j,k)) endif endif endif enddo enddo enddo call CartSymGN(ierr,cctkGH,"adm_constraints::admconstraints") if (CCTK_Equals(bound,"flat") == 1) then call BndFlatGN(ierr,cctkGH,sw,"adm_constraints::admconstraints") end if if(constraint_communication.eq.1) then call CCTK_SyncGroup(cctkGH,"adm_constraints::ADMconstraints") endif end subroutine ADM_Constraints