diff options
Diffstat (limited to 'src/GRHydro_SourceM.F90')
| -rw-r--r-- | src/GRHydro_SourceM.F90 | 30 |
1 files changed, 20 insertions, 10 deletions
diff --git a/src/GRHydro_SourceM.F90 b/src/GRHydro_SourceM.F90 index cf88c62..bcd2b9e 100644 --- a/src/GRHydro_SourceM.F90 +++ b/src/GRHydro_SourceM.F90 @@ -395,7 +395,10 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) halfTdgz = half*(txx*dz_gxx + tyy*dz_gyy + tzz*dz_gzz) +& txy*dz_gxy + txz*dz_gxz + tyz*dz_gyz - sx_source = t00*& + + + + 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) +& @@ -403,7 +406,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) rhohstarW2*invalp*(vlowx*dx_betax + vlowy*dx_betay + vlowz*dx_betaz) -& bt*(bxlow*dx_betax + bylow*dx_betay + bzlow*dx_betaz) - sy_source = t00*& + 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) +& @@ -411,8 +414,8 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) 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 + & + 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) +& @@ -423,20 +426,21 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) 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 + taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source if(clean_divergence.ne.0) then + psidcrhs(i,j,k) = -1.d0 * (kap_dc*alp(i,j,k) + & dx_betax + dy_betay + dz_betaz ) * psidc(i,j,k) + & - Bconsx(i,j,k) * ( dx_alp - half*alp(i,j,k) * & + Bconsx(i,j,k) * (dx_alp - half*alp(i,j,k) * & ( uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + 2.d0*uxy*dx_gxy + & - 2.d0*uxz*dx_gxz + 2.d0*uyz*dx_gyz ) )/ sqrtdet + & + 2.d0*uxz*dx_gxz + 2.d0*uyz*dx_gyz ) )/sqrtdet + & Bconsy(i,j,k) * (dy_alp - half*alp(i,j,k) * & ( uxx*dy_gxx + uyy*dy_gyy + uzz*dy_gzz + 2.d0*uxy*dy_gxy + & - 2.d0*uxz*dy_gxz + 2.d0*uyz*dy_gyz ) )/ sqrtdet + & + 2.d0*uxz*dy_gxz + 2.d0*uyz*dy_gyz ) )/sqrtdet + & Bconsz(i,j,k) * (dz_alp - half*alp(i,j,k) * & - ( uxx*dz_gxx + uyy*dy_gyy + uzz*dy_gzz + 2.d0*uxy*dy_gxy + & - 2.d0*uxz*dz_gxz + 2.d0*uyz*dy_gyz ) )/ sqrtdet + ( uxx*dz_gxx + uyy*dz_gyy + uzz*dz_gzz + 2.d0*uxy*dz_gxy + & + 2.d0*uxz*dz_gxz + 2.d0*uyz*dz_gyz ) )/sqrtdet !!$ g^{jk} d_i g_{kj} = d_i (g) / det dx_det_bydet = uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + & @@ -511,6 +515,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do @@ -525,6 +530,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do @@ -539,6 +545,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do @@ -553,6 +560,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do @@ -567,6 +575,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do @@ -581,6 +590,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do |