From 2651028fcf5d3be767113294ab418e964b7c496a Mon Sep 17 00:00:00 2001 From: rhaas Date: Mon, 17 Sep 2012 16:27:24 +0000 Subject: GRHydro_InitData: Change debug messages; get rid of tabs; From: Bruno Coutinho Mundim git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@164 ac85fae7-cede-4708-beff-ae01c7fa1c26 --- src/GRHydro_Bondi_new.F90 | 59 ++++++++++++++++++++++++++++++++--------------- 1 file changed, 40 insertions(+), 19 deletions(-) diff --git a/src/GRHydro_Bondi_new.F90 b/src/GRHydro_Bondi_new.F90 index b539c5f..63243eb 100755 --- a/src/GRHydro_Bondi_new.F90 +++ b/src/GRHydro_Bondi_new.F90 @@ -51,8 +51,8 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) 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 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) 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 @@ subroutine GRHydro_Bondi_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) @@ -91,7 +95,9 @@ subroutine GRHydro_Bondi_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 @@ -134,10 +140,18 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) 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 @@ subroutine GRHydro_Bondi_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) @@ -170,7 +185,7 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) 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 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) 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 -- cgit v1.2.3