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