GRHydro_InitData: port bondi fixes to magnetized version.

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

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

1 files changed, 41 insertions, 26 deletions

diff --git a/src/GRHydro_BondiM_new.F90 b/src/GRHydro_BondiM_new.F90
index 3b1f47d..30c1a8f 100644
--- a/src/GRHydro_BondiM_new.F90
+++ b/src/GRHydro_BondiM_new.F90
@@ -57,12 +57,16 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
   DECLARE_CCTK_FUNCTIONS
   
   CCTK_INT :: i, j, k, nx, ny, nz, imin, jb,N_points
-  PARAMETER (N_points=2000)
+  CCTK_REAL :: ONEmTINY, tiny
+  PARAMETER (N_points=100000,ONEmTINY=0.999999d0,tiny=1.0d-12)
   CCTK_REAL :: M, Msq, Mdot, rs, gam, rmin_bondi, rmax_bondi, cs_sq,cs,vs_sq,vs,rhos,gtemp,hs, Kval, Qdot
   CCTK_REAL :: logrmin,dlogr,rhotmp,utmp,vtmp,rspher
   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
+
+  character(400) :: debug_message
 
   !!$set_bondi_parameters
   M = bondi_central_mass(1)
@@ -71,7 +75,9 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
   rs = r_sonicpt_bondi
   gam = gl_gamma
 
-  write(*,*) 'Bondi_pars:',M,mdot_sonicpt_bondi,r_sonicpt_bondi,gl_gamma
+  write(debug_message,'(a,4f22.14)') "Bondi_pars: ",M,mdot_sonicpt_bondi, &
+                                                    r_sonicpt_bondi,gl_gamma
+  call CCTK_INFO(debug_message)
 
   rmin_bondi = M * bondi_rmin(1)
   rmax_bondi = M * bondi_rmax(1)
@@ -95,7 +101,9 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
   logrmin = log10(rmin_bondi)
   dlogr = (log10(rmax_bondi) - logrmin)/(1.*(N_points-1))
 
-  write(*,*)'More pars:',cs,vs,rhos,hs,Kval,Qdot,logrmin,dlogr
+  write(debug_message,'(a,8f22.14)') "More pars: ",cs,vs,rhos,hs,Kval,Qdot,&
+                                                   logrmin,dlogr
+  call CCTK_INFO(debug_message)
 
   rhotmp=1.0d30
   imin=1
@@ -137,11 +145,27 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
      u_bondi(i)   = utmp
      v_bondi(i)   = vtmp
   enddo
+
+  open (47,file="bondi.asc",form="formatted")
+  do i=1,N_points   
+    write(47,'(i5,4f22.14)')i,r_bondi(i),rho_bondi(i),&
+                             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)
+!!$  write(debug_message,'(a,4f22.14)') "i=100:",r_bondi(100),rho_bondi(100),&
+!!$                                              u_bondi(100),v_bondi(100)
+!!$  call CCTK_INFO(debug_message)
+!!$  write(debug_message,'(a,4f22.14)') "i=1000:",r_bondi(1000),rho_bondi(1000),&
+!!$                                               u_bondi(1000),v_bondi(1000)
+!!$  call CCTK_INFO(debug_message)
+!!$  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)
   
-  write(*,*)"i=1:",r_bondi(1),rho_bondi(1),u_bondi(1),v_bondi(1)
-  write(*,*)"i=100:",r_bondi(100),rho_bondi(100),u_bondi(100),v_bondi(100)
-  write(*,*)"i=1000:",r_bondi(1000),rho_bondi(1000),u_bondi(1000),v_bondi(1000)
-  write(*,*)"i=1500:",r_bondi(1500),rho_bondi(1500),u_bondi(1500),v_bondi(1500)
 
 !!$    // find the derivative near r=M
   rnew = 2.25 * M
@@ -156,7 +180,8 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
   call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot )
   drhodr = 100.0*(rhocheck2-rhocheck)/M 
 
-  write(6,*)'Rhocheck:',rhocheck,rhocheck2,drhodr
+  write(debug_message,'(a,3f22.14)') 'Rhocheck:',rhocheck,rhocheck2,drhodr
+  call CCTK_INFO(debug_message)
 
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
@@ -174,19 +199,9 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
            xhat = xp/riso
            yhat = yp/riso
            zhat = zp/riso
+           roverm = riso/M
            
-           if(riso < 1.0e-7) then
-              gxx(i,j,k) = 1.0e4
-              gyy(i,j,k) = 1.0e4
-              gzz(i,j,k) = 1.0e4
-              gxy(i,j,k) = 0.0
-              gxz(i,j,k) = 0.0
-              gyz(i,j,k) = 0.0
-           endif
-
-           alp(i,j,k) = 1.0/gxx(i,j,k)**2
-
-           if(riso > M) then
+           if(roverm > ONEmTINY) then
               rsch = 0.25 * ( 2.*riso + M)**2 / riso
               jb = floor( 0.5 +  (log10(rsch) - logrmin) / dlogr ) 
 
@@ -197,7 +212,7 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
               rho(i,j,k) = rhotmp
               uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso)
            else
-              if(riso > 0.5*M) then
+              if(roverm > 0.5d0*ONEmTINY) then
                  rho(i,j,k) = rhocheck+drhodr*riso*(riso-M)/M
               else 
                  rho(i,j,k) = (rhocheck-drhodr*M/4.0)*(1.-cos(2.*M_PI*riso/M))/2.0
@@ -219,11 +234,11 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)
            Bvecz(i,j,k) = bondi_bmag*M**2*zhat/sqrt(det)/riso**2
            
            call Prim2ConGenM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k), &
-		  gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
-		  det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k), &
-		  tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k), &
-                  rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k), &
-		  eps(i,j,k),press(i,j,k), &
+                          gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
+                          det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k), &
+                    tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k), &
+                          rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k), &
+                          eps(i,j,k),press(i,j,k), &
                   Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),w_lorentz(i,j,k))
            
         end do