From a0d787cd49f3782f29b231ae26f51ed1d1c0642a Mon Sep 17 00:00:00 2001 From: rhaas Date: Fri, 9 Nov 2012 01:59:09 +0000 Subject: GRHydro_InitData: port smoother velocity matching to mag. bondi routine. From: Bruno Coutinho Mundim git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@179 ac85fae7-cede-4708-beff-ae01c7fa1c26 --- src/GRHydro_BondiM_new.F90 | 53 +++++++++++++++++++++++++++++++++++----------- src/GRHydro_Bondi_new.F90 | 2 +- 2 files changed, 42 insertions(+), 13 deletions(-) diff --git a/src/GRHydro_BondiM_new.F90 b/src/GRHydro_BondiM_new.F90 index 30c1a8f..325e158 100644 --- a/src/GRHydro_BondiM_new.F90 +++ b/src/GRHydro_BondiM_new.F90 @@ -64,7 +64,7 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS) 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 + CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm,dudr,uisocheck2,auiso,buiso character(400) :: debug_message @@ -118,7 +118,8 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS) endif enddo - rhotmp = -1. !!$ start with guess +!!$ rhotmp = -1. !!$ start with guess + rhotmp=rhos !!$ start with value at sonic point! do i=imin,N_points rspher = r_bondi(i) @@ -132,7 +133,8 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS) v_bondi(i) = vtmp end do - rhotmp = -1. +!!$ rhotmp = -1. + rhotmp=rhos !!$ start with value at sonic point! do i=imin-1,1,-1 rspher = r_bondi(i) @@ -165,27 +167,42 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS) !!$ 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) - -!!$ // find the derivative near r=M +!!$ // find the derivative near r=M in isotropic coords = r=9/4M in schwarzschild; rnew = 2.25 * M - j = floor ( 0.5 + (log10(rnew) - logrmin) / dlogr ) + j = floor ((log10(rnew) - logrmin) / dlogr + 1.0) +!!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck = rho_bondi(j) call find_bondi_solution(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot ) uisocheck = 4.0*vcheck/3.0 +!!$ the previous point was r=M in isotropic coords = r=9/4M in schwarzschild; this one is r=1.01M in isotropic rnew = 0.25 * 3.02**2 * M/1.01 - j = floor( 0.5 + (log10(rnew) - logrmin) / dlogr ) + j = floor((log10(rnew) - logrmin) / dlogr + 1.0) +!!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck2 = rho_bondi(j) call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot ) + uisocheck2 = vcheck2 / (1.0 - 1.0/2.02) / (1.0+ 1.0/2.02) + drhodr = 100.0*(rhocheck2-rhocheck)/M +!!$ Don't divide by M here, to simplify the math + dudr = 100.0*(uisocheck2-uisocheck) + write(debug_message,'(a,3f22.14)') 'Rhocheck:',rhocheck,rhocheck2,drhodr call CCTK_INFO(debug_message) + write(debug_message,'(a,3f22.14)') 'Ucheck:',uisocheck,uisocheck2,dudr + call CCTK_INFO(debug_message) nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) + + + write(debug_message,'(a,3f22.14)') 'Lower bound coordinates:',x(1,1,1),y(1,1,1),z(1,1,1) + call CCTK_INFO(debug_message) + write(debug_message,'(a,3f22.14)') 'Upper bound coordinates:',x(nx,ny,nz),y(nx,ny,nz),z(nx,ny,nz) + call CCTK_INFO(debug_message) do i=1,nx do j=1,ny @@ -203,12 +220,17 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS) if(roverm > ONEmTINY) then rsch = 0.25 * ( 2.*riso + M)**2 / riso - jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) +!!$ jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) + jb = floor( (log10(rsch) - logrmin) / dlogr + 1.0) +!!$ jb = NINT((log10(rsch) - logrmin) / dlogr + 1.0) - if(jb > N_points)jb = N_points + if(jb >= N_points) jb = N_points-1 - rhotmp = rho_bondi(jb) - call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot) + rhotmp = rho_bondi(jb)+(rho_bondi(jb+1)-rho_bondi(jb))*& + (rsch-r_bondi(jb))/(r_bondi(jb+1)-r_bondi(jb)) +!!$ call find_bondi_solution_bracket( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot, & +!!$ rho_bondi(jb),rho_bondi(jb+1)) + call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot) rho(i,j,k) = rhotmp uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso) else @@ -218,7 +240,14 @@ subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS) rho(i,j,k) = (rhocheck-drhodr*M/4.0)*(1.-cos(2.*M_PI*riso/M))/2.0 endif utmp = Kval * rho(i,j,k)**( gam ) / (gam - 1.) - uiso = uisocheck * riso / M + +!!$ match to uiso and dudr at roverm=1 +!!$ a R + b R^3 ---> a+b = uisocheck; a+3b = dudr +!!$ b = (dudr-uisocheck)/2; a=3*uisocheck-dudr)/2 + + auiso = 1.5*uisocheck-0.5*dudr + buiso = 0.5*dudr-0.5*uisocheck + uiso = roverm*(auiso+buiso*roverm**2) endif eps(i,j,k) = utmp/rhotmp diff --git a/src/GRHydro_Bondi_new.F90 b/src/GRHydro_Bondi_new.F90 index a0077c2..b31ba7c 100644 --- a/src/GRHydro_Bondi_new.F90 +++ b/src/GRHydro_Bondi_new.F90 @@ -148,7 +148,7 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) 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) -- cgit v1.2.3