From 91a2ee778dc93e45a759810d790d04e9dfcbdd04 Mon Sep 17 00:00:00 2001 From: rhaas Date: Mon, 17 Sep 2012 16:27:39 +0000 Subject: GRHydro_InitData: port bondi fixes to magnetized version. From: Bruno Coutinho Mundim git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@172 ac85fae7-cede-4708-beff-ae01c7fa1c26 --- src/GRHydro_BondiM_new.F90 | 67 ++++++++++++++++++++++++++++------------------ 1 file changed, 41 insertions(+), 26 deletions(-) diff --git a/src/GRHydro_BondiM_new.F90 b/src/GRHydro_BondiM_new.F90 index 3b1f47d..30c1a8f 100644 --- a/src/GRHydro_BondiM_new.F90 +++ b/src/GRHydro_BondiM_new.F90 @@ -57,12 +57,16 @@ subroutine GRHydro_BondiM_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, 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) @@ -71,7 +75,9 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS) 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 @@ subroutine GRHydro_BondiM_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(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 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS) 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 @@ subroutine GRHydro_BondiM_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(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 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS) 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 @@ subroutine GRHydro_BondiM_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 @@ -219,11 +234,11 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS) 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 -- cgit v1.2.3