From 61eeb17d64b2928122ecbbb410fd3946f20b925a Mon Sep 17 00:00:00 2001 From: rhaas Date: Mon, 17 Sep 2012 16:27:20 +0000 Subject: GRHydro_InitData: Bondi routine: smooth puncture data. From: Bruno Coutinho Mundim git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@162 ac85fae7-cede-4708-beff-ae01c7fa1c26 --- param.ccl | 12 +++++++- src/GRHydro_Bondi_new.F90 | 72 ++++++++++++++++++++++++++++++++++++----------- 2 files changed, 67 insertions(+), 17 deletions(-) diff --git a/param.ccl b/param.ccl index 07b6db2..b3c050f 100644 --- a/param.ccl +++ b/param.ccl @@ -99,10 +99,20 @@ REAL shock_radius "Radius of sperical shock" 0.0:* :: "Anything positive" } 1.0 -BOOLEAN use_smooth_puncture_data "Smooth puncture solution inside R=M using 4th order matching polynomial" +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 index 9e56ef2..9b9d45e 100644 --- a/src/GRHydro_Bondi_new.F90 +++ b/src/GRHydro_Bondi_new.F90 @@ -51,12 +51,14 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) DECLARE_CCTK_FUNCTIONS CCTK_INT :: i, j, k, nx, ny, nz, imin, jb,N_points - PARAMETER (N_points=2000) + CCTK_REAL :: ONEmTINY + PARAMETER (N_points=2000,ONEmTINY=0.999999d0) 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) CCTK_REAL :: drhodr, det, rhocheck, rhocheck2, riso, rnew, rsch, ucheck, psinew CCTK_REAL :: uiso, uisocheck, vcheck, ucheck2, vcheck2, xhat,yhat, zhat, xp, yp, zp + CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm !!$set_bondi_parameters M = bondi_central_mass(1) @@ -65,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) @@ -89,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 @@ -132,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 @@ -150,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) @@ -168,12 +170,45 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) xhat = xp/riso yhat = yp/riso zhat = zp/riso + roverm = riso/M - if(use_smooth_puncture_data.ne.0 .and. riso.lt.M) then - psinew=2.875 - 5*(riso/M)**2 + 6*(riso/M)**4 - gxx(i,j,k) = psinew**4 - gyy(i,j,k) = psinew**4 - gzz(i,j,k) = psinew**4 + 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 @@ -187,9 +222,9 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) gyz(i,j,k) = 0.0 endif - alp(i,j,k) = 1.0/gxx(i,j,k)**2 + alp(i,j,k) = 1.0/sqrt(gxx(i,j,k)) - if(riso > M) then + if(roverm > ONEmTINY) then rsch = 0.25 * ( 2.*riso + M)**2 / riso jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) @@ -200,7 +235,7 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) rho(i,j,k) = rhotmp uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso) else - if(riso > 0.5*M) then + if(roverm > 0.5d0*ONEmTINY) then rho(i,j,k) = rhocheck+drhodr*riso*(riso-M)/M else rho(i,j,k) = (rhocheck-drhodr*M/4.0)*(1.-cos(2.*M_PI*riso/M))/2.0 @@ -224,6 +259,11 @@ 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) + end do end do end do -- cgit v1.2.3