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@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