diff options
Diffstat (limited to 'src/GRHydro_TOV.F90')
| -rw-r--r-- | src/GRHydro_TOV.F90 | 324 |
1 files changed, 324 insertions, 0 deletions
diff --git a/src/GRHydro_TOV.F90 b/src/GRHydro_TOV.F90 new file mode 100644 index 0000000..f40af9a --- /dev/null +++ b/src/GRHydro_TOV.F90 @@ -0,0 +1,324 @@ + /*@@ + @file GRHydro_TOV.F90 + @date Wed Feb 13 02:53:25 2002 + @author Scott Hawley + @desc + Initial data of the TOV type. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + /*@@ + @routine GRHydro_TOV + @date Sat Jan 26 02:53:49 2002 + @author Scott Hawley + @desc + Initial data for TOV stars. + @enddesc + @calls + @calledby + @history + Alteration of the GRHydro_shockwave written by Ian Hawke + @endhistory + +@@*/ + +subroutine GRHydro_TOV(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + real*8 rho0 + common / com_tov / rho0 + + integer nr, n_succ + integer i,j,k,nx,ny,nz + integer CCTK_Equals + CCTK_REAL rmax, dr + CCTK_REAL xlb, xub + CCTK_REAL ylb, yub + CCTK_REAL zlb, zub + CCTK_REAL tol + parameter ( tol = 1e-10 ) + + CCTK_REAL, allocatable, dimension (:) :: gtt, grr + + CCTK_REAL det + integer rc + +#include "EOS_Base.inc" + +#ifdef HAVE_ODEPACK + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +!----------------------------------------------------------- +! Set up the 1D spherical grid +!----------------------------------------------------------- +! decide maximum radius to integrate out to: +! Just use the maximum diagonal width of the grid +! BUG: this may not be the proper procedure call + call CCTK_CoordRange(cctkGH, xlb, xub, -1, "x", "cart3d") + call CCTK_CoordRange(cctkGH, ylb, yub, -1, "y", "cart3d") + call CCTK_CoordRange(cctkGH, zlb, zub, -1, "z", "cart3d") + rmax = sqrt( (xub-xlb)**2 + (yub-ylb)**2 + (zub-zlb)**2 ) + + nr = tov_nr + dr = rmax / (nr - 1) + do j = 1, nr + tov_r(j) = (j-1) * dr + enddo + + +!----------------------------------------------------------- +! Assign data at r=0 +!----------------------------------------------------------- + rho0 = GRHydro_rho_central + tov_rho(1) = rho0 + tov_alpha(1) = 1.0 + tov_m(1) = 0.0 + +!----------------------------------------------------------- +! Get the 1D TOV data +!----------------------------------------------------------- + call tov_lsoda (tov_rho, tov_p, tov_m, tov_alpha, & + tov_r, nr, tol , n_succ, rc) + +!----------------------------------------------------------- +! Handle any errors from integration routine +!----------------------------------------------------------- + if (rc .ne. 0) then + call CCTK_WARN(0,"Fatal error in TOV integration") + endif + if (n_succ .lt. nr) then + call CCTK_WARN(0,"TOV integration did not extend as far out as requested") + endif + nr = n_succ + +!----------------------------------------------------------- +! fill in additional physical info +!----------------------------------------------------------- + do j=1,nr +! BUG: the call looks something like this... + tov_rho(j) = EOS_RestMassDens(GRHydro_eos_handle,tov_p(j)) + end do + +!----------------------------------------------------------- +! Generate 3D data from 1D data +! BUG: These Spin1D functions do not yet exist! +!----------------------------------------------------------- + call Spin1D(tov_rho, tov_r, nr, rho, x, y, z, nx, ny, nz) + call Spin1D(tov_p, tov_r, nr, press, x, y, z, nx, ny, nz) + call Spin1D(tov_rho, tov_r, nr, eps, x, y, z, nx, ny, nz) + + allocate(grr(nr), gtt(nr)) + do i = 1, nr + grr(i) = tov_r(i) / (tov_r(i) - 2*tov_m(i)) + gtt(i) = - tov_alpha(i)**2 + enddo + call Spin1DMetric(tov_m, tov_r, nr, gxx, gxy, gxz, gyy, gyz, gzz, & + x, y, z, nx, ny, nz) + + call Spin1DMetric(grr, gtt, gxx, gxy, gxz, gyy, gyz, gzz, nx, ny, nz, & + x, y, z, r) + deallocate(gtt) + deallocate(grr) + + do i=1,nx + do j=1,ny + do k=1,nz + velx(i,j,k) = 0.0 + vely(i,j,k) = 0.0 + velz(i,j,k) = 0.0 + w_lorentz(i,j,k) = 1.0 + end do + end do + enddo + +!----------------------------------------------------------- +! Set additional variables +!----------------------------------------------------------- +! Conserved Variables + det = 1.0d0 + call prim2con(GRHydro_eos_handle,gxx, gxy, gxz, gyy, gyz, gzz, det, & + dens, sx, sy, sz, tau, rho, & + velx, vely, velz, eps, press, w_lorentz) + +! Copy data to other time steps + do i=1,nx + do j=1,ny + do k=1,nz + +! Primitive variables + rho_p(i,j,k) = rho(i,j,k) + press_p(i,j,k) = press(i,j,k) + eps_p(i,j,k) = eps(i,j,k) + velx_p(i,j,k) = velx(i,j,k) + vely_p(i,j,k) = vely(i,j,k) + velz_p(i,j,k) = velz(i,j,k) + w_lorentz_p(i,j,k) = w_lorentz(i,j,k) + +! Conserved variables + dens_p(i,j,k) = dens(i,j,k) + tau_p(i,j,k) = tau(i,j,k) + sx_p(i,j,k) = sx(i,j,k) + sy_p(i,j,k) = sy(i,j,k) + sz_p(i,j,k) = sz(i,j,k) + + end do + end do + end do + +! BUG: Why am I doing this? + densrhs = 0.d0 + sxrhs = 0.d0 + syrhs = 0.d0 + szrhs = 0.d0 + taurhs = 0.d0 + + return + +end subroutine GRHydro_TOV + + + +! ======================================================================= + + subroutine Spin1D(rad_func, r, nr, cart_func, x, y, z, & + nx, ny, nz) +! ................................................................ +! Takes a "radial grid function" and creates a 3D grid function +! BUG/FEATURE: This doesnt do any fancy interpolation at all! +! + implicit none + + integer nr, nx, ny, nz + CCTK_REAL rad_func(nr), r(nr) + CCTK_REAL cart_func(nx,ny,nz), x(nx,ny,nz), y(nx,ny,nz), z(nx,ny,nz) + + integer i,j,k, ri + real rloc + + do k = 1, nz + do j = 1, ny + do i = 1, nx + +! Given a point in 3D cartesian space, substitute next-further-out data +! value in 1D spherical space + rloc = sqrt(x(i,j,k)**2 + y(i,j,k)**2 + z(i,j,k)**2) + ri = 1 +100 continue + if (r(ri) .lt. rloc) then + ri = ri+1 + endif + + cart_func(i,j,k) = rad_func(ri) + end do + end do + end do + + end subroutine Spin1D + + + +! ======================================================================= + + + subroutine Spin1DMetric(grr,gtt,gxx,gxy, & + gxz,gyy,gyz,gzz,nx,ny,nz,x,y,z,r) + +! ................................................................ +! this subroutine converts spherical metric components to cartesian +! +! This code was totally ripped off of Francisco Guzmans spheretocart +! code. Used with his persmission. + + implicit NONE + integer nx,ny,nz + integer i,j,k + + CCTK_REAL x(nx,ny,nz),y(nx,ny,nz),z(nx,ny,nz),r(nx,ny,nz) + CCTK_REAL grr(nx,ny,nz),gtt(nx,ny,nz) + CCTK_REAL gxx(nx,ny,nz),gxy(nx,ny,nz),gxz(nx,ny,nz) + CCTK_REAL gyy(nx,ny,nz),gyz(nx,ny,nz),gzz(nx,ny,nz) + + CCTK_REAL xp,yp,zp,rp + +! define derivatives drx = (dr/dx) + CCTK_REAL drx,dry,drz + CCTK_REAL dtx,dty,dtz + CCTK_REAL dpx,dpy,dpz + +! this equation for gxx gyy gzz assumes that the metric is diagonal +! and that grr = 1 at the origin, but should compute for all points +! including the origin and the line x=y=0. + + + do k = 1, nz + do j = 1, ny + do i = 1, nx + + xp = x(i,j,k) + yp = y(i,j,k) + zp = z(i,j,k) + rp = r(i,j,k) + + if (rp.eq.0) then + + gxx(i,j,k) = 1. + gyy(i,j,k) = 1. + gzz(i,j,k) = 1. + gxy(i,j,k) = 0. + gxz(i,j,k) = 0. + gyz(i,j,k) = 0. + + else + + gxx(i,j,k) = (xp/rp)**2 * grr(i,j,k) + & + (zp**2 + yp**2)/(rp**2) + gyy(i,j,k) = (yp/rp)**2 * grr(i,j,k) + & + (xp**2 + zp**2)/(rp**2) + gzz(i,j,k) = (zp/rp)**2 * grr(i,j,k) + & + (xp**2 + yp**2)/(rp**2) + + if ((xp.eq.0) .and. (yp.eq.0)) then + + gxy(i,j,k) = 0. + gyz(i,j,k) = 0. + gxz(i,j,k) = 0. + + else + + gxy(i,j,k) = (xp * yp/(rp**2) * grr(i,j,k) & + + (xp*yp*zp**2)/(rp**2*(xp**2+yp**2)) & + - (yp*xp)/(xp**2+yp**2)) + + gxz(i,j,k) = (xp * zp/(rp**2) & + * grr(i,j,k) - (xp * zp/(rp**2))) + + gyz(i,j,k) = (yp * zp/(rp**2) & + * grr(i,j,k) - (yp * zp/(rp**2))) + + endif + + endif + + enddo + enddo + enddo +#else + write(6,*) 'GRHydro_TOV: Does nothing. Recompile with odepack' +#endif + + return + end + + + |