c Mahc_Init_Data: Various Initial Data Routines for Mahc_Evolve c Copyright (C) 2001 Mark Miller /*@@ @file boost_multiple_star.F @date January 2000 @author Mark Miller @desc Multiple TOV stars boosted in arbitrary directions. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_DefineThorn.h" /*@@ @routine boost_multiple_star @date January 2000 @author Mark Miller @desc Multiple TOV stars boosted in arbitrary direction. @enddesc @calls @calledby @history @endhistory @@*/ subroutine boost_multiple_star(CCTK_FARGUMENTS) implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS c LOCAL VARS c TOV SOLVER PARAMETERS CCTK_REAL rmaxtov integer ntov integer i,j,k,imaxtov,nx,ny,nz,CCTK_Equals CCTK_REAL rr,xx,yy,zz,mval,phival,mgrr,dgrr,dx CCTK_REAL alpha,dalpha,pi CCTK_REAL one,zero,two CCTK_REAL star1x,star1y,star1z,star2x,star2y,star2z CCTK_REAL rhotemp,epstemp,gxxtemp,gxytemp,gxztemp,gyytemp,gyztemp CCTK_REAL gzztemp,presstemp CCTK_REAL betaxtemp,betaytemp,betaztemp CCTK_REAL hxxtemp,hxytemp,hxztemp,hyytemp,hyztemp,hzztemp CCTK_REAL alptemp,velxtemp,velytemp,velztemp CCTK_REAL relcoordx,relcoordy,relcoordz,star_mass,star_rad CCTK_REAL w1,w2,r1,r2,temp1,temp2,wtparam,tiny CCTK_REAL rho_temp,m_temp,phi_temp,eps_temp,press_temp CCTK_REAL rho_wt_1,rho_wt_2,vlowx,vlowy,vlowz,det,bms_rho_multifact 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') write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' write(*,*) '*/*/ MAHC */*/' write(*,*) '*/*/ MULTIPLE BOOSTED NEUTRON STARS */*/' write(*,*) '*/*/ ---------------------------------------- */*/' write(*,*) '*/*/ For those times when one is */*/' write(*,*) '*/*/ just not enough! */*/' write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) dx = cctk_delta_space(1) c make sure we only have 3 NS if((boost_multi_star_num .le. 0).or.(bms_centden_1.le.0.0d0)) then write(*,*) 'Error:boost_multiple_star: must have at least 1 NS!' endif if(boost_multi_star_num .gt. 3) then write(*,*) 'Warning:boost_multiple_star: MAHC only supports 3' write(*,*) ' NSs right now!' boost_multi_star_num = 3 endif one = 1.0d0 zero = 0.0d0 two = 2.0d0 tiny = 1.0d-30 pi = 4.0d0 * atan(1.0d0) write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' write(*,*) '*/*/ tov_eos_gamma is ',eos_gamma write(*,*) '*/*/ eos_k is ',eos_k write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' c allocate the 1D arrays rmaxtov = mahc_init_rmaxtov ntov = mahc_init_ntov c solve for the 1st NS call solvetov_iso(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,bms_centden_1,eos_gamma,eos_k,rmaxtov, . mahc_init_eos_handle) write(*,*) '*/*/ boost_multiple_star: star 1: ' write(*,*) '*/*/ central rest mass density:',bms_centden_1 write(*,*) '*/*/ star1 center: ',bms_xpos_1,bms_ypos_1,bms_zpos_1 write(*,*) '*/*/ star1 boost parameters: ', . bms_boost_x_1,bms_boost_y_1,bms_boost_z_1 write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' c do star 1 do k=1,nz do j=1,ny do i=1,nx relcoordx = x(i,j,k) - bms_xpos_1 relcoordy = y(i,j,k) - bms_ypos_1 relcoordz = z(i,j,k) - bms_zpos_1 r1 = sqrt(relcoordx**2 + relcoordy**2 + relcoordz**2) call boost_isons_pt(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . bms_boost_x_1,bms_boost_y_1,bms_boost_z_1,zero,relcoordx, . relcoordy,relcoordz,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(boost_multi_star_num .ge. 2) then if(bms_centden_2 .le. 0.0d0) then write(*,*) 'boost_multiple_star: hey, one of the stars has' write(*,*) ' zero central density!' STOP endif c solve for the 2nd NS call solvetov_iso(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,bms_centden_2,eos_gamma,eos_k,rmaxtov, . mahc_init_eos_handle) write(*,*) '*/*/ boost_multiple_star: star 2: ' write(*,*) '*/*/ central rest mass density:',bms_centden_2 write(*,*) '*/*/ star2 center: ',bms_xpos_2,bms_ypos_2,bms_zpos_2 write(*,*) '*/*/ star2 boost parameters: ', . bms_boost_x_2,bms_boost_y_2,bms_boost_z_2 write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' c for some reason, the T3E really hates passing around variables c that have not been assigned a value. Picky picky picky. rhotemp = 0.0d0 epstemp = 0.0d0 gxxtemp = 0.0d0 gxytemp = 0.0d0 gxztemp = 0.0d0 gyytemp = 0.0d0 gyztemp = 0.0d0 gzztemp = 0.0d0 alptemp = 0.0d0 betaxtemp = 0.0d0 betaytemp = 0.0d0 betaztemp = 0.0d0 velxtemp = 0.0d0 velytemp = 0.0d0 velztemp = 0.0d0 hxxtemp = 0.0d0 hxytemp = 0.0d0 hxztemp = 0.0d0 hyytemp = 0.0d0 hyztemp = 0.0d0 hzztemp = 0.0d0 c do star 2 do k=1,nz do j=1,ny do i=1,nx relcoordx = x(i,j,k) - bms_xpos_2 relcoordy = y(i,j,k) - bms_ypos_2 relcoordz = z(i,j,k) - bms_zpos_2 r1 = sqrt(relcoordx**2 + relcoordy**2 + relcoordz**2) call boost_isons_pt(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . bms_boost_x_2,bms_boost_y_2,bms_boost_z_2,zero,relcoordx, . relcoordy,relcoordz,eos_gamma,eos_k,presstemp,rhotemp,epstemp, . gxxtemp,gxytemp,gxztemp,gyytemp,gyztemp, . gzztemp,alptemp,betaxtemp,betaytemp, . betaztemp,velxtemp,velytemp,velztemp, . hxxtemp,hxytemp,hxztemp,hyytemp,hyztemp, . hzztemp,dx,mahc_init_eos_handle) rho_wt_1 = rho(i,j,k) / (rho(i,j,k) + rhotemp + tiny) rho_wt_2 = rhotemp / (rho(i,j,k) + rhotemp + tiny) c weight average these variables: velx(i,j,k) = rho_wt_1*velx(i,j,k) + rho_wt_2*velxtemp vely(i,j,k) = rho_wt_1*vely(i,j,k) + rho_wt_2*velytemp velz(i,j,k) = rho_wt_1*velz(i,j,k) + rho_wt_2*velztemp betax(i,j,k) =rho_wt_1*betax(i,j,k)+ rho_wt_2*betaxtemp betay(i,j,k) =rho_wt_1*betay(i,j,k)+ rho_wt_2*betaytemp betaz(i,j,k) =rho_wt_1*betaz(i,j,k)+ rho_wt_2*betaztemp c add these variables: rho(i,j,k) = rho(i,j,k) + rhotemp eps(i,j,k) = eps(i,j,k) + epstemp press(i,j,k) = press(i,j,k) + presstemp gxy(i,j,k) = gxy(i,j,k)+gxytemp gxz(i,j,k) = gxz(i,j,k)+gxztemp gyz(i,j,k) = gyz(i,j,k)+gyztemp kxx(i,j,k) = kxx(i,j,k)+hxxtemp kxy(i,j,k) = kxy(i,j,k)+hxytemp kxz(i,j,k) = kxz(i,j,k)+hxztemp kyy(i,j,k) = kyy(i,j,k)+hyytemp kyz(i,j,k) = kyz(i,j,k)+hyztemp kzz(i,j,k) = kzz(i,j,k)+hzztemp c add deviations from 1 for these variables: gxx(i,j,k) = gxx(i,j,k) + gxxtemp - 1.0d0 gyy(i,j,k) = gyy(i,j,k) + gyytemp - 1.0d0 gzz(i,j,k) = gzz(i,j,k) + gzztemp - 1.0d0 alp(i,j,k) = alp(i,j,k) + alptemp - 1.0d0 enddo enddo enddo endif if(boost_multi_star_num .ge. 3) then write(*,*) 'boost_multiple_star: bug mark to add 3 NS functionality.' write(*,*) 'Remind him that it is not too bad, things are' write(*,*) 'already set up for it!' call CCTK_WARN(0,"boost_multi_star: need to put in 3 NSs!") endif c 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 if(.not. 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_multiple_star") #endif endif c set up consistent evolved variables, though you are c going to use thorn_IVP on this initial data if you c want to satisfy the constraints! do k=1,nz do j=1,ny do i=1,nx 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_multiple_star") #endif endif 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) w_lorentz(i,j,k) = 1.0d0/sqrt(1.0d0 - . velx(i,j,k)*vlowx - vely(i,j,k)*vlowy -velz(i,j,k)*vlowz) 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) dens(i,j,k) = sqrt(det) * rho(i,j,k) * w_lorentz(i,j,k) tau(i,j,k) = sqrt(det) * ((w_lorentz(i,j,k)**2)* . (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k)) - . press(i,j,k)) - dens(i,j,k) sx(i,j,k) = sqrt(det) * (w_lorentz(i,j,k)**2)* . (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k))*vlowx sy(i,j,k) = sqrt(det) * (w_lorentz(i,j,k)**2)* . (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k))*vlowy sz(i,j,k) = sqrt(det) * (w_lorentz(i,j,k)**2)* . (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k))*vlowz enddo enddo enddo c set shift to zero here betax = 0.0d0 betay = 0.0d0 betaz = 0.0d0 return end