[Commits] [svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 162)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Mon Sep 17 11:27:20 CDT 2012
User: rhaas
Date: 2012/09/17 11:27 AM
Modified:
/trunk/
param.ccl
/trunk/src/
GRHydro_Bondi_new.F90
Log:
GRHydro_InitData: Bondi routine: smooth puncture data.
From: Bruno Coutinho Mundim <bcmsma at astro.rit.edu>
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Bondi_new.F90
Delta lines: +56 -16
===================================================================
--- trunk/src/GRHydro_Bondi_new.F90 2012-09-17 16:27:18 UTC (rev 161)
+++ trunk/src/GRHydro_Bondi_new.F90 2012-09-17 16:27:20 UTC (rev 162)
@@ -51,12 +51,14 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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
Directory: /trunk/
==================
File [modified]: param.ccl
Delta lines: +11 -1
===================================================================
--- trunk/param.ccl 2012-09-17 16:27:18 UTC (rev 161)
+++ trunk/param.ccl 2012-09-17 16:27:20 UTC (rev 162)
@@ -99,10 +99,20 @@
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"
More information about the Commits
mailing list