aboutsummaryrefslogtreecommitdiff
path: root/src/GRHydro_SourceAM.F90
diff options
context:
space:
mode:
Diffstat (limited to 'src/GRHydro_SourceAM.F90')
-rw-r--r--src/GRHydro_SourceAM.F90679
1 files changed, 679 insertions, 0 deletions
diff --git a/src/GRHydro_SourceAM.F90 b/src/GRHydro_SourceAM.F90
new file mode 100644
index 0000000..1e0b190
--- /dev/null
+++ b/src/GRHydro_SourceAM.F90
@@ -0,0 +1,679 @@
+ /*@@
+ @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) vup(i,j,k,1)
+#define vely(i,j,k) vup(i,j,k,2)
+#define velz(i,j,k) vup(i,j,k,3)
+#define Bvecx(i,j,k) Bprim(i,j,k,1)
+#define Bvecy(i,j,k) Bprim(i,j,k,2)
+#define Bvecz(i,j,k) Bprim(i,j,k,3)
+#define Avecx(i,j,k) Avec(i,j,k,1)
+#define Avecy(i,j,k) Avec(i,j,k,2)
+#define Avecz(i,j,k) Avec(i,j,k,3)
+#define Avecrhsx(i,j,k) Avecrhs(i,j,k,1)
+#define Avecrhsy(i,j,k) Avecrhs(i,j,k,2)
+#define Avecrhsz(i,j,k) Avecrhs(i,j,k,3)
+
+ /*@@
+ @routine SourceTermsAM
+ @date Aug 30, 2010
+ @author Tanja Bode, 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 SourceTermsAM(CCTK_ARGUMENTS)
+
+ implicit none
+
+ ! save memory when MP is not used
+ ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1
+ TARGET gaa, gab, gac, gbb, gbc, gcc
+ TARGET gxx, gxy, gxz, gyy, gyz, gzz
+ TARGET kaa, kab, kac, kbb, kbc, kcc
+ TARGET kxx, kxy, kxz, kyy, kyz, kzz
+ TARGET betaa, betab, betac
+ TARGET betax, betay, betaz
+ TARGET lvel, vel
+ TARGET lBvec, Bvec
+
+ 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 :: Avcx_source, Avcy_source, Avcz_source
+ CCTK_REAL :: dx_det_bydet, dy_det_bydet, dz_det_bydet
+ CCTK_REAL :: gdg_x, gdg_y, gdg_z !! g^{ik} d_k g_{ij}
+
+ CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow,bdotv,b2,dum1,dum2,bxlow,bylow,bzlow
+ CCTK_REAL :: bt,bx,by,bz,rhohstarW2,pstar
+
+ logical, allocatable, dimension (:,:,:) :: force_spatial_second_order
+
+ ! save memory when MP is not used
+ CCTK_INT :: GRHydro_UseGeneralCoordinates
+ CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33
+ CCTK_REAL, DIMENSION(:,:,:), POINTER :: k11, k12, k13, k22, k23, k33
+ CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3
+ CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim
+
+ if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
+ g11 => gaa
+ g12 => gab
+ g13 => gac
+ g22 => gbb
+ g23 => gbc
+ g33 => gcc
+ k11 => kaa
+ k12 => kab
+ k13 => kac
+ k22 => kbb
+ k23 => kbc
+ k33 => kcc
+ beta1 => betaa
+ beta2 => betab
+ beta3 => betac
+ vup => lvel
+ Bprim => lBvec
+ else
+ g11 => gxx
+ g12 => gxy
+ g13 => gxz
+ g22 => gyy
+ g23 => gyz
+ g33 => gzz
+ k11 => kxx
+ k12 => kxy
+ k13 => kxz
+ k22 => kyy
+ k23 => kyz
+ k33 => kzz
+ beta1 => betax
+ beta2 => betay
+ beta3 => betaz
+ vup => vel
+ Bprim => Bvec
+ end if
+#define gxx faulty_gxx
+#define gxy faulty_gxy
+#define gxz faulty_gxz
+#define gyy faulty_gyy
+#define gyz faulty_gyz
+#define gzz faulty_gzz
+#define kxx faulty_kxx
+#define kxy faulty_kxy
+#define kxz faulty_kxz
+#define kyy faulty_kyy
+#define kyz faulty_kyz
+#define kzz faulty_kzz
+#define betax faulty_betax
+#define betay faulty_betay
+#define betaz faulty_betaz
+#define gaa faulty_gaa
+#define gab faulty_gab
+#define gac faulty_gac
+#define gbb faulty_gbb
+#define gbc faulty_gbc
+#define gcc faulty_gcc
+#define kaa faulty_kaa
+#define kab faulty_kab
+#define kac faulty_kac
+#define kbb faulty_kbb
+#define kbc faulty_kbc
+#define kcc faulty_kcc
+#define betaa faulty_betaa
+#define betab faulty_betab
+#define betac faulty_betac
+#define vel faulty_vel
+#define Bvec faulty_Bvec
+
+ 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
+ Avecrhs = 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
+
+ if (track_divB .ne. 0) then
+ divB=0.d0
+ endif
+
+ if (transport_constraints .ne. 0) then
+ Evec = 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, k)
+ 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, k, 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,dum1,dum2,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,&
+ !$OMP tau_source,sx_source,sy_source,sz_source,&
+ !$OMP dx_det_bydet,dy_det_bydet,dz_det_bydet,&
+ !$OMP gdg_x,gdg_y,gdg_z,&
+ !$OMP Avcx_source,Avcy_source,Avcz_source,evolve_Lorenz_gge,&
+ !$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 = g11(i,j,k)
+ localgxy = g12(i,j,k)
+ localgxz = g13(i,j,k)
+ localgyy = g22(i,j,k)
+ localgyz = g23(i,j,k)
+ localgzz = g33(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)
+
+
+ shiftx = beta1(i,j,k)
+ shifty = beta2(i,j,k)
+ shiftz = beta3(i,j,k)
+
+ if (local_spatial_order .eq. 2) then
+
+ dx_betax = DIFF_X_2(beta1)
+ dx_betay = DIFF_X_2(beta2)
+ dx_betaz = DIFF_X_2(beta3)
+
+ dy_betax = DIFF_Y_2(beta1)
+ dy_betay = DIFF_Y_2(beta2)
+ dy_betaz = DIFF_Y_2(beta3)
+
+ dz_betax = DIFF_Z_2(beta1)
+ dz_betay = DIFF_Z_2(beta2)
+ dz_betaz = DIFF_Z_2(beta3)
+
+ else
+
+ dx_betax = DIFF_X_4(beta1)
+ dx_betay = DIFF_X_4(beta2)
+ dx_betaz = DIFF_X_4(beta3)
+
+ dy_betax = DIFF_Y_4(beta1)
+ dy_betay = DIFF_Y_4(beta2)
+ dy_betaz = DIFF_Y_4(beta3)
+
+ dz_betax = DIFF_Z_4(beta1)
+ dz_betay = DIFF_Z_4(beta2)
+ dz_betaz = DIFF_Z_4(beta3)
+
+ end if
+
+ 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,dum1,dum2,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)+0.5d0*b2
+
+!!$ 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)
+
+ 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(g11)
+ dx_gxy = DIFF_X_2(g12)
+ dx_gxz = DIFF_X_2(g13)
+ dx_gyy = DIFF_X_2(g22)
+ dx_gyz = DIFF_X_2(g23)
+ dx_gzz = DIFF_X_2(g33)
+ dy_gxx = DIFF_Y_2(g11)
+ dy_gxy = DIFF_Y_2(g12)
+ dy_gxz = DIFF_Y_2(g13)
+ dy_gyy = DIFF_Y_2(g22)
+ dy_gyz = DIFF_Y_2(g23)
+ dy_gzz = DIFF_Y_2(g33)
+ dz_gxx = DIFF_Z_2(g11)
+ dz_gxy = DIFF_Z_2(g12)
+ dz_gxz = DIFF_Z_2(g13)
+ dz_gyy = DIFF_Z_2(g22)
+ dz_gyz = DIFF_Z_2(g23)
+ dz_gzz = DIFF_Z_2(g33)
+
+ else
+
+ dx_gxx = DIFF_X_4(g11)
+ dx_gxy = DIFF_X_4(g12)
+ dx_gxz = DIFF_X_4(g13)
+ dx_gyy = DIFF_X_4(g22)
+ dx_gyz = DIFF_X_4(g23)
+ dx_gzz = DIFF_X_4(g33)
+ dy_gxx = DIFF_Y_4(g11)
+ dy_gxy = DIFF_Y_4(g12)
+ dy_gxz = DIFF_Y_4(g13)
+ dy_gyy = DIFF_Y_4(g22)
+ dy_gyz = DIFF_Y_4(g23)
+ dy_gzz = DIFF_Y_4(g33)
+ dz_gxx = DIFF_Z_4(g11)
+ dz_gxy = DIFF_Z_4(g12)
+ dz_gxz = DIFF_Z_4(g13)
+ dz_gyy = DIFF_Z_4(g22)
+ dz_gyz = DIFF_Z_4(g23)
+ dz_gzz = DIFF_Z_4(g33)
+
+ end if
+
+!!$ Contract the shift with the extrinsic curvature
+
+ shiftshiftk = shiftx*shiftx*k11(i,j,k) + &
+ shifty*shifty*k22(i,j,k) + &
+ shiftz*shiftz*k33(i,j,k) + &
+ two*(shiftx*shifty*k12(i,j,k) + &
+ shiftx*shiftz*k13(i,j,k) + &
+ shifty*shiftz*k23(i,j,k))
+
+ shiftkx = shiftx*k11(i,j,k) + shifty*k12(i,j,k) + shiftz*k13(i,j,k)
+ shiftky = shiftx*k12(i,j,k) + shifty*k22(i,j,k) + shiftz*k23(i,j,k)
+ shiftkz = shiftx*k13(i,j,k) + shifty*k23(i,j,k) + shiftz*k33(i,j,k)
+
+!!$ Contract the matter terms with the extrinsic curvature
+
+ sumTK = txx*k11(i,j,k) + tyy*k22(i,j,k) + tzz*k33(i,j,k) &
+ + two*(txy*k12(i,j,k) + txz*k13(i,j,k) + tyz*k23(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) + halfTdgz + &
+ 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)
+
+ !! B^i and A^i both live in cell centers currently
+ Avcx_source = sqrtdet*(vely(i,j,k)*Bvecz(i,j,k) - velz(i,j,k)*Bvecy(i,j,k))
+ Avcy_source = sqrtdet*(velz(i,j,k)*Bvecx(i,j,k) - velx(i,j,k)*Bvecz(i,j,k))
+ Avcz_source = sqrtdet*(velx(i,j,k)*Bvecy(i,j,k) - vely(i,j,k)*Bvecx(i,j,k))
+
+ if ( evolve_Lorenz_gge.gt.0 ) then
+ Aphi(i,j,k) = 0.d0
+ end if
+
+ 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
+ Avecrhsx(i,j,k) = Avcx_source
+ Avecrhsy(i,j,k) = Avcy_source
+ Avecrhsz(i,j,k) = Avcz_source
+
+ enddo
+ enddo
+ enddo
+ !$OMP END PARALLEL DO
+
+ deallocate(force_spatial_second_order)
+
+#if(0) /* poison edges of domain */
+ if(last_iteration_seen .ne. cctk_iteration .or. reflevel .ne. grhydro_reflevel) then
+ last_iteration_seen = cctk_iteration
+ reflevel = grhydro_reflevel
+ mol_substep = 0
+ else
+ mol_substep = mol_substep + 1
+ end if
+ do k = 1, GRHydro_stencil*mol_substep
+ do j = 1, cctk_lsh(2)
+ do i = 1, cctk_lsh(1)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Avecx(i,j,k) = -1d100
+ Avecy(i,j,k) = -1d100
+ Avecz(i,j,k) = -1d100
+ if ( evolve_Lorenz_gge.gt.0 ) then
+ Aphi(i,j,k) = -1d100
+ end if
+ end do
+ end do
+ end do
+ do k = cctk_lsh(3)-GRHydro_stencil*mol_substep+1, cctk_lsh(3)
+ do j = 1, cctk_lsh(2)
+ do i = 1, cctk_lsh(1)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Avecx(i,j,k) = -1d100
+ Avecy(i,j,k) = -1d100
+ Avecz(i,j,k) = -1d100
+ if ( evolve_Lorenz_gge.gt.0 ) then
+ Aphi(i,j,k) = -1d100
+ end if
+ end do
+ end do
+ end do
+ do i = 1, GRHydro_stencil*mol_substep
+ do k = 1, cctk_lsh(3)
+ do j = 1, cctk_lsh(2)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Avecx(i,j,k) = -1d100
+ Avecy(i,j,k) = -1d100
+ Avecz(i,j,k) = -1d100
+ if ( evolve_Lorenz_gge.gt.0 ) then
+ Aphi(i,j,k) = -1d100
+ end if
+ end do
+ end do
+ end do
+ do i = cctk_lsh(1)-GRHydro_stencil*mol_substep+1, cctk_lsh(1)
+ do k = 1, cctk_lsh(3)
+ do j = 1, cctk_lsh(2)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Avecx(i,j,k) = -1d100
+ Avecy(i,j,k) = -1d100
+ Avecz(i,j,k) = -1d100
+ if ( evolve_Lorenz_gge.gt.0 ) then
+ Aphi(i,j,k) = -1d100
+ end if
+ end do
+ end do
+ end do
+ do j = 1, GRHydro_stencil*mol_substep
+ do i = 1, cctk_lsh(1)
+ do k = 1, cctk_lsh(3)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Avecx(i,j,k) = -1d100
+ Avecy(i,j,k) = -1d100
+ Avecz(i,j,k) = -1d100
+ if ( evolve_Lorenz_gge.gt.0 ) then
+ Aphi(i,j,k) = -1d100
+ end if
+ end do
+ end do
+ end do
+ do j = cctk_lsh(2)-GRHydro_stencil*mol_substep+1, cctk_lsh(2)
+ do i = 1, cctk_lsh(1)
+ do k = 1, cctk_lsh(3)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Avecx(i,j,k) = -1d100
+ Avecy(i,j,k) = -1d100
+ Avecz(i,j,k) = -1d100
+ if ( evolve_Lorenz_gge.gt.0 ) then
+ Aphi(i,j,k) = -1d100
+ end if
+ end do
+ end do
+ end do
+#endif
+
+#undef faulty_gxx
+#undef faulty_gxy
+#undef faulty_gxz
+#undef faulty_gyy
+#undef faulty_gyz
+#undef faulty_gzz
+#undef faulty_vel
+#undef faulty_Bvec
+#undef faulty_gxx_p
+#undef faulty_gxy_p
+#undef faulty_gxz_p
+#undef faulty_gyy_p
+#undef faulty_gyz_p
+#undef faulty_gzz_p
+#undef faulty_vel_p
+#undef faulty_Bvec_p
+#undef faulty_gxx_p_p
+#undef faulty_gxy_p_p
+#undef faulty_gxz_p_p
+#undef faulty_gyy_p_p
+#undef faulty_gyz_p_p
+#undef faulty_gzz_p_p
+#undef faulty_vel_p_p
+#undef faulty_Bvec_p_p
+end subroutine SourceTermsAM
+
+
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-26 15:04:21 +0000