GRHydro_InitData: Move the metric fudge from Bondi routine to SmoothedPuncture thorn

From: Bruno Coutinho Mundim <bcmsma@astro.rit.edu>

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@163 ac85fae7-cede-4708-beff-ae01c7fa1c26

2 files changed, 11 insertions, 77 deletions

diff --git a/param.ccl b/param.ccl
index b3c050f..466fbe0 100644..100755
--- a/param.ccl
+++ b/param.ccl
@@ -99,20 +99,6 @@ REAL shock_radius "Radius of sperical shock"
   0.0:*   :: "Anything positive"
 } 1.0
 
-BOOLEAN bondi_use_smooth_puncture_data "Smooth puncture solution inside R=M using 4th order matching polynomial"
-{
-} "no"
-
-
-BOOLEAN bondi_smooth_gxx "Smooth puncture data by matching to gxx not psi"
-{
-} "yes"
-
-REAL bondi_smooth_puncture_radius "Radius inside which to mathc functions, in units of R/M"
-{
- 0:0.5 :: "Less than/equal to 0.5"
-} 0.5
-
 BOOLEAN change_shock_direction "Change the shock direction"
 { 
 } "no"
diff --git a/src/GRHydro_Bondi_new.F90 b/src/GRHydro_Bondi_new.F90
index 9b9d45e..b539c5f 100644..100755
--- a/src/GRHydro_Bondi_new.F90
+++ b/src/GRHydro_Bondi_new.F90
@@ -67,7 +67,7 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
   rs = r_sonicpt_bondi
   gam = gl_gamma
 
-  !!write(*,*) 'Bondi_pars:',M,mdot_sonicpt_bondi,r_sonicpt_bondi,gl_gamma
+  write(*,*) 'Bondi_pars:',M,mdot_sonicpt_bondi,r_sonicpt_bondi,gl_gamma
 
   rmin_bondi = M * bondi_rmin(1)
   rmax_bondi = M * bondi_rmax(1)
@@ -91,7 +91,7 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
   logrmin = log10(rmin_bondi)
   dlogr = (log10(rmax_bondi) - logrmin)/(1.*(N_points-1))
 
-  !!write(*,*)'More pars:',cs,vs,rhos,hs,Kval,Qdot,logrmin,dlogr
+  write(*,*)'More pars:',cs,vs,rhos,hs,Kval,Qdot,logrmin,dlogr
 
   rhotmp=1.0d30
   imin=1
@@ -134,10 +134,10 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
      v_bondi(i)   = vtmp
   enddo
   
-  !!write(*,*)"i=1:",r_bondi(1),rho_bondi(1),u_bondi(1),v_bondi(1)
-  !!write(*,*)"i=100:",r_bondi(100),rho_bondi(100),u_bondi(100),v_bondi(100)
-  !!write(*,*)"i=1000:",r_bondi(1000),rho_bondi(1000),u_bondi(1000),v_bondi(1000)
-  !!write(*,*)"i=1500:",r_bondi(1500),rho_bondi(1500),u_bondi(1500),v_bondi(1500)
+  write(*,*)"i=1:",r_bondi(1),rho_bondi(1),u_bondi(1),v_bondi(1)
+  write(*,*)"i=100:",r_bondi(100),rho_bondi(100),u_bondi(100),v_bondi(100)
+  write(*,*)"i=1000:",r_bondi(1000),rho_bondi(1000),u_bondi(1000),v_bondi(1000)
+  write(*,*)"i=1500:",r_bondi(1500),rho_bondi(1500),u_bondi(1500),v_bondi(1500)
 
 !!$    // find the derivative near r=M
   rnew = 2.25 * M
@@ -152,7 +152,7 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
   call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot )
   drhodr = 100.0*(rhocheck2-rhocheck)/M 
 
-  !!write(6,*)'Rhocheck:',rhocheck,rhocheck2,drhodr
+  write(6,*)'Rhocheck:',rhocheck,rhocheck2,drhodr
 
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
@@ -172,58 +172,6 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
            zhat = zp/riso
 	   roverm = riso/M
            
-           if(bondi_use_smooth_puncture_data.ne.0 .and. roverm.lt.(bondi_smooth_puncture_radius*ONEmTINY)) then
-
-	      rsm=bondi_smooth_puncture_radius
-		!!write(6,*)'smoothing at a radius',rsm
-
-	      if(bondi_smooth_gxx.ne.0) then
-	
-  		f=(1.d0+0.5d0/rsm)**4
-  		df=-1.d0*(2.d0*rsm+1.d0)**3/4.d0/rsm**5
-  		ddf=(2.d0*rsm+1.d0)**2*(4.d0*rsm+5.d0)/4.d0/rsm**6
-
-  		a=f-0.625d0*rsm*df+0.125d0*rsm**2*ddf
-		b=0.75d0/rsm*df-0.25d0*ddf
-		c=-0.125d0/rsm**3*df+0.125d0/rsm**2*ddf
-     		
-		!!write(6,*)'smooth gxx:',rsm,f,df,ddf,a,b,c
-
-  	        gxx(i,j,k) = a + b*roverm**2 + c*roverm**4
-	        gyy(i,j,k) = a + b*roverm**2 + c*roverm**4
-	        gzz(i,j,k) = a + b*roverm**2 + c*roverm**4
-
-	      else
-		f=1.d0+0.5d0/rsm
-		df=-0.5d0/rsm**2
-		ddf=1.d0/rsm**3
-
-		a=f-0.625d0*rsm*df+0.125d0*rsm**2*ddf
-		b=0.75d0/rsm*df-0.25d0*ddf
-		c=-0.125d0/rsm**3*df+0.125d0/rsm**2*ddf
-
-		!!write(6,*)'smooth psi:',rsm,f,df,ddf,a,b,c
-
-                psinew=a + b*roverm**2 + c*roverm**4
-                gxx(i,j,k) = psinew**4
-                gyy(i,j,k) = psinew**4
-                gzz(i,j,k) = psinew**4
-	      endif
-              gxy(i,j,k) = 0.0
-              gxz(i,j,k) = 0.0
-              gyz(i,j,k) = 0.0
-
-           else if(riso < 1.0e-7) then
-              gxx(i,j,k) = 1.0e4
-              gyy(i,j,k) = 1.0e4
-              gzz(i,j,k) = 1.0e4
-              gxy(i,j,k) = 0.0
-              gxz(i,j,k) = 0.0
-              gyz(i,j,k) = 0.0
-           endif
-
-           alp(i,j,k) = 1.0/sqrt(gxx(i,j,k))
-
            if(roverm > ONEmTINY) then
               rsch = 0.25 * ( 2.*riso + M)**2 / riso
               jb = floor( 0.5 +  (log10(rsch) - logrmin) / dlogr ) 
@@ -259,10 +207,10 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
 		  velx(i,j,k),vely(i,j,k),velz(i,j,k), &
 		  eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
            
-	  !! if(riso.gt.1.014d0.and.riso.lt.1.015)write(6,*)'Point to check:', &
-	  !!        x(i,j,k),y(i,j,k),z(i,j,k),riso,gxx(i,j,k),dens(i,j,k),tau(i,j,k),&
-	  !!        sx(i,j,k),sy(i,j,k),sz(i,j,k),rho(i,j,k),eps(i,j,k),&
-	  !!        velx(i,j,k),vely(i,j,k),velz(i,j,k)	
+	   if(riso.gt.1.014d0.and.riso.lt.1.015)write(6,*)'Point to check:', &
+	          x(i,j,k),y(i,j,k),z(i,j,k),riso,gxx(i,j,k),dens(i,j,k),tau(i,j,k),&
+	          sx(i,j,k),sy(i,j,k),sz(i,j,k),rho(i,j,k),eps(i,j,k),&
+	          velx(i,j,k),vely(i,j,k),velz(i,j,k)	
 
         end do
      end do