From 2f159f13dfb98c6c906118cd5eb5b861ec9d0cd8 Mon Sep 17 00:00:00 2001 From: rhaas Date: Mon, 17 Sep 2012 16:27:22 +0000 Subject: GRHydro_InitData: Move the metric fudge from Bondi routine to SmoothedPuncture thorn From: Bruno Coutinho Mundim git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@163 ac85fae7-cede-4708-beff-ae01c7fa1c26 --- param.ccl | 14 --------- src/GRHydro_Bondi_new.F90 | 74 +++++++---------------------------------------- 2 files changed, 11 insertions(+), 77 deletions(-) mode change 100644 => 100755 param.ccl mode change 100644 => 100755 src/GRHydro_Bondi_new.F90 diff --git a/param.ccl b/param.ccl old mode 100644 new mode 100755 index b3c050f..466fbe0 --- a/param.ccl +++ b/param.ccl @@ -99,20 +99,6 @@ REAL shock_radius "Radius of sperical shock" 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" diff --git a/src/GRHydro_Bondi_new.F90 b/src/GRHydro_Bondi_new.F90 old mode 100644 new mode 100755 index 9b9d45e..b539c5f --- a/src/GRHydro_Bondi_new.F90 +++ b/src/GRHydro_Bondi_new.F90 @@ -67,7 +67,7 @@ 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(*,*) '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 @@ 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(*,*)'More pars:',cs,vs,rhos,hs,Kval,Qdot,logrmin,dlogr rhotmp=1.0d30 imin=1 @@ -134,10 +134,10 @@ 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(*,*)"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 @@ 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(6,*)'Rhocheck:',rhocheck,rhocheck2,drhodr nx = cctk_lsh(1) ny = cctk_lsh(2) @@ -172,58 +172,6 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) 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 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) 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 -- cgit v1.2.3