User: rhaas Date: 2012/09/17 11:27 AM
Modified: /trunk/ param.ccl /trunk/src/ GRHydro_Bondi_new.F90
Log: GRHydro_InitData: Move the metric fudge from Bondi routine to SmoothedPuncture thorn
From: Bruno Coutinho Mundim bcmsma@astro.rit.edu
Directory Changes:
Directory: /svn:executable/ ===========================
+ *
File Changes:
Directory: /trunk/src/ ======================
File [modified]: GRHydro_Bondi_new.F90 Delta lines: +11 -63 =================================================================== --- trunk/src/GRHydro_Bondi_new.F90 2012-09-17 16:27:20 UTC (rev 162) +++ trunk/src/GRHydro_Bondi_new.F90 2012-09-17 16:27:22 UTC (rev 163) @@ -67,7 +67,7 @@ rs = r_sonicpt_bondi gam = gl_gamma
- !!write(*,*) 'Bondi_pars:',M,mdot_sonicpt_bondi,r_sonicpt_bondi,gl_gamma + write(*,*) 'Bondi_pars:',M,mdot_sonicpt_bondi,r_sonicpt_bondi,gl_gamma
rmin_bondi = M * bondi_rmin(1) rmax_bondi = M * bondi_rmax(1) @@ -91,7 +91,7 @@ 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(*,*)'More pars:',cs,vs,rhos,hs,Kval,Qdot,logrmin,dlogr
rhotmp=1.0d30 imin=1 @@ -134,10 +134,10 @@ 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(*,*)"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 @@ -152,7 +152,7 @@ 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(6,*)'Rhocheck:',rhocheck,rhocheck2,drhodr
nx = cctk_lsh(1) ny = cctk_lsh(2) @@ -172,58 +172,6 @@ zhat = zp/riso roverm = riso/M
- if(bondi_use_smooth_puncture_data.ne.0 .and. roverm.lt.(bondi_smooth_puncture_radius*ONEmTINY)) then - - rsm=bondi_smooth_puncture_radius - !!write(6,*)'smoothing at a radius',rsm - - if(bondi_smooth_gxx.ne.0) then - - f=(1.d0+0.5d0/rsm)**4 - df=-1.d0*(2.d0*rsm+1.d0)**3/4.d0/rsm**5 - ddf=(2.d0*rsm+1.d0)**2*(4.d0*rsm+5.d0)/4.d0/rsm**6 - - a=f-0.625d0*rsm*df+0.125d0*rsm**2*ddf - b=0.75d0/rsm*df-0.25d0*ddf - c=-0.125d0/rsm**3*df+0.125d0/rsm**2*ddf - - !!write(6,*)'smooth gxx:',rsm,f,df,ddf,a,b,c - - gxx(i,j,k) = a + b*roverm**2 + c*roverm**4 - gyy(i,j,k) = a + b*roverm**2 + c*roverm**4 - gzz(i,j,k) = a + b*roverm**2 + c*roverm**4 - - else - f=1.d0+0.5d0/rsm - df=-0.5d0/rsm**2 - ddf=1.d0/rsm**3 - - a=f-0.625d0*rsm*df+0.125d0*rsm**2*ddf - b=0.75d0/rsm*df-0.25d0*ddf - c=-0.125d0/rsm**3*df+0.125d0/rsm**2*ddf - - !!write(6,*)'smooth psi:',rsm,f,df,ddf,a,b,c - - psinew=a + b*roverm**2 + c*roverm**4 - gxx(i,j,k) = psinew**4 - gyy(i,j,k) = psinew**4 - gzz(i,j,k) = psinew**4 - endif - gxy(i,j,k) = 0.0 - gxz(i,j,k) = 0.0 - gyz(i,j,k) = 0.0 - - else 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/sqrt(gxx(i,j,k)) - if(roverm > ONEmTINY) then rsch = 0.25 * ( 2.*riso + M)**2 / riso jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) @@ -259,10 +207,10 @@ 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)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)
end do end do
Property changes on: trunk/src/GRHydro_Bondi_new.F90 ___________________________________________________________________
Directory: /trunk/ ==================
File [modified]: param.ccl Delta lines: +0 -14 =================================================================== --- trunk/param.ccl 2012-09-17 16:27:20 UTC (rev 162) +++ trunk/param.ccl 2012-09-17 16:27:22 UTC (rev 163) @@ -99,20 +99,6 @@ 0.0:* :: "Anything positive" } 1.0
-BOOLEAN bondi_use_smooth_puncture_data "Smooth puncture solution inside R=M using 4th order matching polynomial" -{ -} "no" - - -BOOLEAN bondi_smooth_gxx "Smooth puncture data by matching to gxx not psi" -{ -} "yes" - -REAL bondi_smooth_puncture_radius "Radius inside which to mathc functions, in units of R/M" -{ - 0:0.5 :: "Less than/equal to 0.5" -} 0.5 - BOOLEAN change_shock_direction "Change the shock direction" { } "no"
Property changes on: trunk/param.ccl ___________________________________________________________________