c Mahc_Init_Data: Various Initial Data Routines for Mahc_Evolve c Copyright (C) 2001 Mark Miller /*@@ @file boost_star.F @date January 2000 @author Mark Miller @desc Calculates initial data that corresponds to a boosted TOV solution. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_DefineThorn.h" /*@@ @routine boost_star @date January 2000 @author Mark Miller @desc Calculates initial data that corresponds to a boosted TOV solution. @enddesc @calls @calledby @history @endhistory @@*/ subroutine boost_star(CCTK_FARGUMENTS) implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS c LOCAL VARS c TOV SOLVER PARAMETERS CCTK_REAL rmaxtov,evol_eos_gamma,tiny,dx c I have checked that the hamiltonian constraint converges c correctly up to dx = .1 for drtov = 1/512 CCTK_REAL star_mass,star_rad CCTK_REAL zero,one integer i,j,k,imaxtov,CCTK_Equals,ntov,nx,ny,nz c metric derivatives CCTK_REAL dx_gxx,dx_gxy,dx_gxz,dx_gyy,dx_gyz,dx_gzz CCTK_REAL dy_gxx,dy_gxy,dy_gxz,dy_gyy,dy_gyz,dy_gzz CCTK_REAL dz_gxx,dz_gxy,dz_gxz,dz_gyy,dz_gyz,dz_gzz CCTK_REAL dx_betax,dx_betay,dx_betaz,dy_betax,dy_betay,dy_betaz CCTK_REAL dz_betax,dz_betay,dz_betaz CCTK_REAL liexx,liexy,liexz,lieyy,lieyz,liezz CCTK_REAL vlowx,vlowy,vlowz,det,my_time logical mahc_init_eos_tables #ifdef CACTUSEOS_EOS_1D_BASE #include "EOS_1D_Base.inc" #endif mahc_init_eos_tables = CCTK_EQUALS(mahc_init_eos,'tables') nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) dx = cctk_delta_space(1) write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' write(*,*) '*/*/ MAHC */*/' write(*,*) '*/*/ BOOSTED TOV INITIAL DATA */*/' write(*,*) '*/*/ ---------------- */*/' write(*,*) '*/*/ More fun than humans should */*/' write(*,*) '*/*/ be allowed to have!!! */*/' write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' zero = 0.0d0 one = 1.0d0 rmaxtov = mahc_init_rmaxtov ntov = mahc_init_ntov c we have got to have a shift if(CCTK_EQUALS(shift,"none")) then write(*,*) '*/*/ Hey, if you want a boosted TOV solution, you had' write(*,*) '*/*/ better turn on the shift. That is, unless you' write(*,*) '*/*/ really like solving some nasty integrability' write(*,*) '*/*/ conditions.....' write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' call CCTK_WARN(0,"boost_star: turn on shift!") endif c now, get the boost parameters write(*,*) '*/*/ boost parameters: ',boost_x,boost_y,boost_z write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' write(*,*) '*/*/ Solving TOV solution in the rest frame....' write(*,*) '*/*/ ...' c get the TOV solution in the unprimed (rest) frame call solvetov_iso(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,centden,eos_gamma,eos_k,rmaxtov, . mahc_init_eos_handle) write(*,*) '*/*/ Done with TOV solve.' write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' c Special initial data for boost_star_SL_test #ifdef THORN_ADM if(CONTAINS('boost_star_SL_test','yes')) then write(*,*) '*/*/ BOOST_STAR: special initial data for ' write(*,*) '*/*/ boost_star_SL_test ' write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' c level n-1 my_time = zero - dt do k=1,nz do j=1,ny do i=1,nx call boost_isons_pt(presstov,rhotov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,my_time,x(i,j,k),y(i,j,k),z(i,j,k), . eos_gamma,eos_k,rho(i,j,k),eps(i,j,k),ADM_gxx_p(i,j,k), . ADM_gxy_p(i,j,k),ADM_gxz_p(i,j,k),ADM_gyy_p(i,j,k), . ADM_gyz_p(i,j,k),ADM_gzz_p(i,j,k), . ADM_alp_p(i,j,k),betax(i,j,k),betay(i,j,k),betaz(i,j,k), . velx(i,j,k),vely(i,j,k),velz(i,j,k),kxx(i,j,k), . kxy(i,j,k),kxz(i,j,k),kyy(i,j,k),kyz(i,j,k),kzz(i,j,k),dx, . mahc_init_eos_handle) enddo enddo enddo c level n - 1/2 my_time = zero - dt/2.0d0 do k=1,nz do j=1,ny do i=1,nx call boost_isons_pt(presstov,rhotov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,my_time,x(i,j,k),y(i,j,k),z(i,j,k), . eos_gamma,eos_k,rho(i,j,k),eps(i,j,k),gxx(i,j,k), . gxy(i,j,k),gxz(i,j,k),gyy(i,j,k), . gyz(i,j,k),gzz(i,j,k), . alp(i,j,k),betax(i,j,k),betay(i,j,k),betaz(i,j,k), . velx(i,j,k),vely(i,j,k),velz(i,j,k),ADM_kxx_p(i,j,k), . ADM_kxy_p(i,j,k),ADM_kxz_p(i,j,k),ADM_kyy_p(i,j,k), . ADM_kyz_p(i,j,k),ADM_kzz_p(i,j,k),dx,mahc_init_eos_handle) enddo enddo enddo c level n + 1/2 my_time = zero + dt/2.0d0 do k=1,nz do j=1,ny do i=1,nx call boost_isons_pt(presstov,rhotov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,my_time,x(i,j,k),y(i,j,k),z(i,j,k), . eos_gamma,eos_k,rho(i,j,k),eps(i,j,k),gxx(i,j,k), . gxy(i,j,k),gxz(i,j,k),gyy(i,j,k), . gyz(i,j,k),gzz(i,j,k), . alp(i,j,k),betax(i,j,k),betay(i,j,k),betaz(i,j,k), . velx(i,j,k),vely(i,j,k),velz(i,j,k),ADM_kxx_stag(i,j,k), . ADM_kxy_stag(i,j,k),ADM_kxz_stag(i,j,k),ADM_kyy_stag(i,j,k), . ADM_kyz_stag(i,j,k),ADM_kzz_stag(i,j,k),dx,mahc_init_eos_handle) enddo enddo enddo endif #endif do k=1,nz do j=1,ny do i=1,nx call boost_isons_pt(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,zero,x(i,j,k),y(i,j,k),z(i,j,k), . eos_gamma,eos_k,press(i,j,k),rho(i,j,k),eps(i,j,k),gxx(i,j,k), . gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), . alp(i,j,k),betax(i,j,k),betay(i,j,k),betaz(i,j,k), . velx(i,j,k),vely(i,j,k),velz(i,j,k),kxx(i,j,k), . kxy(i,j,k),kxz(i,j,k),kyy(i,j,k),kyz(i,j,k),kzz(i,j,k),dx, . mahc_init_eos_handle) enddo enddo enddo if(zero_shift.eq.1) then betax = 0.0d0 betay = 0.0d0 betaz = 0.0d0 endif c now, since the extrinsic curvature only contains the qdot c terms, we must adjust by the Lie derivative wrt beta c do k=2,nz-1 c do j=2,ny-1 c do i=2,nx-1 c calculate the derivative terms needed. c dx_gxx = (gxx(i+1,j,k) - gxx(i-1,j,k))/(2.0D0*dx) c dx_gxy = (gxy(i+1,j,k) - gxy(i-1,j,k))/(2.0D0*dx) c dx_gxz = (gxz(i+1,j,k) - gxz(i-1,j,k))/(2.0D0*dx) c dx_gyy = (gyy(i+1,j,k) - gyy(i-1,j,k))/(2.0D0*dx) c dx_gyz = (gyz(i+1,j,k) - gyz(i-1,j,k))/(2.0D0*dx) c dx_gzz = (gzz(i+1,j,k) - gzz(i-1,j,k))/(2.0D0*dx) c dy_gxx = (gxx(i,j+1,k) - gxx(i,j-1,k))/(2.0D0*dy) c dy_gxy = (gxy(i,j+1,k) - gxy(i,j-1,k))/(2.0D0*dy) c dy_gxz = (gxz(i,j+1,k) - gxz(i,j-1,k))/(2.0D0*dy) c dy_gyy = (gyy(i,j+1,k) - gyy(i,j-1,k))/(2.0D0*dy) c dy_gyz = (gyz(i,j+1,k) - gyz(i,j-1,k))/(2.0D0*dy) c dy_gzz = (gzz(i,j+1,k) - gzz(i,j-1,k))/(2.0D0*dy) c dz_gxx = (gxx(i,j,k+1) - gxx(i,j,k-1))/(2.0D0*dz) c dz_gxy = (gxy(i,j,k+1) - gxy(i,j,k-1))/(2.0D0*dz) c dz_gxz = (gxz(i,j,k+1) - gxz(i,j,k-1))/(2.0D0*dz) c dz_gyy = (gyy(i,j,k+1) - gyy(i,j,k-1))/(2.0D0*dz) c dz_gyz = (gyz(i,j,k+1) - gyz(i,j,k-1))/(2.0D0*dz) c dz_gzz = (gzz(i,j,k+1) - gzz(i,j,k-1))/(2.0D0*dz) c dx_betax = (betax(i+1,j,k) - betax(i-1,j,k))/(2.0D0*dx) c dx_betay = (betay(i+1,j,k) - betay(i-1,j,k))/(2.0D0*dx) c dx_betaz = (betaz(i+1,j,k) - betaz(i-1,j,k))/(2.0D0*dx) c dy_betax = (betax(i,j+1,k) - betax(i,j-1,k))/(2.0D0*dy) c dy_betay = (betay(i,j+1,k) - betay(i,j-1,k))/(2.0D0*dy) c dy_betaz = (betaz(i,j+1,k) - betaz(i,j-1,k))/(2.0D0*dy) c dz_betax = (betax(i,j,k+1) - betax(i,j,k-1))/(2.0D0*dz) c dz_betay = (betay(i,j,k+1) - betay(i,j,k-1))/(2.0D0*dz) c dz_betaz = (betaz(i,j,k+1) - betaz(i,j,k-1))/(2.0D0*dz) c liexx = c . betax(i,j,k)*dx_gxx + betay(i,j,k)*dy_gxx + betaz(i,j,k)*dz_gxx + c . gxx(i,j,k)*dx_betax + gxy(i,j,k)*dx_betay + gxz(i,j,k)*dx_betaz + c . gxx(i,j,k)*dx_betax + gxy(i,j,k)*dx_betay + gxz(i,j,k)*dx_betaz c liexy = c . betax(i,j,k)*dx_gxy + betay(i,j,k)*dy_gxy + betaz(i,j,k)*dz_gxy + c . gxy(i,j,k)*dx_betax + gyy(i,j,k)*dx_betay + gyz(i,j,k)*dx_betaz + c . gxx(i,j,k)*dy_betax + gxy(i,j,k)*dy_betay + gxz(i,j,k)*dy_betaz c liexz = c . betax(i,j,k)*dx_gxz + betay(i,j,k)*dy_gxz + betaz(i,j,k)*dz_gxz + c . gxz(i,j,k)*dx_betax + gyz(i,j,k)*dx_betay + gzz(i,j,k)*dx_betaz + c . gxx(i,j,k)*dz_betax + gxy(i,j,k)*dz_betay + gxz(i,j,k)*dz_betaz c lieyy = c . betax(i,j,k)*dx_gyy + betay(i,j,k)*dy_gyy + betaz(i,j,k)*dz_gyy + c . gxy(i,j,k)*dy_betax + gyy(i,j,k)*dy_betay + gyz(i,j,k)*dy_betaz + c . gxy(i,j,k)*dy_betax + gyy(i,j,k)*dy_betay + gyz(i,j,k)*dy_betaz c lieyz = c . betax(i,j,k)*dx_gyz + betay(i,j,k)*dy_gyz + betaz(i,j,k)*dz_gyz + c . gxz(i,j,k)*dy_betax + gyz(i,j,k)*dy_betay + gzz(i,j,k)*dy_betaz + c . gxy(i,j,k)*dz_betax + gyy(i,j,k)*dz_betay + gyz(i,j,k)*dz_betaz c liezz = c . betax(i,j,k)*dx_gzz + betay(i,j,k)*dy_gzz + betaz(i,j,k)*dz_gzz + c . gxz(i,j,k)*dz_betax + gyz(i,j,k)*dz_betay + gzz(i,j,k)*dz_betaz + c . gxz(i,j,k)*dz_betax + gyz(i,j,k)*dz_betay + gzz(i,j,k)*dz_betaz c hxx(i,j,k) = hxx(i,j,k) + liexx/(2.0d0 * alp(i,j,k)) c hxy(i,j,k) = hxy(i,j,k) + liexy/(2.0d0 * alp(i,j,k)) c hxz(i,j,k) = hxz(i,j,k) + liexz/(2.0d0 * alp(i,j,k)) c hyy(i,j,k) = hyy(i,j,k) + lieyy/(2.0d0 * alp(i,j,k)) c hyz(i,j,k) = hyz(i,j,k) + lieyz/(2.0d0 * alp(i,j,k)) c hzz(i,j,k) = hzz(i,j,k) + liezz/(2.0d0 * alp(i,j,k)) c enddo c enddo c enddo c do not need to sync anymore c call syncallfuncs(GH) c now fix up the hydro: remember that we do not yet have c specific internal energy density.... c first, add atmosphere: do k=1,nz do j=1,ny do i=1,nx if(rho(i,j,k) .lt. centden/ns_atmos_fact) then rho(i,j,k) = centden/ns_atmos_fact endif enddo enddo enddo c fix up specific internal energy density if(.not.(CCTK_EQUALS(mahc_init_eos,'tables'))) then eps = (eos_k * rho**(eos_gamma - 1.0d0)) / . (eos_gamma - 1.0d0 + 1.e-30) else #ifdef CACTUSEOS_EOS_1D_BASE do k=1,nz do j=1,ny do i=1,nx eps(i,j,k) = EOS_1D_epsfromrho(mahc_init_eos_handle,rho(i,j,k)) enddo enddo enddo #else write(*,*) 'if you want to use tables in mahc_init_data, you' write(*,*) 'need to compile in with EOS_1D_Base' call CCTK_WARN(0,"boost_star") #endif endif c now, set evolved variables do k=1,nz do j=1,ny do i=1,nx vlowx = gxx(i,j,k)*velx(i,j,k) + gxy(i,j,k)*vely(i,j,k) + . gxz(i,j,k)*velz(i,j,k) vlowy = gxy(i,j,k)*velx(i,j,k) + gyy(i,j,k)*vely(i,j,k) + . gyz(i,j,k)*velz(i,j,k) vlowz = gxz(i,j,k)*velx(i,j,k) + gyz(i,j,k)*vely(i,j,k) + . gzz(i,j,k)*velz(i,j,k) det = ( gxx(i,j,k)*gyy(i,j,k)*gzz(i,j,k) + . 2.0d0*gxy(i,j,k)*gxz(i,j,k)*gyz(i,j,k) - . gxx(i,j,k) * (gyz(i,j,k)**2) - . gyy(i,j,k) * (gxz(i,j,k)**2) - . gzz(i,j,k) * (gxy(i,j,k)**2) ) w_lorentz(i,j,k) = 1.0d0 / sqrt(1.0d0 - vlowx*velx(i,j,k) - . vlowy*vely(i,j,k) - vlowz*velz(i,j,k)) c set the pressure if(.not. mahc_init_eos_tables) then press(i,j,k) = (eos_gamma - 1.0d0) * rho(i,j,k) * eps(i,j,k) else #ifdef CACTUSEOS_EOS_1D_BASE press(i,j,k) = EOS_1D_pressfromrho(mahc_init_eos_handle,rho(i,j,k)) #else write(*,*) 'if you want to use tables in mahc_init_data, you' write(*,*) 'need to compile in with EOS_1D_Base' call CCTK_WARN(0,"boost_star") #endif endif dens(i,j,k) = sqrt(det) * w_lorentz(i,j,k) * rho(i,j,k) tau(i,j,k) = sqrt(det) * ((rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + . press(i,j,k)) * w_lorentz(i,j,k)**2 - press(i,j,k)) . - dens(i,j,k) sx(i,j,k) = sqrt(det)* . (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k)) * . w_lorentz(i,j,k)**2 * vlowx sy(i,j,k) = sqrt(det)* . (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k)) * . w_lorentz(i,j,k)**2 * vlowy sz(i,j,k) = sqrt(det)* . (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k)) * . w_lorentz(i,j,k)**2 * vlowz enddo enddo enddo c this is taken care of in the scheduler now c call get_prim_varr(CCTK_FARGUMENTS) return end