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