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

[rhaas at tapir.caltech.edu]

User: rhaas
Date: 2012/09/17 11:27 AM

Modified:
 /trunk/src/
  GRHydro_Bondi_new.F90

Log:
 GRHydro_InitData: Change debug messages; get rid of tabs;
 
 From: Bruno Coutinho Mundim <bcmsma at astro.rit.edu>

File Changes:

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

File [modified]: GRHydro_Bondi_new.F90
Delta lines: +40 -19
===================================================================
--- trunk/src/GRHydro_Bondi_new.F90	2012-09-17 16:27:22 UTC (rev 163)
+++ trunk/src/GRHydro_Bondi_new.F90	2012-09-17 16:27:24 UTC (rev 164)
@@ -51,8 +51,8 @@
   DECLARE_CCTK_FUNCTIONS
   
   CCTK_INT :: i, j, k, nx, ny, nz, imin, jb,N_points
-  CCTK_REAL :: ONEmTINY
-  PARAMETER (N_points=2000,ONEmTINY=0.999999d0)
+  CCTK_REAL :: ONEmTINY, tiny
+  PARAMETER (N_points=2000,ONEmTINY=0.999999d0,tiny=1.0d-14)
   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)
@@ -60,6 +60,8 @@
   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)
   Msq  = M*M
@@ -67,7 +69,9 @@
   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)
@@ -91,7 +95,9 @@
   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
@@ -134,10 +140,18 @@
      v_bondi(i)   = vtmp
   enddo
   
-  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)
+  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)
 
 !!$    // find the derivative near r=M
   rnew = 2.25 * M
@@ -152,7 +166,8 @@
   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)
@@ -170,7 +185,7 @@
            xhat = xp/riso
            yhat = yp/riso
            zhat = zp/riso
-	   roverm = riso/M
+           roverm = riso/M
            
            if(roverm > ONEmTINY) then
               rsch = 0.25 * ( 2.*riso + M)**2 / riso
@@ -201,16 +216,22 @@
            det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
 
            call Prim2ConGen(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),rho(i,j,k), &
-		  velx(i,j,k),vely(i,j,k),velz(i,j,k), &
-		  eps(i,j,k),press(i,j,k),w_lorentz(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),rho(i,j,k), &
+                velx(i,j,k),vely(i,j,k),velz(i,j,k), &
+                eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
            
-	   if(riso.gt.1.014d0.and.riso.lt.1.015)write(6,*)'Point to check:', &
-	          x(i,j,k),y(i,j,k),z(i,j,k),riso,gxx(i,j,k),dens(i,j,k),tau(i,j,k),&
-	          sx(i,j,k),sy(i,j,k),sz(i,j,k),rho(i,j,k),eps(i,j,k),&
-	          velx(i,j,k),vely(i,j,k),velz(i,j,k)	
+!           if(riso.gt.1.014d0.and.riso.lt.1.015) then 
+           if(abs(z(i,j,k)).lt.tiny.and. &
+              x(i,j,k).gt.6.7d0.and.x(i,j,k).lt.6.9.and. &
+              y(i,j,k).gt.1.34d0.and.y(i,j,k).lt.1.38 ) then 
+             write(debug_message,'(a,15f22.14)') 'Point to check:', &
+               x(i,j,k),y(i,j,k),z(i,j,k),riso,gxx(i,j,k),dens(i,j,k),&
+               tau(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),rho(i,j,k),eps(i,j,k),&
+               velx(i,j,k),vely(i,j,k),velz(i,j,k)
+             call CCTK_INFO(debug_message)
+           end if 
 
         end do
      end do