From 338d9065455d7a73b219f03c4d4c49d19cf8a556 Mon Sep 17 00:00:00 2001 From: rhaas Date: Fri, 9 Nov 2012 01:59:00 +0000 Subject: Add find_bondi_solution_bracket and velocity smoothing from Josh's patch. From: Bruno Coutinho Mundim git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@176 ac85fae7-cede-4708-beff-ae01c7fa1c26 --- src/GRHydro_Bondi_new.F90 | 200 +++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 188 insertions(+), 12 deletions(-) diff --git a/src/GRHydro_Bondi_new.F90 b/src/GRHydro_Bondi_new.F90 index 6327080..e161043 100644 --- a/src/GRHydro_Bondi_new.F90 +++ b/src/GRHydro_Bondi_new.F90 @@ -58,7 +58,7 @@ subroutine GRHydro_Bondi_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 @@ -112,7 +112,8 @@ subroutine GRHydro_Bondi_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) @@ -126,7 +127,8 @@ subroutine GRHydro_Bondi_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) @@ -162,23 +164,39 @@ subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) !!$ // find the derivative near r=M rnew = 2.25 * M - j = floor ( 0.5 + (log10(rnew) - logrmin) / dlogr ) + j = floor ((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 ) + call find_bondi_solution_bracket(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot, & + rho_bondi(j),rho_bondi(j+1)) 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 rhocheck2 = rho_bondi(j) - call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot ) + call find_bondi_solution_bracket( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot, & + rho_bondi(j),rho_bondi(j+1)) + 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 @@ -196,12 +214,14 @@ subroutine GRHydro_Bondi_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 - 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)) rho(i,j,k) = rhotmp uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso) else @@ -211,7 +231,14 @@ subroutine GRHydro_Bondi_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 @@ -382,3 +409,152 @@ subroutine find_bondi_solution(r, rho, u, v, rs, rhos, M, Mdot, Kval, gam, Qdot end subroutine find_bondi_solution + subroutine find_bondi_solution_bracket(r, rho, u, v, rs, rhos, M, Mdot, Kval, gam, Qdot,rho1,rho2 ) + implicit none + + CCTK_REAL rho1,rho2 + CCTK_REAL :: r, rho, u, v, rs, rhos, M, Mdot, Kval, gam, Qdot + CCTK_REAL :: ur,r_sol, rho_old + CCTK_REAL :: f, df, dx, x_old, resid, jac + CCTK_REAL :: errx, x_orig, x + CCTK_INT :: n_iter, i_extra, doing_extra, keep_iterating, i_increase + CCTK_REAL :: vp, h, hp, term + + CCTK_REAL :: newt_tol_b,small_bondi + CCTK_INT :: max_newt_iter_b, extra_newt_iter_b + + if(rho.gt.rho1.or.rho.lt.rho2) then + write(6,*)'Very Bad rho! ',rho1,rho2,rho + stop + endif + + max_newt_iter_b = 30 + newt_tol_b = 1.0e-15 + extra_newt_iter_b = 2 + small_bondi = 1.0e-20 + +!!$ if(r>8.1043 .and. r<8.1044)write(*,*)'init guess:',r,rho +!!$ write(*,*)'init guess:',r,rho + +!!$ if (rho < 0.) then +!!$ if( r > 0.9*rs .and. r < 1.1*rs ) then +!!$ rho = rhos +!!$ else +!!$ if(r < rs) then +!!$ ur = r**(-0.5) +!!$ else +!!$ ur = 0.5*r**(-1.5) +!!$ endif +!!$ rho = Mdot / (4.*M_PI * r * r * ur) +!!$ endif +!!$ endif + +!!$ if(r>8.1043 .and. r<8.1044) +!!$ write(*,*)'init guess:',r,ur,rho,rs,rhos,Mdot +!!$ if(r<1.0001e-10)write(*,*)'init guess:',r,ur,rho,rs,rhos,Mdot +!!$ write(*,*)'init guess:',r,ur,rho,rs,rhos,Mdot + +!!$ set global variables needed by residual function: + + r_sol = r + + +!!$ Use Newton's method to find rho: +!!$ gnr_bondi( rho, NEWT_DIM_B, bondi_resid) + errx = 1.0 + df=1.0 + f=1.0 + doing_extra = 0 + + rho_old=rho + x=rho + + n_iter = 0 + +!!$ Start the Newton-Raphson iterations : + keep_iterating = 1 + do while( keep_iterating == 1 ) + + hp = Kval * gam * x**(gam - 2.) !!$ // dh/drho + h = 1. + hp * x / ( gam - 1. ) + v = Mdot / ( 4. * M_PI * r_sol * r_sol * x ) + vp = -v / x !!$ // dv/drho + term = 1. - 2.*M/r_sol + v*v + resid = -Qdot + h * h * term + jac = 2. * h *( hp*term + h*v*vp ) + dx = -resid / jac + f = 0.5*resid**2 + df = -2. * f + + errx = 0. + x_old = x + + x=x+dx + + if(x.gt.rho1.or.x.lt.rho2) then +!!$ write(6,*)'Bad rho! ',rho1,rho2,rho + if(x.gt.rho1)x=0.5*(x_old+rho1) + if(x.lt.rho2)x=0.5*(x_old+rho2) + endif + + if(x.gt.rho1.or.x.lt.rho2) then + write(6,*)'Bad rho, bad! ',rho1,rho2,rho + stop + endif + +!!$ if(r>8.1043 .and. r<8.1044)write(*,*)'iter:',x,dx,resid,jac + +!!$ /****************************************/ +!!$ /* Calculate the convergence criterion */ +!!$ /****************************************/ + +!!$ /* For the new criterion, always look at relative error in indep. variable: */ +!!$ // METHOD specific: + if(x==0) then + errx = abs(dx) + x = small_bondi + else + errx = abs(dx/x) + endif + +!!$ /*****************************************************************************/ +!!$ /* If we've reached the tolerance level, then just do a few extra iterations */ +!!$ /* before stopping */ +!!$ /*****************************************************************************/ + +!!$ if(r>8.1043 .and. r<8.1044)write(*,*)'iter3:',errx,newt_tol_b,keep_iterating, & +!!$ doing_extra,i_extra + + if((abs(errx)<=newt_tol_b) .and. (doing_extra == 0) .and. (extra_newt_iter_b > 0)) & + doing_extra=1 + + if( doing_extra == 1 ) i_extra=i_extra+1 + + if( ((abs(errx) <= newt_tol_b).and.(doing_extra == 0)) .or. & + (i_extra > extra_newt_iter_b) .or. (n_iter >= (max_newt_iter_b-1)) ) then + keep_iterating = 0 + if(n_iter>=(max_newt_iter_b-1))write(6,*)'Extra iterations!' + endif + +!!$ if(r>8.1043 .and. r<8.1044)write(*,*)'iter4:',errx,newt_tol_b,keep_iterating, & +!!$ doing_extra,i_extra + + + n_iter=n_iter+1 + + end do + + rho=x + +!!$ Calculate other quantities: + u = Kval * rho**( gam ) / (gam - 1.) + v = Mdot / ( 4. * M_PI * r * r * rho ) + +!!$ if(r>8.1043 .and. r<8.1044) +!!$ write(*,*)'final:',r,rho,u,v + + + return + + end subroutine find_bondi_solution_bracket + -- cgit v1.2.3