GRHydro_InitData: port smoother velocity matching to mag. bondi routine.

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

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

2 files changed, 42 insertions, 13 deletions

diff --git a/src/GRHydro_BondiM_new.F90 b/src/GRHydro_BondiM_new.F90
index 30c1a8f..325e158 100644
--- a/src/GRHydro_BondiM_new.F90
+++ b/src/GRHydro_BondiM_new.F90
@@ -64,7 +64,7 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
   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
   CCTK_REAL :: uiso, uisocheck, vcheck, ucheck2, vcheck2, xhat,yhat, zhat, xp, yp, zp
-  CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm
+  CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm,dudr,uisocheck2,auiso,buiso
 
   character(400) :: debug_message
 
@@ -118,7 +118,8 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
      endif
   enddo
 
-  rhotmp = -1.  !!$ start with guess
+!!$  rhotmp = -1.  !!$ start with guess
+  rhotmp=rhos      !!$ start with value at sonic point!
 
   do i=imin,N_points
      rspher = r_bondi(i) 
@@ -132,7 +133,8 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
      v_bondi(i)   = vtmp
   end do
   
-  rhotmp = -1.
+!!$  rhotmp = -1.
+  rhotmp=rhos      !!$ start with value at sonic point!
   
   do i=imin-1,1,-1
      rspher = r_bondi(i)
@@ -165,27 +167,42 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
 !!$  write(debug_message,'(a,4f22.14)') "i=1500:",r_bondi(1500),rho_bondi(1500),&
 !!$                                               u_bondi(1500),v_bondi(1500)
 !!$  call CCTK_INFO(debug_message)
-  
 
-!!$    // find the derivative near r=M
+!!$    // find the derivative near r=M in isotropic coords = r=9/4M in schwarzschild; 
   rnew = 2.25 * M
-  j = floor ( 0.5 +  (log10(rnew) - logrmin) / dlogr )  
+  j = floor ((log10(rnew) - logrmin) / dlogr + 1.0)
+!!$  j = NINT((log10(rnew) - logrmin) / dlogr + 1.0)
   rhocheck = rho_bondi(j)
   call find_bondi_solution(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot )
   uisocheck = 4.0*vcheck/3.0
   
+!!$ the previous point was r=M in isotropic coords = r=9/4M in schwarzschild;  this one is r=1.01M in isotropic
   rnew = 0.25 * 3.02**2 * M/1.01
-  j = floor( 0.5 +  (log10(rnew) - logrmin) / dlogr ) 
+  j = floor((log10(rnew) - logrmin) / dlogr + 1.0)
+!!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0)
   rhocheck2 = rho_bondi(j)
   call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot )
+  uisocheck2 = vcheck2 / (1.0 - 1.0/2.02) / (1.0+ 1.0/2.02)
+
   drhodr = 100.0*(rhocheck2-rhocheck)/M 
 
+!!$ Don't divide by M here, to simplify the math
+  dudr   = 100.0*(uisocheck2-uisocheck)
+
   write(debug_message,'(a,3f22.14)') 'Rhocheck:',rhocheck,rhocheck2,drhodr
   call CCTK_INFO(debug_message)
+  write(debug_message,'(a,3f22.14)') 'Ucheck:',uisocheck,uisocheck2,dudr
+  call CCTK_INFO(debug_message)
 
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
   nz = cctk_lsh(3)
+
+
+  write(debug_message,'(a,3f22.14)') 'Lower bound coordinates:',x(1,1,1),y(1,1,1),z(1,1,1)
+  call CCTK_INFO(debug_message)
+  write(debug_message,'(a,3f22.14)') 'Upper bound coordinates:',x(nx,ny,nz),y(nx,ny,nz),z(nx,ny,nz)
+  call CCTK_INFO(debug_message)
   
   do i=1,nx
      do j=1,ny
@@ -203,12 +220,17 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
            
            if(roverm > ONEmTINY) then
               rsch = 0.25 * ( 2.*riso + M)**2 / riso
-              jb = floor( 0.5 +  (log10(rsch) - logrmin) / dlogr ) 
+!!$              jb = floor( 0.5 +  (log10(rsch) - logrmin) / dlogr ) 
+              jb = floor( (log10(rsch) - logrmin) / dlogr  + 1.0)
+!!$              jb = NINT((log10(rsch) - logrmin) / dlogr + 1.0)
 
-              if(jb > N_points)jb = N_points
+              if(jb >= N_points) jb = N_points-1
 
-              rhotmp = rho_bondi(jb)
-              call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot) 
+              rhotmp = rho_bondi(jb)+(rho_bondi(jb+1)-rho_bondi(jb))*&
+                  (rsch-r_bondi(jb))/(r_bondi(jb+1)-r_bondi(jb))
+!!$              call find_bondi_solution_bracket( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot, &
+!!$                   rho_bondi(jb),rho_bondi(jb+1)) 
+              call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot)
               rho(i,j,k) = rhotmp
               uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso)
            else
@@ -218,7 +240,14 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
                  rho(i,j,k) = (rhocheck-drhodr*M/4.0)*(1.-cos(2.*M_PI*riso/M))/2.0
               endif
               utmp = Kval * rho(i,j,k)**( gam ) / (gam - 1.)
-              uiso = uisocheck * riso / M
+
+!!$ match to uiso and dudr at roverm=1
+!!$ a R + b R^3 --->   a+b = uisocheck;  a+3b = dudr
+!!$ b = (dudr-uisocheck)/2;  a=3*uisocheck-dudr)/2
+
+              auiso = 1.5*uisocheck-0.5*dudr
+              buiso = 0.5*dudr-0.5*uisocheck
+              uiso =  roverm*(auiso+buiso*roverm**2)
            endif
            eps(i,j,k) = utmp/rhotmp
            
diff --git a/src/GRHydro_Bondi_new.F90 b/src/GRHydro_Bondi_new.F90
index a0077c2..b31ba7c 100644
--- a/src/GRHydro_Bondi_new.F90
+++ b/src/GRHydro_Bondi_new.F90
@@ -148,7 +148,7 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
                              u_bondi(i),v_bondi(i)
   end do
   close(47)
-  
+
 !!$  write(debug_message,'(a,4f22.14)') "i=1:",r_bondi(1),rho_bondi(1),&
 !!$                                            u_bondi(1),v_bondi(1)
 !!$  call CCTK_INFO(debug_message)