GRHydro_InitData: Bondi routine: smooth puncture data.

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

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

1 files changed, 56 insertions, 16 deletions

diff --git a/src/GRHydro_Bondi_new.F90 b/src/GRHydro_Bondi_new.F90
index 9e56ef2..9b9d45e 100644
--- a/src/GRHydro_Bondi_new.F90
+++ b/src/GRHydro_Bondi_new.F90
@@ -51,12 +51,14 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
   DECLARE_CCTK_FUNCTIONS
   
   CCTK_INT :: i, j, k, nx, ny, nz, imin, jb,N_points
-  PARAMETER (N_points=2000)
+  CCTK_REAL :: ONEmTINY
+  PARAMETER (N_points=2000,ONEmTINY=0.999999d0)
   CCTK_REAL :: M, Msq, Mdot, rs, gam, rmin_bondi, rmax_bondi, cs_sq,cs,vs_sq,vs,rhos,gtemp,hs, Kval, Qdot
   CCTK_REAL :: logrmin,dlogr,rhotmp,utmp,vtmp,rspher
   CCTK_REAL :: r_bondi(N_points), logr_bondi(N_points), rho_bondi(N_points), u_bondi(N_points), v_bondi(N_points)
   CCTK_REAL :: drhodr, det, rhocheck, rhocheck2, riso, rnew, rsch, ucheck, psinew
   CCTK_REAL :: uiso, uisocheck, vcheck, ucheck2, vcheck2, xhat,yhat, zhat, xp, yp, zp
+  CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm
 
   !!$set_bondi_parameters
   M = bondi_central_mass(1)
@@ -65,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)
@@ -89,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
@@ -132,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
@@ -150,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)
@@ -168,12 +170,45 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
            xhat = xp/riso
            yhat = yp/riso
            zhat = zp/riso
+	   roverm = riso/M
            
-           if(use_smooth_puncture_data.ne.0 .and. riso.lt.M) then
-              psinew=2.875 - 5*(riso/M)**2 + 6*(riso/M)**4
-              gxx(i,j,k) = psinew**4
-              gyy(i,j,k) = psinew**4
-              gzz(i,j,k) = psinew**4
+           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
@@ -187,9 +222,9 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
               gyz(i,j,k) = 0.0
            endif
 
-           alp(i,j,k) = 1.0/gxx(i,j,k)**2
+           alp(i,j,k) = 1.0/sqrt(gxx(i,j,k))
 
-           if(riso > M) then
+           if(roverm > ONEmTINY) then
               rsch = 0.25 * ( 2.*riso + M)**2 / riso
               jb = floor( 0.5 +  (log10(rsch) - logrmin) / dlogr ) 
 
@@ -200,7 +235,7 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
               rho(i,j,k) = rhotmp
               uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso)
            else
-              if(riso > 0.5*M) then
+              if(roverm > 0.5d0*ONEmTINY) then
                  rho(i,j,k) = rhocheck+drhodr*riso*(riso-M)/M
               else 
                  rho(i,j,k) = (rhocheck-drhodr*M/4.0)*(1.-cos(2.*M_PI*riso/M))/2.0
@@ -224,6 +259,11 @@ 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)	
+
         end do
      end do
   end do