ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c MAHC_IVP_CONFORMAL ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Makes the conformal transformation calculated in thorn_IVP on c the hydro variables. c (for now, uses P = (gamma-1) rho eps) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_DefineThorn.h" subroutine mahc_ivp_conformal(CCTK_FARGUMENTS) implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS c LOCAL VARS integer i,j,k,newt_count,CCTK_Equals,reduce_handle,ierr CCTK_REAL rhoyork,jx_york,jy_york,jz_york,j2 CCTK_REAL rhonew, rhoold,f,df,mass_scale CCTK_REAL vlowx,vlowy,vlowz,det,rho_global_max CCTK_REAL rho_high,rho_low,rho_local_max,rhoh logical ivp_eos_tables #ifdef CACTUSEOS_EOS_1D_BASE #include "EOS_1D_Base.inc" #endif ivp_eos_tables = CCTK_EQUALS(IVP_eos,"tables") if(ivp_eos_tables) then #ifdef CACTUSEOS_EOS_1D_BASE c first, get rho_global_max so we know how to scale this c problem rho_local_max = maxval(rho) call CCTK_ReductionArrayHandle(reduce_handle,"maximum") call CCTK_ReduceLocalScalar(ierr ,cctkGH,-1,reduce_handle, . rho_local_max,rho_global_max,CCTK_VARIABLE_REAL) if (ierr.ne.0) then call CCTK_WARN(1,"Reduction of max failed in ivp!"); write(*,*)'ierr = ',ierr endif do k=1,cctk_lsh(3) do j=1,cctk_lsh(2) do i=1,cctk_lsh(1) c get physical york variables rhoyork = (ivp_psi(i,j,k)**(ivp_matterdens_pow)) * . ( (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)) jx_york = (ivp_psi(i,j,k)**(ivp_matterdens_pow - 2.0d0)) * . ( (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k))* . (w_lorentz(i,j,k)**2) * velx(i,j,k)) jy_york = (ivp_psi(i,j,k)**(ivp_matterdens_pow - 2.0d0)) * . ( (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k))* . (w_lorentz(i,j,k)**2) * vely(i,j,k)) jz_york = (ivp_psi(i,j,k)**(ivp_matterdens_pow - 2.0d0)) * . ( (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k))* . (w_lorentz(i,j,k)**2) * velz(i,j,k)) c get j2 (g is now the physical metric) j2 = gxx(i,j,k)*jx_york**2 + gyy(i,j,k)*jy_york**2 + . gzz(i,j,k)*jz_york**2 + 2.0d0*gxy(i,j,k)*jx_york*jy_york + . 2.0d0*gxz(i,j,k)*jx_york*jz_york + . 2.0d0*gyz(i,j,k)*jy_york*jz_york c physical check if(rhoyork .le. 0.0d0) then write(*,*) 'mahc_ivp_conformal: physical error: ' write(*,*) ' total energy density less than 0!!!' write(*,*) ' position: (x,y,z,r) = ',x(i,j,k),y(i,j,k), . z(i,j,k),r(i,j,k) write(*,*) ' rhoyork = ', rhoyork write(*,*) ' j2 = ',j2 STOP endif rhonew = rho(i,j,k) eps(i,j,k) = EOS_1D_epsfromrho(IVP_EOS_1D_handle,rhonew) press(i,j,k) = EOS_1D_pressfromrho(IVP_EOS_1D_handle,rhonew) rhoh = rhonew + rhonew*eps(i,j,k) + press(i,j,k) f = rhoyork - 0.5d0*(rhoh + sqrt(rhoh**2 + 4.0d0*j2))+press(i,j,k) f = f/rho_global_max c if f is above 0, then we have our low, get the high if(f .ge. 0.0d0) then rho_low = rhonew rho_high = rho_low do while (f .ge. 0.0d0) rho_high = rho_high * 1.1d0 eps(i,j,k) = EOS_1D_epsfromrho(IVP_EOS_1D_handle,rho_high) press(i,j,k) = EOS_1D_pressfromrho(IVP_EOS_1D_handle,rho_high) rhoh = rho_high + rho_high*eps(i,j,k) + press(i,j,k) f = rhoyork - 0.5d0*(rhoh + sqrt(rhoh**2 + 4.0d0*j2))+ . press(i,j,k) f = f/rho_global_max enddo endif c if f is below 0, then we have our high, get the low if(f .lt. 0.0d0) then rho_high = rhonew rho_low = rho_high do while (f .lt. 0.0d0) rho_low = rho_low / 1.1d0 eps(i,j,k) = EOS_1D_epsfromrho(IVP_EOS_1D_handle,rho_low) press(i,j,k) = EOS_1D_pressfromrho(IVP_EOS_1D_handle,rho_low) rhoh = rho_low + rho_low*eps(i,j,k) + press(i,j,k) f = rhoyork - 0.5d0*(rhoh + sqrt(rhoh**2 + 4.0d0*j2))+ . press(i,j,k) f = f/rho_global_max enddo endif c now, bisect away! rhonew = 0.5d0*(rho_low + rho_high) eps(i,j,k) = EOS_1D_epsfromrho(IVP_EOS_1D_handle,rhonew) press(i,j,k) = EOS_1D_pressfromrho(IVP_EOS_1D_handle,rhonew) rhoh = rhonew + rhonew*eps(i,j,k) + press(i,j,k) f = rhoyork - 0.5d0*(rhoh + sqrt(rhoh**2 + 4.0d0*j2))+press(i,j,k) f = f/rho_global_max do while ( (abs(f) .gt. 1.0d-14) .and. . (rho_high-rho_low .gt. 1.0d-11*rho_global_max)) rhonew = 0.5d0*(rho_low + rho_high) eps(i,j,k) = EOS_1D_epsfromrho(IVP_EOS_1D_handle,rhonew) press(i,j,k) = EOS_1D_pressfromrho(IVP_EOS_1D_handle,rhonew) rhoh = rhonew + rhonew*eps(i,j,k) + press(i,j,k) f = rhoyork - . 0.5d0*(rhoh + sqrt(rhoh**2 + 4.0d0*j2))+press(i,j,k) f = f/rho_global_max if(f .ge. 0.0d0) then rho_low = rhonew else rho_high = rhonew endif enddo c set all the primative variables! rho(i,j,k) = rhonew eps(i,j,k) = EOS_1D_epsfromrho(IVP_EOS_1D_handle,rhonew) press(i,j,k) = EOS_1D_pressfromrho(IVP_EOS_1D_handle,rhonew) velx(i,j,k) = jx_york / (rhoyork + press(i,j,k)) vely(i,j,k) = jy_york / (rhoyork + press(i,j,k)) velz(i,j,k) = jz_york / (rhoyork + press(i,j,k)) c set all of the evolved variables! 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)) 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 return #else write(*,*) 'if you want to use eos tables in ivp, you' write(*,*) 'need to compile with EOS_1D_Base' call CCTK_WARN(0,"mahc_ivp_conformal") #endif endif c regular assumption P=(gamma-1)*rho*eps = k*rho**gamma do k=1,cctk_lsh(3) do j=1,cctk_lsh(2) do i=1,cctk_lsh(1) c get physical york variables rhoyork = (ivp_psi(i,j,k)**(ivp_matterdens_pow)) * . ( (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)) jx_york = (ivp_psi(i,j,k)**(ivp_matterdens_pow - 2.0d0)) * . ( (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k))* . (w_lorentz(i,j,k)**2) * velx(i,j,k)) jy_york = (ivp_psi(i,j,k)**(ivp_matterdens_pow - 2.0d0)) * . ( (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k))* . (w_lorentz(i,j,k)**2) * vely(i,j,k)) jz_york = (ivp_psi(i,j,k)**(ivp_matterdens_pow - 2.0d0)) * . ( (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k))* . (w_lorentz(i,j,k)**2) * velz(i,j,k)) c get j2 (g is now the physical metric) j2 = gxx(i,j,k)*jx_york**2 + gyy(i,j,k)*jy_york**2 + . gzz(i,j,k)*jz_york**2 + 2.0d0*gxy(i,j,k)*jx_york*jy_york + . 2.0d0*gxz(i,j,k)*jx_york*jz_york + . 2.0d0*gyz(i,j,k)*jy_york*jz_york c physical check if(rhoyork .le. 0.0d0) then write(*,*) 'mahc_ivp_conformal: physical error: ' write(*,*) ' total energy density less than 0!!!' write(*,*) ' position: (x,y,z,r) = ',x(i,j,k),y(i,j,k), . z(i,j,k),r(i,j,k) write(*,*) ' rhoyork = ', rhoyork write(*,*) ' j2 = ',j2 STOP endif c now do Newtons method on density (this method assumes c P = (gamma-1) rho eps = k rho**gamma) rhoold = rho(i,j,k) rhonew = rhoold newt_count = 0 f = & -(eos_k*rhonew**eos_gamma)-rhoyork+5.d-1*(rhonew+eos_gamma*eos_k*rho & new**eos_gamma/(-1.d0+eos_gamma))*(1.d0+sqrt(1.d0+4.d0*j2/(rhonew+ & eos_gamma*eos_k*rhonew**eos_gamma/(-1.d0+eos_gamma))**2)) mass_scale = max(rhoyork,sqrt(j2)) do while (abs(f/mass_scale) .gt. mahc_ivp_conformal_tol) newt_count = newt_count + 1 rhoold = rhonew df = & -(eos_gamma*eos_k*rhonew**(-1.d0+eos_gamma))-2.d0*j2*(1.d0+eos_gamma & **2*eos_k*rhonew**(-1.d0+eos_gamma)/(-1.d0+eos_gamma))/((rhonew+eos_ & gamma*eos_k*rhonew**eos_gamma/(-1.d0+eos_gamma))**2*sqrt(1.d0+4.d0* & j2/(rhonew+eos_gamma*eos_k*rhonew**eos_gamma/(-1.d0+eos_gamma))**2)) & +5.d-1*(1.d0+eos_gamma**2*eos_k*rhonew**(-1.d0+eos_gamma)/(-1.d0+eo & s_gamma))*(1.d0+sqrt(1.d0+4.d0*j2/(rhonew+eos_gamma*eos_k*rhonew**e & os_gamma/(-1.d0+eos_gamma))**2)) rhonew = rhoold - f/df f = & -(eos_k*rhonew**eos_gamma)-rhoyork+5.d-1*(rhonew+eos_gamma*eos_k*rho & new**eos_gamma/(-1.d0+eos_gamma))*(1.d0+sqrt(1.d0+4.d0*j2/(rhonew+ & eos_gamma*eos_k*rhonew**eos_gamma/(-1.d0+eos_gamma))**2)) if(newt_count .gt. 50) then write(*,*) 'mahc_ivp_conformal: newtons method not converging' write(*,*) ' in 50 steps' write(*,*) ' position: (x,y,z,r) = ',x(i,j,k),y(i,j,k), . z(i,j,k),r(i,j,k) write(*,*) ' rhonew = ',rhonew write(*,*) ' rhoold = ',rhoold write(*,*) ' minimizing function f = ',f write(*,*) ' df/drho = ',df STOP endif enddo c now, do a couple more Newtons method sweeps just for fun! do newt_count=1,3 rhoold = rhonew df = & -(eos_gamma*eos_k*rhonew**(-1.d0+eos_gamma))-2.d0*j2*(1.d0+eos_gamma & **2*eos_k*rhonew**(-1.d0+eos_gamma)/(-1.d0+eos_gamma))/((rhonew+eos_ & gamma*eos_k*rhonew**eos_gamma/(-1.d0+eos_gamma))**2*sqrt(1.d0+4.d0* & j2/(rhonew+eos_gamma*eos_k*rhonew**eos_gamma/(-1.d0+eos_gamma))**2)) & +5.d-1*(1.d0+eos_gamma**2*eos_k*rhonew**(-1.d0+eos_gamma)/(-1.d0+eo & s_gamma))*(1.d0+sqrt(1.d0+4.d0*j2/(rhonew+eos_gamma*eos_k*rhonew**e & os_gamma/(-1.d0+eos_gamma))**2)) rhonew = rhoold - f/df f = & -(eos_k*rhonew**eos_gamma)-rhoyork+5.d-1*(rhonew+eos_gamma*eos_k*rho & new**eos_gamma/(-1.d0+eos_gamma))*(1.d0+sqrt(1.d0+4.d0*j2/(rhonew+ & eos_gamma*eos_k*rhonew**eos_gamma/(-1.d0+eos_gamma))**2)) enddo c set all the primative variables! rho(i,j,k) = rhonew eps(i,j,k) = eos_k * (rhonew**(eos_gamma-1.0d0))/(eos_gamma-1.0d0) press(i,j,k) = (eos_gamma-1.0d0) * rhonew * eps(i,j,k) velx(i,j,k) = jx_york / (rhoyork + press(i,j,k)) vely(i,j,k) = jy_york / (rhoyork + press(i,j,k)) velz(i,j,k) = jz_york / (rhoyork + press(i,j,k)) c set all of the evolved variables! 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)) 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 now set the evolved variables return end