[Commits] [svn:einsteintoolkit] GRHydro/branches/divclean/src/ (Rev. 231)

[tanja.bode at physics.gatech.edu]

User: tbode
Date: 2011/04/27 03:30 PM

Modified:
 /branches/divclean/src/
  GRHydro_SourceM.F90

Log:
 Modify rhs for divergence cleaning.

File Changes:

Directory: /branches/divclean/src/
==================================

File [modified]: GRHydro_SourceM.F90
Delta lines: +50 -3
===================================================================
--- branches/divclean/src/GRHydro_SourceM.F90	2011-04-27 17:51:47 UTC (rev 230)
+++ branches/divclean/src/GRHydro_SourceM.F90	2011-04-27 20:30:04 UTC (rev 231)
@@ -33,6 +33,12 @@
 #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)
+#define Bconsx(i,j,k) Bcons(i,j,k,1)
+#define Bconsy(i,j,k) Bcons(i,j,k,2)
+#define Bconsz(i,j,k) Bcons(i,j,k,3)
+#define Bconsrhsx(i,j,k) Bconsrhs(i,j,k,1)
+#define Bconsrhsy(i,j,k) Bconsrhs(i,j,k,2)
+#define Bconsrhsz(i,j,k) Bconsrhs(i,j,k,3)
 
  /*@@
    @routine    SourceTermsM
@@ -76,8 +82,9 @@
   CCTK_REAL :: halfshiftdgx, halfshiftdgy, halfshiftdgz
   CCTK_REAL :: halfTdgx, halfTdgy, halfTdgz
   CCTK_REAL :: invalp, invalp2
-  CCTK_REAL :: dx_psidc, dy_psidc, dz_psidc
+  !!CCTK_REAL :: dx_psidc, dy_psidc, dz_psidc
   CCTK_REAL :: bvcx_source, bvcy_source, bvcz_source
+  CCTK_REAL :: dx_det_bydet, dy_det_bydet, dz_det_bydet
 
   CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow,bdotv,b2,dum,bxlow,bylow,bzlow
   CCTK_REAL :: bt,bx,by,bz,rhohstarW2,pstar
@@ -103,7 +110,7 @@
   densrhs = 0.d0
   srhs = 0.d0
   taurhs = 0.d0
-  Bvecrhs=0.d0
+  Bconsrhs=0.d0
 
   if (evolve_tracer .ne. 0) then
     cons_tracerrhs = 0.d0
@@ -436,7 +443,47 @@
         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
-           psidcrhs(i,j,k) = -1.d0*ch_dc*ch_dc/cp_dc/cp_dc*psidc(i,j,k)
+           psidcrhs(i,j,k) = -1.d0 * ( ch_dc*ch_dc/cp_dc/cp_dc*alp(i,j,k) + &
+                dx_shiftx + dy_shifty + dz_shiftz ) * psidc(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 + &
+                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 + &
+                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
+
+           !!$ g^{jk} d_i g_{kj} = d_i (g) / det
+           dx_det_bydet = uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + &
+                2.d0*(uxy*dx_gxy+uxz*dx_gxz+uyz*dx_gyz)
+           dy_det_bydet = uxx*dy_gxx + uyy*dy_gyy + uzz*dy_gzz + &
+                2.d0*(uxy*dy_gxy+uxz*dy_gxz+uyz*dy_gyz)
+           dz_det_bydet = uxx*dz_gxx + uyy*dz_gyy + uzz*dz_gzz + &
+                2.d0*(uxy*dz_gxy+uxz*dz_gxz+uyz*dz_gyz)
+
+           bvcx_source = -1.d0 * ( Bconsx(i,j,k)*dx_shiftx + &
+                Bconsy(i,j,k)*dy_shiftx + Bconsz(i,j,k)*dz_shiftx ) + &
+                psidc(i,j,k)*sqrtdet*( uxx*dx_alp+uxy*dy_alp+uxz*dz_alp ) + &
+                psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxx*dx_det_bydet + &
+                  uxy*dy_det_bydet + uxz*dz_det_bydet )
+
+           bvcy_source = -1.d0 * ( Bconsx(i,j,k)*dx_shifty + &
+                Bconsy(i,j,k)*dy_shifty + Bconsz(i,j,k)*dz_shifty ) + &
+                psidc(i,j,k)*sqrtdet*( uxy*dx_alp+uyy*dy_alp+uyz*dz_alp ) + &
+                psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxy*dx_det_bydet + &
+                  uyy*dy_det_bydet + uyz*dz_det_bydet )
+
+           bvcz_source = -1.d0 * ( Bconsx(i,j,k)*dx_shiftz + &
+                Bconsy(i,j,k)*dy_shiftz + Bconsz(i,j,k)*dz_shiftz ) + &
+                psidc(i,j,k)*sqrtdet*( uxz*dx_alp+uyz*dy_alp+uzz*dz_alp ) + &
+                psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxz*dx_det_bydet + &
+                  uyz*dy_det_bydet + uzz*dz_det_bydet )
+
+           Bconsrhsx(i,j,k) = bvcx_source
+           Bconsrhsy(i,j,k) = bvcy_source
+           Bconsrhsz(i,j,k) = bvcz_source
         endif
 
       enddo