Add find_bondi_solution_bracket and velocity smoothing from Josh's patch.

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

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

1 files changed, 188 insertions, 12 deletions

diff --git a/src/GRHydro_Bondi_new.F90 b/src/GRHydro_Bondi_new.F90
index 6327080..e161043 100644
--- a/src/GRHydro_Bondi_new.F90
+++ b/src/GRHydro_Bondi_new.F90
@@ -58,7 +58,7 @@ subroutine GRHydro_Bondi_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
 
@@ -112,7 +112,8 @@ subroutine GRHydro_Bondi_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) 
@@ -126,7 +127,8 @@ subroutine GRHydro_Bondi_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)
@@ -162,23 +164,39 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS)
 
 !!$    // find the derivative near r=M
   rnew = 2.25 * M
-  j = floor ( 0.5 +  (log10(rnew) - logrmin) / dlogr )  
+  j = floor ((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 )
+  call find_bondi_solution_bracket(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot, &
+       rho_bondi(j),rho_bondi(j+1))
   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
   rhocheck2 = rho_bondi(j)
-  call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot )
+  call find_bondi_solution_bracket( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot, &
+       rho_bondi(j),rho_bondi(j+1))
+  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
@@ -196,12 +214,14 @@ subroutine GRHydro_Bondi_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
 
-              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)) 
               rho(i,j,k) = rhotmp
               uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso)
            else
@@ -211,7 +231,14 @@ subroutine GRHydro_Bondi_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
            
@@ -382,3 +409,152 @@ subroutine find_bondi_solution(r, rho, u, v, rs, rhos, M, Mdot, Kval, gam, Qdot
  
 end subroutine find_bondi_solution
 
+    subroutine find_bondi_solution_bracket(r, rho, u, v, rs, rhos, M, Mdot, Kval, gam, Qdot,rho1,rho2 )
+      implicit none
+    
+      CCTK_REAL rho1,rho2
+      CCTK_REAL :: r, rho, u, v, rs, rhos, M, Mdot, Kval, gam, Qdot
+      CCTK_REAL :: ur,r_sol, rho_old
+      CCTK_REAL :: f, df, dx, x_old, resid, jac
+      CCTK_REAL :: errx, x_orig, x
+      CCTK_INT :: n_iter, i_extra, doing_extra, keep_iterating, i_increase
+      CCTK_REAL :: vp, h, hp, term
+      
+      CCTK_REAL :: newt_tol_b,small_bondi
+      CCTK_INT :: max_newt_iter_b, extra_newt_iter_b
+
+      if(rho.gt.rho1.or.rho.lt.rho2) then
+         write(6,*)'Very Bad rho! ',rho1,rho2,rho
+         stop
+      endif
+    
+      max_newt_iter_b = 30
+      newt_tol_b = 1.0e-15
+      extra_newt_iter_b = 2
+      small_bondi = 1.0e-20
+      
+!!$  if(r>8.1043 .and. r<8.1044)write(*,*)'init guess:',r,rho
+!!$  write(*,*)'init guess:',r,rho
+      
+!!$      if (rho < 0.) then  
+!!$         if( r > 0.9*rs .and. r < 1.1*rs ) then 
+!!$            rho = rhos
+!!$         else
+!!$            if(r < rs) then
+!!$               ur = r**(-0.5)
+!!$            else     
+!!$               ur = 0.5*r**(-1.5)
+!!$            endif
+!!$            rho = Mdot / (4.*M_PI * r * r * ur)
+!!$         endif
+!!$      endif
+      
+!!$  if(r>8.1043 .and. r<8.1044)
+!!$  write(*,*)'init guess:',r,ur,rho,rs,rhos,Mdot
+!!$  if(r<1.0001e-10)write(*,*)'init guess:',r,ur,rho,rs,rhos,Mdot
+!!$  write(*,*)'init guess:',r,ur,rho,rs,rhos,Mdot
+      
+!!$ set global variables needed by residual function:
+
+      r_sol = r 
+      
+      
+!!$ Use Newton's method to find rho:
+!!$  gnr_bondi( rho, NEWT_DIM_B, bondi_resid)     
+      errx = 1.0
+      df=1.0
+      f=1.0
+      doing_extra = 0
+      
+      rho_old=rho
+      x=rho
+      
+      n_iter = 0
+      
+!!$ Start the Newton-Raphson iterations : 
+      keep_iterating = 1
+      do while( keep_iterating == 1 ) 
+         
+         hp  =  Kval * gam *  x**(gam - 2.)   !!$   //   dh/drho 
+         h   =  1. +  hp * x / ( gam - 1. )
+         v   =  Mdot / ( 4. * M_PI * r_sol * r_sol * x )    
+         vp  =  -v / x   !!$ //  dv/drho
+         term = 1. - 2.*M/r_sol + v*v
+         resid  =  -Qdot  +  h * h * term
+         jac =  2. * h *( hp*term + h*v*vp )
+         dx = -resid / jac
+         f = 0.5*resid**2
+         df = -2. * f
+         
+         errx = 0.
+         x_old = x
+         
+         x=x+dx
+
+         if(x.gt.rho1.or.x.lt.rho2) then
+!!$            write(6,*)'Bad rho! ',rho1,rho2,rho
+            if(x.gt.rho1)x=0.5*(x_old+rho1)
+            if(x.lt.rho2)x=0.5*(x_old+rho2)
+         endif
+
+         if(x.gt.rho1.or.x.lt.rho2) then
+            write(6,*)'Bad rho, bad! ',rho1,rho2,rho
+            stop
+         endif
+         
+!!$     if(r>8.1043 .and. r<8.1044)write(*,*)'iter:',x,dx,resid,jac
+         
+!!$    /****************************************/
+!!$    /* Calculate the convergence criterion */
+!!$    /****************************************/
+         
+!!$    /* For the new criterion, always look at relative error in indep. variable: */
+!!$    // METHOD specific:
+         if(x==0) then
+            errx = abs(dx)
+            x = small_bondi
+         else
+            errx = abs(dx/x)
+         endif
+         
+!!$    /*****************************************************************************/
+!!$    /* If we've reached the tolerance level, then just do a few extra iterations */
+!!$    /*  before stopping                                                          */
+!!$    /*****************************************************************************/
+         
+!!$     if(r>8.1043 .and. r<8.1044)write(*,*)'iter3:',errx,newt_tol_b,keep_iterating, &
+!!$          doing_extra,i_extra
+         
+         if((abs(errx)<=newt_tol_b) .and. (doing_extra == 0) .and. (extra_newt_iter_b > 0)) &
+              doing_extra=1
+         
+         if( doing_extra == 1 ) i_extra=i_extra+1 
+         
+         if( ((abs(errx) <= newt_tol_b).and.(doing_extra == 0)) .or. &
+              (i_extra > extra_newt_iter_b) .or. (n_iter >= (max_newt_iter_b-1)) ) then
+            keep_iterating = 0
+            if(n_iter>=(max_newt_iter_b-1))write(6,*)'Extra iterations!'
+         endif
+         
+!!$     if(r>8.1043 .and. r<8.1044)write(*,*)'iter4:',errx,newt_tol_b,keep_iterating, &
+!!$          doing_extra,i_extra
+         
+         
+         n_iter=n_iter+1
+         
+      end do
+      
+      rho=x
+      
+!!$ Calculate other quantities:
+      u = Kval * rho**( gam ) / (gam - 1.)
+      v = Mdot / ( 4. * M_PI * r * r * rho )
+      
+!!$ if(r>8.1043 .and. r<8.1044)
+!!$ write(*,*)'final:',r,rho,u,v
+      
+      
+      return
+      
+    end subroutine find_bondi_solution_bracket
+