/*@@ @file GRHydro_Source.F90 @date Sat Jan 26 02:03:56 2002 @author 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) /*@@ @routine SourceTerms @date Sat Jan 26 02:04:21 2002 @author 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 SourceTerms(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, rhoenthalpyW2 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 :: psi4pt, dx_psi4, dy_psi4, dz_psi4 CCTK_REAL :: shiftshiftk, shiftkx, shiftky, shiftkz CCTK_REAL :: sumTK CCTK_REAL :: halfshiftdgx, halfshiftdgy, halfshiftdgz CCTK_REAL :: halfTdgx, halfTdgy, halfTdgz CCTK_REAL :: invalp, invalp2 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 if (evolve_tracer .ne. 0) then cons_tracerrhs = 0.0d0 end if if (evolve_Y_e .ne. 0) then y_e_con_rhs = 0.0d0 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 psi4pt,det,sqrtdet,rhoenthalpyW2,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,t00,t0x,t0y,t0z,txx,txy,txz,tyy,tyz,tzz,& !$OMP dx_alp,dy_alp,dz_alp,tau_source,sx_source,sy_source,sz_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 dx_psi4,dy_psi4,dz_psi4,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) !!$ All the matter variables (except velocity) always appear !!$ together in this form rhoenthalpyW2 = (rho(i,j,k)*(one + eps(i,j,k)) + press(i,j,k))*& w_lorentz(i,j,k)**2 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 vlowx = velx(i,j,k)*localgxx + vely(i,j,k)*localgxy +& velz(i,j,k)*localgxz vlowy = velx(i,j,k)*localgxy + vely(i,j,k)*localgyy +& velz(i,j,k)*localgyz vlowz = velx(i,j,k)*localgxz + vely(i,j,k)*localgyz +& velz(i,j,k)*localgzz !!$ For a change, these are T^{ij} t00 = (rhoenthalpyW2 - press(i,j,k))*invalp2 t0x = rhoenthalpyW2*velxshift/alp(i,j,k) +& press(i,j,k)*shiftx*invalp2 t0y = rhoenthalpyW2*velyshift/alp(i,j,k) +& press(i,j,k)*shifty*invalp2 t0z = rhoenthalpyW2*velzshift/alp(i,j,k) +& press(i,j,k)*shiftz*invalp2 txx = rhoenthalpyW2*velxshift*velxshift +& press(i,j,k)*(uxx - shiftx*shiftx*invalp2) txy = rhoenthalpyW2*velxshift*velyshift +& press(i,j,k)*(uxy - shiftx*shifty*invalp2) txz = rhoenthalpyW2*velxshift*velzshift +& press(i,j,k)*(uxz - shiftx*shiftz*invalp2) tyy = rhoenthalpyW2*velyshift*velyshift +& press(i,j,k)*(uyy - shifty*shifty*invalp2) tyz = rhoenthalpyW2*velyshift*velzshift +& press(i,j,k)*(uyz - shifty*shiftz*invalp2) tzz = rhoenthalpyW2*velzshift*velzshift +& press(i,j,k)*(uzz - shiftz*shiftz*invalp2) !!$ 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) else dx_alp = DIFF_X_4(alp) dy_alp = DIFF_Y_4(alp) dz_alp = DIFF_Z_4(alp) 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) +& 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) +& halfTdgx + rhoenthalpyW2*& (vlowx*dx_betax + vlowy*dx_betay + vlowz*dx_betaz)*& invalp sy_source = t00*& (halfshiftdgy - alp(i,j,k)*dy_alp) +& 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) +& halfTdgy + rhoenthalpyW2*& (vlowx*dy_betax + vlowy*dy_betay + vlowz*dy_betaz)*& invalp sz_source = t00*& (halfshiftdgz - alp(i,j,k)*dz_alp) +& 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) +& halfTdgz + rhoenthalpyW2*& (vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz)*& invalp 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 enddo enddo enddo !$OMP END PARALLEL DO deallocate(force_spatial_second_order) end subroutine SourceTerms