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

[rhaas at tapir.caltech.edu]

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

Modified:
 /trunk/src/
  GRHydro_BondiM_new.F90, GRHydro_Bondi_new.F90

Log:
 GRHydro_InitData: port smoother velocity matching to mag. bondi routine.
 
 From: Bruno Coutinho Mundim <bcmsma at astro.rit.edu>

File Changes:

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

File [modified]: GRHydro_BondiM_new.F90
Delta lines: +41 -12
===================================================================
--- trunk/src/GRHydro_BondiM_new.F90	2012-11-09 01:59:06 UTC (rev 178)
+++ trunk/src/GRHydro_BondiM_new.F90	2012-11-09 01:59:09 UTC (rev 179)
@@ -64,7 +64,7 @@
   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 @@
      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 @@
      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 @@
 !!$  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 @@
            
            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 @@
                  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
            

File [modified]: GRHydro_Bondi_new.F90
Delta lines: +1 -1
===================================================================
--- trunk/src/GRHydro_Bondi_new.F90	2012-11-09 01:59:06 UTC (rev 178)
+++ trunk/src/GRHydro_Bondi_new.F90	2012-11-09 01:59:09 UTC (rev 179)
@@ -148,7 +148,7 @@
                              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)