[Commits] [svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 177)

[rhaas at tapir.caltech.edu]

User: rhaas
Date: 2012/11/08 07:59 PM

Modified:
 /trunk/src/
  GRHydro_Bondi_new.F90

Log:
 GRHydro_InitData: improve error messages for find_bondi_solution_bracket func.
 
 From: Bruno Coutinho Mundim <bcmsma at astro.rit.edu>

File Changes:

Directory: /trunk/src/
======================

File [modified]: GRHydro_Bondi_new.F90
Delta lines: +14 -8
===================================================================
--- trunk/src/GRHydro_Bondi_new.F90	2012-11-09 01:59:00 UTC (rev 176)
+++ trunk/src/GRHydro_Bondi_new.F90	2012-11-09 01:59:02 UTC (rev 177)
@@ -164,7 +164,8 @@
 
 !!$    // find the derivative near r=M
   rnew = 2.25 * M
-  j = floor ((log10(rnew) - logrmin) / dlogr )+1.0
+  j = floor ((log10(rnew) - logrmin) / dlogr + 1.0)
+!!$  j = NINT((log10(rnew) - logrmin) / dlogr + 1.0)
   rhocheck = rho_bondi(j)
   call find_bondi_solution_bracket(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot, &
        rho_bondi(j),rho_bondi(j+1))
@@ -172,7 +173,8 @@
   
 !!$ 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((log10(rnew) - logrmin) / dlogr )+1.0
+  j = floor((log10(rnew) - logrmin) / dlogr + 1.0)
+!!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0)
   rhocheck2 = rho_bondi(j)
   call find_bondi_solution_bracket( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot, &
        rho_bondi(j),rho_bondi(j+1))
@@ -215,11 +217,13 @@
            if(roverm > ONEmTINY) then
               rsch = 0.25 * ( 2.*riso + M)**2 / riso
 !!$              jb = floor( 0.5 +  (log10(rsch) - logrmin) / dlogr ) 
-              jb = floor( (log10(rsch) - logrmin) / dlogr ) +1.0
+              jb = floor( (log10(rsch) - logrmin) / dlogr  + 1.0)
+!!$              jb = NINT((log10(rsch) - logrmin) / dlogr + 1.0)
 
-              if(jb >= N_points)jb = N_points-1
+              if(jb >= N_points) jb = N_points-1
 
-              rhotmp = rho_bondi(jb)+(rho_bondi(jb+1)-rho_bondi(jb))*(rsch-r_bondi(jb))/(r_bondi(jb+1)-r_bondi(jb))
+              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
@@ -424,9 +428,11 @@
       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
+         write(6,*)'find_bondi_solution_bracket: Very bad rho! (rholow,rhoup,rho) = ',rho1,rho2,rho
          stop
       endif
+
+         write(6,*)'find_bondi_solution_bracket: (rholow,rhoup,rho) = ',rho1,rho2,rho
     
       max_newt_iter_b = 30
       newt_tol_b = 1.0e-15
@@ -498,7 +504,7 @@
          endif
 
          if(x.gt.rho1.or.x.lt.rho2) then
-            write(6,*)'Bad rho, bad! ',rho1,rho2,rho
+            write(6,*)'find_bondi_solution_bracket: Bad rho, bad! (rholow,rhoup,rho) = ',rho1,rho2,rho
             stop
          endif
          
@@ -533,7 +539,7 @@
          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!'
+            if(n_iter>=(max_newt_iter_b-1))write(6,*)'find_bondi_solution_bracket: Extra iterations!'
          endif
          
 !!$     if(r>8.1043 .and. r<8.1044)write(*,*)'iter4:',errx,newt_tol_b,keep_iterating, &