diff options
Diffstat (limited to 'src/GRHydro_Bondi_new.F90')
| -rw-r--r-- | src/GRHydro_Bondi_new.F90 | 359 |
1 files changed, 359 insertions, 0 deletions
diff --git a/src/GRHydro_Bondi_new.F90 b/src/GRHydro_Bondi_new.F90 new file mode 100644 index 0000000..31e7c43 --- /dev/null +++ b/src/GRHydro_Bondi_new.F90 @@ -0,0 +1,359 @@ + /*@@ + @file GRHydro_Bondi.F90 + @date Wed Jan 13 13:00:49 EST 2010 + @author Scott C. Noble + @desc + Hydro initial data for the relativistic Bondi solution about + a single Schwarzschild black hole. + @enddesc + @@*/ + +/* + Calculates the Bondi solution, or the spherically symmetric hydrostationary + solution to a fluid on a static fixed background spacetime. We assume that one can + calculate a radius "r" from the grid and that with respect to this radial coordinate, + the solution satisfies + + d (\rho u^r) / dr = 0 + + Assumes that the equation of state is P = K \rho^\Gamma and K is set by + the location of the sonic point. + + + -- Implicitly assumes that there is no spin in the geometry as there is no Bondi + solution for spinning black holes. If a spin is specified, a spherically symmetric + is still assumed but the 4-velocity is set consistently with the spinning spacetime. +*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" +#include "GRHydro_Macros.h" + +# define M_PI 3.14159265358979323846 /* pi */ + +!!$Newton-Raphson parameters: + +#define velx(i,j,k) vel(i,j,k,1) +#define vely(i,j,k) vel(i,j,k,2) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) + +subroutine GRHydro_Bondi_Iso(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, nx, ny, nz, imin, jb,N_points + PARAMETER (N_points=2000) + 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 + + !!$set_bondi_parameters + M = bondi_central_mass(1) + Msq = M*M + Mdot = mdot_sonicpt_bondi + rs = r_sonicpt_bondi + gam = 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) + + cs_sq = M / ( 2.*rs - 3.*M ) + + if( cs_sq > (gam - 1.)) then + cs_sq = gam - 1. + rs = 0.5 * M * ( 3. + 1./cs_sq ) + endif + + cs = sqrt(cs_sq) + vs_sq = M / ( 2. * rs ) + vs = sqrt(vs_sq) + rhos = Mdot / ( 4. * M_PI * vs * rs * rs ) + gtemp = gam - 1. + hs = 1. / ( 1. - cs_sq / (gam - 1.) ) + Kval = hs * cs_sq * rhos**(-gtemp) / gam + Qdot = hs * hs * ( 1. - 3. * vs_sq ) + + logrmin = log10(rmin_bondi) + dlogr = (log10(rmax_bondi) - logrmin)/(1.*(N_points-1)) + + write(*,*)'More pars:',cs,vs,rhos,hs,Kval,Qdot,logrmin,dlogr + + rhotmp=1.0d30 + imin=1 + + do i=1,N_points + logr_bondi(i) = logrmin + dlogr*(i-1) + r_bondi(i) = 10.**(logr_bondi(i)) + utmp = abs(r_bondi(i) - r_sonicpt_bondi) + if (utmp < rhotmp) then + rhotmp = utmp + imin = i + endif + enddo + + rhotmp = -1. !!$ start with guess + + do i=imin,N_points + rspher = r_bondi(i) + call find_bondi_solution( rspher, rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot ) + if(rhotmp < initial_rho_abs_min) then + rhotmp = initial_rho_abs_min + utmp = Kval * rhotmp**gl_gamma / (gl_gamma - 1.) + endif + rho_bondi(i) = rhotmp + u_bondi(i) = utmp + v_bondi(i) = vtmp + end do + + rhotmp = -1. + + do i=imin-1,1,-1 + rspher = r_bondi(i) + call find_bondi_solution( rspher, rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot ) + if(rhotmp < initial_rho_abs_min) then + rhotmp = initial_rho_abs_min + utmp = K * rhotmp**gl_gamma / (gl_gamma - 1.) + endif + rho_bondi(i) = rhotmp + u_bondi(i) = utmp + 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) + +!!$ // find the derivative near r=M + rnew = 2.25 * M + j = floor ( 0.5 + (log10(rnew) - logrmin) / dlogr ) + 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 + + rnew = 0.25 * 3.02**2 * M/1.01 + j = floor( 0.5 + (log10(rnew) - logrmin) / dlogr ) + rhocheck2 = rho_bondi(j) + 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 + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + do i=1,nx + do j=1,ny + do k=1,nz + + xp=x(i,j,k) + yp=y(i,j,k) + zp=z(i,j,k) + + riso = sqrt(xp*xp + yp*yp + zp*zp +1.0e-16) + xhat = xp/riso + yhat = yp/riso + zhat = zp/riso + + 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 + rsch = 0.25 * ( 2.*riso + M)**2 / riso + jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) + + if(jb > N_points)jb = N_points + + rhotmp = rho_bondi(jb) + 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 + if(riso > 0.5*M) 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 + endif + utmp = Kval * rho(i,j,k)**( gam ) / (gam - 1.) + uiso = uisocheck * riso / M + endif + eps(i,j,k) = utmp/rhotmp + + w_lorentz(i,j,k) = sqrt(1.0+gxx(i,j,k) * uiso**2) + velx(i,j,k) = -1.0*uiso * xhat / w_lorentz(i,j,k) + vely(i,j,k) = -1.0*uiso * yhat / w_lorentz(i,j,k) + velz(i,j,k) = -1.0*uiso * zhat / w_lorentz(i,j,k) + + det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + + call Prim2ConGen(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),rho(i,j,k), & + 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)) + + end do + end do + end do + + densrhs = 0.d0 + srhs = 0.d0 + taurhs = 0.d0 + + return + +end subroutine GRHydro_Bondi_Iso + +subroutine find_bondi_solution(r, rho, u, v, rs, rhos, M, Mdot, Kval, gam, Qdot ) + implicit none + + 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 + + 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 + + /* Save old values before calculating the new: */ + errx = 0. + x_old = x + + x=x+dx + +!!$ 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)) ) & + keep_iterating = 0 + +!!$ 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 + |