[Commits] [svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 172)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Mon Sep 17 11:27:39 CDT 2012
User: rhaas
Date: 2012/09/17 11:27 AM
Modified:
/trunk/src/
GRHydro_BondiM_new.F90
Log:
GRHydro_InitData: port bondi fixes to magnetized version.
From: Bruno Coutinho Mundim <bcmsma at astro.rit.edu>
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_BondiM_new.F90
Delta lines: +41 -26
===================================================================
--- trunk/src/GRHydro_BondiM_new.F90 2012-09-17 16:27:37 UTC (rev 171)
+++ trunk/src/GRHydro_BondiM_new.F90 2012-09-17 16:27:39 UTC (rev 172)
@@ -57,13 +57,17 @@
DECLARE_CCTK_FUNCTIONS
CCTK_INT :: i, j, k, nx, ny, nz, imin, jb,N_points
- PARAMETER (N_points=2000)
+ CCTK_REAL :: ONEmTINY, tiny
+ PARAMETER (N_points=100000,ONEmTINY=0.999999d0,tiny=1.0d-12)
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
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
@@ -71,7 +75,9 @@
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)
@@ -95,7 +101,9 @@
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
@@ -137,11 +145,27 @@
u_bondi(i) = utmp
v_bondi(i) = vtmp
enddo
+
+ open (47,file="bondi.asc",form="formatted")
+ do i=1,N_points
+ write(47,'(i5,4f22.14)')i,r_bondi(i),rho_bondi(i),&
+ u_bondi(i),v_bondi(i)
+ end do
+ close(47)
+
+!!$ 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)
- 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
@@ -156,7 +180,8 @@
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)
@@ -174,19 +199,9 @@
xhat = xp/riso
yhat = yp/riso
zhat = zp/riso
+ roverm = riso/M
- 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/gxx(i,j,k)**2
-
- if(riso > M) then
+ if(roverm > ONEmTINY) then
rsch = 0.25 * ( 2.*riso + M)**2 / riso
jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr )
@@ -197,7 +212,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
@@ -219,11 +234,11 @@
Bvecz(i,j,k) = bondi_bmag*M**2*zhat/sqrt(det)/riso**2
call Prim2ConGenM(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),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(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), &
+ 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),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(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), &
Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),w_lorentz(i,j,k))
end do
More information about the Commits
mailing list