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

[rhaas at tapir.caltech.edu]

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

Modified:
 /trunk/src/
  GRHydro_Bondi_new.F90

Log:
 GRHydro_InitData: exchange find_bondi_solution_bracket for find_bondi_solution
 and improve smoothing of v field.
 
 From: Bruno Coutinho Mundim <bcmsma at astro.rit.edu>

File Changes:

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

File [modified]: GRHydro_Bondi_new.F90
Delta lines: +10 -7
===================================================================
--- trunk/src/GRHydro_Bondi_new.F90	2012-11-09 01:59:02 UTC (rev 177)
+++ trunk/src/GRHydro_Bondi_new.F90	2012-11-09 01:59:06 UTC (rev 178)
@@ -162,13 +162,14 @@
 !!$                                               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 ((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))
+!!$  call find_bondi_solution_bracket(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot, &
+!!$       rho_bondi(j),rho_bondi(j+1))
+  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
@@ -176,8 +177,9 @@
   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))
+!!$  call find_bondi_solution_bracket( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot, &
+!!$       rho_bondi(j),rho_bondi(j+1))
+  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 
@@ -224,8 +226,9 @@
 
               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_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