cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccc ADD_DIV_ASTAR_VEC_SOURCE cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccc This routine adds in the divergence of the astar tensor ccc to the source for the momentum equation. Usually, ccc Yorks formulation assumes astar is divergence free, but ccc if it is not, we have to account for it in the momentum ccc equation. (M.M. 5-27-98) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc #include "cctk.h" #include "cctk_Arguments.h" subroutine add_div_astar_vec_source(CCTK_FARGUMENTS) implicit none DECLARE_CCTK_FARGUMENTS c local vars integer i,j,k CCTK_REAL delxgxx,delxgxy,delxgxz,delxgyy,delxgyz,delxgzz CCTK_REAL delygxx,delygxy,delygxz,delygyy,delygyz,delygzz CCTK_REAL delzgxx,delzgxy,delzgxz,delzgyy,delzgyz,delzgzz CCTK_REAL astarxx,astarxy,astarxz,astaryy,astaryz,astarzz CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz, dx, dy, dz c dx settings dx = cctk_delta_space(1) dy = cctk_delta_space(2) dz = cctk_delta_space(3) do k=2,cctk_lsh(3)-1 do j=2,cctk_lsh(2)-1 do i=2,cctk_lsh(1)-1 c first add in the "divergence of a vector" part ivp_wx_source(i,j,k) = ivp_wx_source(i,j,k) + . ( (sqrt(ivp_det(i+1,j,k))*ivp_astarxx(i+1,j,k) - . sqrt(ivp_det(i-1,j,k))*ivp_astarxx(i-1,j,k))/(2.0d0*dx) + . (sqrt(ivp_det(i,j+1,k))*ivp_astarxy(i,j+1,k) - . sqrt(ivp_det(i,j-1,k))*ivp_astarxy(i,j-1,k))/(2.0d0*dy) + . (sqrt(ivp_det(i,j,k+1))*ivp_astarxz(i,j,k+1) - . sqrt(ivp_det(i,j,k-1))*ivp_astarxz(i,j,k-1))/(2.0d0*dz))/ . sqrt(ivp_det(i,j,k)) ivp_wy_source(i,j,k) = ivp_wy_source(i,j,k) + . ( (sqrt(ivp_det(i+1,j,k))*ivp_astarxy(i+1,j,k) - . sqrt(ivp_det(i-1,j,k))*ivp_astarxy(i-1,j,k))/(2.0d0*dx) + . (sqrt(ivp_det(i,j+1,k))*ivp_astaryy(i,j+1,k) - . sqrt(ivp_det(i,j-1,k))*ivp_astaryy(i,j-1,k))/(2.0d0*dy) + . (sqrt(ivp_det(i,j,k+1))*ivp_astaryz(i,j,k+1) - . sqrt(ivp_det(i,j,k-1))*ivp_astaryz(i,j,k-1))/(2.0d0*dz))/ . sqrt(ivp_det(i,j,k)) ivp_wz_source(i,j,k) = ivp_wz_source(i,j,k) + . ( (sqrt(ivp_det(i+1,j,k))*ivp_astarxz(i+1,j,k) - . sqrt(ivp_det(i-1,j,k))*ivp_astarxz(i-1,j,k))/(2.0d0*dx) + . (sqrt(ivp_det(i,j+1,k))*ivp_astaryz(i,j+1,k) - . sqrt(ivp_det(i,j-1,k))*ivp_astaryz(i,j-1,k))/(2.0d0*dy) + . (sqrt(ivp_det(i,j,k+1))*ivp_astarzz(i,j,k+1) - . sqrt(ivp_det(i,j,k-1))*ivp_astarzz(i,j,k-1))/(2.0d0*dz))/ . sqrt(ivp_det(i,j,k)) c now take care of that other index! astarxx = ivp_astarxx(i,j,k) astarxy = ivp_astarxy(i,j,k) astarxz = ivp_astarxz(i,j,k) astaryy = ivp_astaryy(i,j,k) astaryz = ivp_astaryz(i,j,k) astarzz = ivp_astarzz(i,j,k) delxgxx = (gxx(i+1,j,k) - gxx(i-1,j,k))/(2.0D0*dx) delxgxy = (gxy(i+1,j,k) - gxy(i-1,j,k))/(2.0D0*dx) delxgxz = (gxz(i+1,j,k) - gxz(i-1,j,k))/(2.0D0*dx) delxgyy = (gyy(i+1,j,k) - gyy(i-1,j,k))/(2.0D0*dx) delxgyz = (gyz(i+1,j,k) - gyz(i-1,j,k))/(2.0D0*dx) delxgzz = (gzz(i+1,j,k) - gzz(i-1,j,k))/(2.0D0*dx) delygxx = (gxx(i,j+1,k) - gxx(i,j-1,k))/(2.0D0*dy) delygxy = (gxy(i,j+1,k) - gxy(i,j-1,k))/(2.0D0*dy) delygxz = (gxz(i,j+1,k) - gxz(i,j-1,k))/(2.0D0*dy) delygyy = (gyy(i,j+1,k) - gyy(i,j-1,k))/(2.0D0*dy) delygyz = (gyz(i,j+1,k) - gyz(i,j-1,k))/(2.0D0*dy) delygzz = (gzz(i,j+1,k) - gzz(i,j-1,k))/(2.0D0*dy) delzgxx = (gxx(i,j,k+1) - gxx(i,j,k-1))/(2.0D0*dz) delzgxy = (gxy(i,j,k+1) - gxy(i,j,k-1))/(2.0D0*dz) delzgxz = (gxz(i,j,k+1) - gxz(i,j,k-1))/(2.0D0*dz) delzgyy = (gyy(i,j,k+1) - gyy(i,j,k-1))/(2.0D0*dz) delzgyz = (gyz(i,j,k+1) - gyz(i,j,k-1))/(2.0D0*dz) delzgzz = (gzz(i,j,k+1) - gzz(i,j,k-1))/(2.0D0*dz) uxx = ivp_uxx(i,j,k) uxy = ivp_uxy(i,j,k) uxz = ivp_uxz(i,j,k) uyy = ivp_uyy(i,j,k) uyz = ivp_uyz(i,j,k) uzz = ivp_uzz(i,j,k) ivp_wx_source(i,j,k) = ivp_wx_source(i,j,k) + & 5.d-1*astarxx*delxgxx*uxx-5.d-1*astaryy*delxgyy*uxx-astaryz*delx & gyz*uxx-5.d-1*astarzz*delxgzz*uxx+astarxy*delygxx*uxx+astaryy*de & lygxy*uxx+astaryz*delygxz*uxx+astarxz*delzgxx*uxx+astaryz*delzgx & y*uxx+astarzz*delzgxz*uxx+astarxx*delxgxy*uxy+astarxy*delxgyy*ux & y+astarxz*delxgyz*uxy-5.d-1*astarxx*delygxx*uxy-astarxz*delygxz* & uxy+5.d-1*astaryy*delygyy*uxy-5.d-1*astarzz*delygzz*uxy+astarxz* & delzgxy*uxy+astaryz*delzgyy*uxy+astarzz*delzgyz*uxy+astarxx*delx & gxz*uxz+astarxy*delxgyz*uxz+astarxz*delxgzz*uxz+astarxy*delygxz* & uxz+astaryy*delygyz*uxz+astaryz*delygzz*uxz-5.d-1*astarxx*delzgx & x*uxz-astarxy*delzgxy*uxz-5.d-1*astaryy*delzgyy*uxz+5.d-1*astarz & z*delzgzz*uxz ivp_wy_source(i,j,k) = ivp_wy_source(i,j,k) + & 5.d-1*astarxx*delxgxx*uxy-5.d-1*astaryy*delxgyy*uxy-astaryz*delx & gyz*uxy-5.d-1*astarzz*delxgzz*uxy+astarxy*delygxx*uxy+astaryy*de & lygxy*uxy+astaryz*delygxz*uxy+astarxz*delzgxx*uxy+astaryz*delzgx & y*uxy+astarzz*delzgxz*uxy+astarxx*delxgxy*uyy+astarxy*delxgyy*uy & y+astarxz*delxgyz*uyy-5.d-1*astarxx*delygxx*uyy-astarxz*delygxz* & uyy+5.d-1*astaryy*delygyy*uyy-5.d-1*astarzz*delygzz*uyy+astarxz* & delzgxy*uyy+astaryz*delzgyy*uyy+astarzz*delzgyz*uyy+astarxx*delx & gxz*uyz+astarxy*delxgyz*uyz+astarxz*delxgzz*uyz+astarxy*delygxz* & uyz+astaryy*delygyz*uyz+astaryz*delygzz*uyz-5.d-1*astarxx*delzgx & x*uyz-astarxy*delzgxy*uyz-5.d-1*astaryy*delzgyy*uyz+5.d-1*astarz & z*delzgzz*uyz ivp_wz_source(i,j,k) = ivp_wz_source(i,j,k) + & 5.d-1*astarxx*delxgxx*uxz-5.d-1*astaryy*delxgyy*uxz-astaryz*delx & gyz*uxz-5.d-1*astarzz*delxgzz*uxz+astarxy*delygxx*uxz+astaryy*de & lygxy*uxz+astaryz*delygxz*uxz+astarxz*delzgxx*uxz+astaryz*delzgx & y*uxz+astarzz*delzgxz*uxz+astarxx*delxgxy*uyz+astarxy*delxgyy*uy & z+astarxz*delxgyz*uyz-5.d-1*astarxx*delygxx*uyz-astarxz*delygxz* & uyz+5.d-1*astaryy*delygyy*uyz-5.d-1*astarzz*delygzz*uyz+astarxz* & delzgxy*uyz+astaryz*delzgyy*uyz+astarzz*delzgyz*uyz+astarxx*delx & gxz*uzz+astarxy*delxgyz*uzz+astarxz*delxgzz*uzz+astarxy*delygxz* & uzz+astaryy*delygyz*uzz+astaryz*delygzz*uzz-5.d-1*astarxx*delzgx & x*uzz-astarxy*delzgxy*uzz-5.d-1*astaryy*delzgyy*uzz+5.d-1*astarz & z*delzgzz*uzz enddo enddo enddo return end