path: root/src/GRHydro_SourceAM.F90

 /*@@
   @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, 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 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,&
  !$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)

        sqrtdet = sdetg(i,j,k)
        call UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, sqrtdet*sqrtdet, 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