/*@@ @file GRHydro_SourceM.F90 @date Oct 11, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke @desc The geometric source terms for the matter evolution @enddesc @@*/ ! Second order f.d. #define DIFF_X_2(q) (0.5d0 * (q(i+1,j,k) - q(i-1,j,k)) * idx) #define DIFF_Y_2(q) (0.5d0 * (q(i,j+1,k) - q(i,j-1,k)) * idy) #define DIFF_Z_2(q) (0.5d0 * (q(i,j,k+1) - q(i,j,k-1)) * idz) ! Fourth order f.d. #define DIFF_X_4(q) ((-q(i+2,j,k) + 8.d0 * q(i+1,j,k) - 8.d0 * q(i-1,j,k) + \ q(i-2,j,k)) / 12.d0 * idx) #define DIFF_Y_4(q) ((-q(i,j+2,k) + 8.d0 * q(i,j+1,k) - 8.d0 * q(i,j-1,k) + \ q(i,j-2,k)) / 12.d0 * idy) #define DIFF_Z_4(q) ((-q(i,j,k+2) + 8.d0 * q(i,j,k+1) - 8.d0 * q(i,j,k-1) + \ q(i,j,k-2)) / 12.d0 * idz) #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "GRHydro_Macros.h" #define velx(i,j,k) vel(i,j,k,1) #define vely(i,j,k) vel(i,j,k,2) #define velz(i,j,k) vel(i,j,k,3) #define Bvecx(i,j,k) Bvec(i,j,k,1) #define Bvecy(i,j,k) Bvec(i,j,k,2) #define Bvecz(i,j,k) Bvec(i,j,k,3) /*@@ @routine SourceTermsM @date Aug 30, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke @desc Calculate the geometric source terms and add to the update GFs @enddesc @calls @calledby @history Minor alterations of routine from GR3D. @endhistory @@*/ subroutine SourceTermsM(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k, nx, ny, nz CCTK_REAL :: one, two, half CCTK_REAL :: t00, t0x, t0y, t0z, txx, txy, txz, tyy, tyz, tzz CCTK_REAL :: sqrtdet, det, uxx, uxy, uxz, uyy, uyz, uzz CCTK_REAL :: shiftx, shifty, shiftz, velxshift, velyshift, velzshift CCTK_REAL :: vlowx, vlowy, vlowz CCTK_REAL :: dx_betax, dx_betay, dx_betaz, dy_betax, dy_betay,& dy_betaz, dz_betax, dz_betay, dz_betaz CCTK_REAL :: dx_alp, dy_alp, dz_alp CCTK_REAL :: tau_source, sx_source, sy_source, sz_source CCTK_REAL :: localgxx,localgxy,localgxz,localgyy,localgyz,localgzz CCTK_REAL :: dx_gxx, dx_gxy, dx_gxz, dx_gyy, dx_gyz, dx_gzz CCTK_REAL :: dy_gxx, dy_gxy, dy_gxz, dy_gyy, dy_gyz, dy_gzz CCTK_REAL :: dz_gxx, dz_gxy, dz_gxz, dz_gyy, dz_gyz, dz_gzz CCTK_REAL :: dx, dy, dz, idx, idy, idz CCTK_REAL :: shiftshiftk, shiftkx, shiftky, shiftkz CCTK_REAL :: sumTK CCTK_REAL :: halfshiftdgx, halfshiftdgy, halfshiftdgz CCTK_REAL :: halfTdgx, halfTdgy, halfTdgz CCTK_REAL :: invalp, invalp2 CCTK_REAL :: dx_psidc, dy_psidc, dz_psidc CCTK_REAL :: bvcx_source, bvcy_source, bvcz_source CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow,bdotv,b2,dum,bxlow,bylow,bzlow CCTK_REAL :: bt,bx,by,bz,rhohstarW2,pstar logical, allocatable, dimension (:,:,:) :: force_spatial_second_order one = 1.0d0 two = 2.0d0 half = 0.5d0 nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) dx = CCTK_DELTA_SPACE(1) dy = CCTK_DELTA_SPACE(2) dz = CCTK_DELTA_SPACE(3) idx = 1.d0/dx idy = 1.d0/dy idz = 1.d0/dz !!$ Initialize the update terms to be zero. !!$ This will guarantee that no garbage in the boundaries is updated. densrhs = 0.d0 srhs = 0.d0 taurhs = 0.d0 Bvecrhs=0.d0 if (evolve_tracer .ne. 0) then cons_tracerrhs = 0.d0 end if if (evolve_Y_e .ne. 0) then Y_e_con_rhs = 0.0d0 endif if (clean_divergence .ne. 0) then psidcrhs=0.d0 endif !!$ Set up the array for checking the order. We switch to second order !!$ differencing at boundaries and near excision regions. !!$ Copied straight from BSSN. allocate (force_spatial_second_order(nx,ny,nz)) force_spatial_second_order = .FALSE. if (spatial_order > 2) then !$OMP PARALLEL DO PRIVATE(i, j) do k = 1 + GRHydro_stencil, nz - GRHydro_stencil do j = 1 + GRHydro_stencil, ny - GRHydro_stencil do i = 1 + GRHydro_stencil, nx - GRHydro_stencil 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 force_spatial_second_order(i,j,k) = .TRUE. else if ( use_mask > 0 ) then if (minval(emask(i-2:i+2,j-2:j+2,k-2:k+2)) < 0.75d0) then force_spatial_second_order(i,j,k) = .TRUE. end if end if end do end do end do !$OMP END PARALLEL DO end if !$OMP PARALLEL DO PRIVATE(i, j, local_spatial_order,& !$OMP localgxx,localgxy,localgxz,localgyy,localgyz,localgzz,& !$OMP det,sqrtdet,shiftx,shifty,shiftz,& !$OMP dx_betax,dx_betay,dx_betaz,dy_betax,dy_betay,dy_betaz,& !$OMP dz_betax,dz_betay,dz_betaz,velxshift,velyshift,velzshift,& !$OMP vlowx,vlowy,vlowz,Bvecxlow,Bvecylow,Bveczlow, & !$OMP bdotv,b2,dum,bxlow,bylow,bzlow,bt,bx,by,bz,rhohstarW2,pstar,& !$OMP t00,t0x,t0y,t0z,txx,txy,txz,tyy,tyz,tzz,& !$OMP dx_alp,dy_alp,dz_alp,dx_psidc,dy_psidc,dz_psidc,& !$OMP tau_source,sx_source,sy_source,sz_source,& !$OMP bvcx_source,bvcy_source,bvcz_source,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz,& !$OMP dx_gxx, dx_gxy, dx_gxz, dx_gyy, dx_gyz, dx_gzz,& !$OMP dy_gxx, dy_gxy, dy_gxz, dy_gyy, dy_gyz, dy_gzz,& !$OMP dz_gxx, dz_gxy, dz_gxz, dz_gyy, dz_gyz, dz_gzz,& !$OMP shiftshiftk,shiftkx,shiftky,shiftkz,& !$OMP sumTK,halfshiftdgx,halfshiftdgy,halfshiftdgz,& !$OMP halfTdgx,halfTdgy,halfTdgz,invalp,invalp2) do k=1 + GRHydro_stencil,nz - GRHydro_stencil do j=1 + GRHydro_stencil,ny - GRHydro_stencil do i=1 + GRHydro_stencil,nx - GRHydro_stencil local_spatial_order = spatial_order if (force_spatial_second_order(i,j,k)) then local_spatial_order = 2 end if !!$ Set the metric terms. localgxx = gxx(i,j,k) localgxy = gxy(i,j,k) localgxz = gxz(i,j,k) localgyy = gyy(i,j,k) localgyz = gyz(i,j,k) localgzz = gzz(i,j,k) det = SPATIAL_DETERMINANT(localgxx, localgxy, localgxz,\ localgyy, localgyz, localgzz) sqrtdet = sqrt(det) call UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, det, localgxx,& localgxy, localgxz, localgyy, localgyz, localgzz) if (shift_state .ne. 0) then shiftx = betax(i,j,k) shifty = betay(i,j,k) shiftz = betaz(i,j,k) if (local_spatial_order .eq. 2) then dx_betax = DIFF_X_2(betax) dx_betay = DIFF_X_2(betay) dx_betaz = DIFF_X_2(betaz) dy_betax = DIFF_Y_2(betax) dy_betay = DIFF_Y_2(betay) dy_betaz = DIFF_Y_2(betaz) dz_betax = DIFF_Z_2(betax) dz_betay = DIFF_Z_2(betay) dz_betaz = DIFF_Z_2(betaz) else dx_betax = DIFF_X_4(betax) dx_betay = DIFF_X_4(betay) dx_betaz = DIFF_X_4(betaz) dy_betax = DIFF_Y_4(betax) dy_betay = DIFF_Y_4(betay) dy_betaz = DIFF_Y_4(betaz) dz_betax = DIFF_Z_4(betax) dz_betay = DIFF_Z_4(betay) dz_betaz = DIFF_Z_4(betaz) end if else shiftx = 0.0d0 shifty = 0.0d0 shiftz = 0.0d0 dx_betax = 0.0d0 dx_betay = 0.0d0 dx_betaz = 0.0d0 dy_betax = 0.0d0 dy_betay = 0.0d0 dy_betaz = 0.0d0 dz_betax = 0.0d0 dz_betay = 0.0d0 dz_betaz = 0.0d0 endif invalp = 1.0d0 / alp(i,j,k) invalp2 = invalp**2 velxshift = velx(i,j,k) - shiftx*invalp velyshift = vely(i,j,k) - shifty*invalp velzshift = velz(i,j,k) - shiftz*invalp call calc_vlow_blow(localgxx,localgxy,localgxz,localgyy,localgyz,localgzz, & velx(i,j,k),vely(i,j,k),velz(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k), & vlowx,vlowy,vlowz,Bvecxlow,Bvecylow,Bveczlow, & bdotv,b2,dum,dum,bxlow,bylow,bzlow) !!$ These are the contravariant components bt = w_lorentz(i,j,k)/alp(i,j,k)*bdotv bx = Bvecx(i,j,k)/w_lorentz(i,j,k)+w_lorentz(i,j,k)*bdotv*velxshift by = Bvecy(i,j,k)/w_lorentz(i,j,k)+w_lorentz(i,j,k)*bdotv*velyshift bz = Bvecz(i,j,k)/w_lorentz(i,j,k)+w_lorentz(i,j,k)*bdotv*velzshift rhohstarW2 = (rho(i,j,k)*(one + eps(i,j,k)) + press(i,j,k)+ b2)*& w_lorentz(i,j,k)**2 pstar = press(i,j,k)+b2/2.d0 !!$ For a change, these are T^{ij} t00 = (rhohstarW2 - pstar)*invalp2-bt**2 t0x = rhohstarW2*velxshift*invalp +& pstar*shiftx*invalp2-bt*bx t0y = rhohstarW2*velyshift*invalp +& pstar*shifty*invalp2-bt*by t0z = rhohstarW2*velzshift*invalp +& pstar*shiftz*invalp2-bt*bz txx = rhohstarW2*velxshift*velxshift +& pstar*(uxx - shiftx*shiftx*invalp2)-bx**2 txy = rhohstarW2*velxshift*velyshift +& pstar*(uxy - shiftx*shifty*invalp2)-bx*by txz = rhohstarW2*velxshift*velzshift +& pstar*(uxz - shiftx*shiftz*invalp2)-bx*bz tyy = rhohstarW2*velyshift*velyshift +& pstar*(uyy - shifty*shifty*invalp2)-by**2 tyz = rhohstarW2*velyshift*velzshift +& pstar*(uyz - shifty*shiftz*invalp2)-by*bz tzz = rhohstarW2*velzshift*velzshift +& pstar*(uzz - shiftz*shiftz*invalp2)-bz**2 !!$ Derivatives of the lapse, metric and shift if (local_spatial_order .eq. 2) then dx_alp = DIFF_X_2(alp) dy_alp = DIFF_Y_2(alp) dz_alp = DIFF_Z_2(alp) if(clean_divergence.ne.0) then dx_psidc = DIFF_X_2(psidc) dy_psidc = DIFF_Y_2(psidc) dz_psidc = DIFF_Z_2(psidc) endif else dx_alp = DIFF_X_4(alp) dy_alp = DIFF_Y_4(alp) dz_alp = DIFF_Z_4(alp) if(clean_divergence.ne.0) then dx_psidc = DIFF_X_4(psidc) dy_psidc = DIFF_Y_4(psidc) dz_psidc = DIFF_Z_4(psidc) endif end if if (local_spatial_order .eq. 2) then dx_gxx = DIFF_X_2(gxx) dx_gxy = DIFF_X_2(gxy) dx_gxz = DIFF_X_2(gxz) dx_gyy = DIFF_X_2(gyy) dx_gyz = DIFF_X_2(gyz) dx_gzz = DIFF_X_2(gzz) dy_gxx = DIFF_Y_2(gxx) dy_gxy = DIFF_Y_2(gxy) dy_gxz = DIFF_Y_2(gxz) dy_gyy = DIFF_Y_2(gyy) dy_gyz = DIFF_Y_2(gyz) dy_gzz = DIFF_Y_2(gzz) dz_gxx = DIFF_Z_2(gxx) dz_gxy = DIFF_Z_2(gxy) dz_gxz = DIFF_Z_2(gxz) dz_gyy = DIFF_Z_2(gyy) dz_gyz = DIFF_Z_2(gyz) dz_gzz = DIFF_Z_2(gzz) else dx_gxx = DIFF_X_4(gxx) dx_gxy = DIFF_X_4(gxy) dx_gxz = DIFF_X_4(gxz) dx_gyy = DIFF_X_4(gyy) dx_gyz = DIFF_X_4(gyz) dx_gzz = DIFF_X_4(gzz) dy_gxx = DIFF_Y_4(gxx) dy_gxy = DIFF_Y_4(gxy) dy_gxz = DIFF_Y_4(gxz) dy_gyy = DIFF_Y_4(gyy) dy_gyz = DIFF_Y_4(gyz) dy_gzz = DIFF_Y_4(gzz) dz_gxx = DIFF_Z_4(gxx) dz_gxy = DIFF_Z_4(gxy) dz_gxz = DIFF_Z_4(gxz) dz_gyy = DIFF_Z_4(gyy) dz_gyz = DIFF_Z_4(gyz) dz_gzz = DIFF_Z_4(gzz) end if !!$ Contract the shift with the extrinsic curvature shiftshiftk = shiftx*shiftx*kxx(i,j,k) + & shifty*shifty*kyy(i,j,k) + & shiftz*shiftz*kzz(i,j,k) + & two*(shiftx*shifty*kxy(i,j,k) + & shiftx*shiftz*kxz(i,j,k) + & shifty*shiftz*kyz(i,j,k)) shiftkx = shiftx*kxx(i,j,k) + shifty*kxy(i,j,k) + shiftz*kxz(i,j,k) shiftky = shiftx*kxy(i,j,k) + shifty*kyy(i,j,k) + shiftz*kyz(i,j,k) shiftkz = shiftx*kxz(i,j,k) + shifty*kyz(i,j,k) + shiftz*kzz(i,j,k) !!$ Contract the matter terms with the extrinsic curvature sumTK = txx*kxx(i,j,k) + tyy*kyy(i,j,k) + tzz*kzz(i,j,k) & + two*(txy*kxy(i,j,k) + txz*kxz(i,j,k) + tyz*kyz(i,j,k)) !!$ Update term for tau tau_source = t00* & (shiftshiftk - (shiftx*dx_alp + shifty*dy_alp + shiftz*dz_alp) )& + t0x*(-dx_alp + two*shiftkx) & + t0y*(-dy_alp + two*shiftky) & + t0z*(-dz_alp + two*shiftkz) & + sumTK !!$ The following looks very little like the terms in the !!$ standard papers. Take a look in the ThornGuide to see why !!$ it is really the same thing. !!$ Contract the shift with derivatives of the metric halfshiftdgx = half*(shiftx*shiftx*dx_gxx + & shifty*shifty*dx_gyy + shiftz*shiftz*dx_gzz) + & shiftx*shifty*dx_gxy + shiftx*shiftz*dx_gxz + & shifty*shiftz*dx_gyz halfshiftdgy = half*(shiftx*shiftx*dy_gxx + & shifty*shifty*dy_gyy + shiftz*shiftz*dy_gzz) + & shiftx*shifty*dy_gxy + shiftx*shiftz*dy_gxz + & shifty*shiftz*dy_gyz halfshiftdgz = half*(shiftx*shiftx*dz_gxx + & shifty*shifty*dz_gyy + shiftz*shiftz*dz_gzz) + & shiftx*shifty*dz_gxy + shiftx*shiftz*dz_gxz + & shifty*shiftz*dz_gyz !!$ Contract the matter with derivatives of the metric halfTdgx = half*(txx*dx_gxx + tyy*dx_gyy + tzz*dx_gzz) +& txy*dx_gxy + txz*dx_gxz + tyz*dx_gyz halfTdgy = half*(txx*dy_gxx + tyy*dy_gyy + tzz*dy_gzz) +& txy*dy_gxy + txz*dy_gxz + tyz*dy_gyz halfTdgz = half*(txx*dz_gxx + tyy*dz_gyy + tzz*dz_gzz) +& txy*dz_gxy + txz*dz_gxz + tyz*dz_gyz sx_source = t00*& (halfshiftdgx - alp(i,j,k)*dx_alp) + halfTdgx + & t0x*(shiftx*dx_gxx + shifty*dx_gxy + shiftz*dx_gxz) +& t0y*(shiftx*dx_gxy + shifty*dx_gyy + shiftz*dx_gyz) +& t0z*(shiftx*dx_gxz + shifty*dx_gyz + shiftz*dx_gzz) +& rhohstarW2*invalp*(vlowx*dx_betax + vlowy*dx_betay + vlowz*dx_betaz) -& bt*(bxlow*dx_betax + bylow*dx_betay + bzlow*dx_betaz) sy_source = t00*& (halfshiftdgy - alp(i,j,k)*dy_alp) + halfTdgy + & t0x*(shiftx*dy_gxx + shifty*dy_gxy + shiftz*dy_gxz) +& t0y*(shiftx*dy_gxy + shifty*dy_gyy + shiftz*dy_gyz) +& t0z*(shiftx*dy_gxz + shifty*dy_gyz + shiftz*dy_gzz) +& rhohstarW2*invalp*(vlowx*dy_betax + vlowy*dy_betay + vlowz*dy_betaz) -& bt*(bxlow*dy_betax + bylow*dy_betay + bzlow*dy_betaz) sz_source = t00*& (halfshiftdgz - alp(i,j,k)*dz_alp) + halfTdgy + & t0x*(shiftx*dz_gxx + shifty*dz_gxy + shiftz*dz_gxz) +& t0y*(shiftx*dz_gxy + shifty*dz_gyy + shiftz*dz_gyz) +& t0z*(shiftx*dz_gxz + shifty*dz_gyz + shiftz*dz_gzz) +& rhohstarW2*invalp*(vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz) -& bt*(bxlow*dz_betax + bylow*dz_betay + bzlow*dz_betaz) if(clean_divergence.ne.0) then bvcx_source = uxx*dx_psidc+uxy*dy_psidc+uxz*dz_psidc bvcy_source = uxy*dx_psidc+uyy*dy_psidc+uyz*dz_psidc bvcz_source = uxz*dx_psidc+uyz*dy_psidc+uzz*dz_psidc endif densrhs(i,j,k) = 0.d0 srhs(i,j,k,1) = alp(i,j,k)*sqrtdet*sx_source srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sy_source srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sz_source taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source if(clean_divergence.ne.0) then Bvecrhs(i,j,k,1) = alp(i,j,k)*sqrtdet*bvcx_source Bvecrhs(i,j,k,2) = alp(i,j,k)*sqrtdet*bvcy_source Bvecrhs(i,j,k,3) = alp(i,j,k)*sqrtdet*bvcz_source psidcrhs(i,j,k) = -1.d0*ch_dc*ch_dc/cp_dc/cp_dc*psidc(i,j,k) endif enddo enddo enddo !$OMP END PARALLEL DO deallocate(force_spatial_second_order) end subroutine SourceTermsM