[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 451)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Mon Jan 14 08:23:23 CST 2013
User: rhaas
Date: 2013/01/14 08:23 AM
Modified:
/trunk/
param.ccl
/trunk/src/
GRHydro_Con2PrimM.F90
Log:
GRHydro: call polytrope EOS in a region dominated by mag. field.
From: Bruno Coutinho Mundim <bcmsma at astro.rit.edu>
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Con2PrimM.F90
Delta lines: +60 -12
===================================================================
--- trunk/src/GRHydro_Con2PrimM.F90 2013-01-14 14:23:21 UTC (rev 450)
+++ trunk/src/GRHydro_Con2PrimM.F90 2013-01-14 14:23:23 UTC (rev 451)
@@ -93,6 +93,8 @@
CCTK_REAL :: velx_tmp, vely_tmp, velz_tmp, w_lorentz_tmp
CCTK_REAL :: Bvecx_tmp, Bvecy_tmp, Bvecz_tmp
+ CCTK_REAL :: bdotv, magpress
+
! Assume 3-metric is positive definite. Check deep inside the horizon
! if this is actually satisfied and if it is not then cast the metric
!as conformally flat only for con2prim inversion purposes.
@@ -166,7 +168,7 @@
!$OMP sdet,d2,s2,oob,bscon,bxhat,byhat,bzhat, &
!$OMP bhatscon,Wm,Wm0,Wm_resid,Wmold,s2m,s2m0,s2m_resid,s2mold,s2max, &
!$OMP taum,niter,rho_tmp,press_tmp,eps_tmp,velx_tmp,vely_tmp,velz_tmp, &
- !$OMP w_lorentz_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp)
+ !$OMP w_lorentz_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,bdotv,magpress)
do k = 1, nz
do j = 1, ny
do i = 1, nx
@@ -345,6 +347,18 @@
epsnegative,GRHydro_C2P_failed(i,j,k))
if(GRHydro_C2P_failed(i,j,k).ne.0) then
+
+ xrho=1.0d0; xtemp=0.0d0; xeps=1.0d0
+ call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ xrho,xeps,xtemp,xye,xpress,keyerr,anyerr)
+ local_K = xpress(1);
+
+ xrho=10.0d0; xeps=1.0d0
+ call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ xrho,xeps,xtemp,xye,xpress,keyerr,anyerr)
+ local_pgam=log(xpress(1)/local_K)/log(xrho(1))
+ sc = local_K*dens(i,j,k)
+
if(sdet.ge.sqrtdet_thr) then
GRHydro_C2P_failed(i,j,k) = 0
@@ -358,18 +372,7 @@
Bvecx_tmp = Bprim(i,j,k,1)
Bvecy_tmp = Bprim(i,j,k,2)
Bvecz_tmp = Bprim(i,j,k,3)
-
- xrho=1.0d0; xtemp=0.0d0; xeps=1.0d0
- call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
- xrho,xeps,xtemp,xye,xpress,keyerr,anyerr)
- local_K = xpress(1);
- xrho=10.0d0; xeps=1.0d0
- call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
- xrho,xeps,xtemp,xye,xpress,keyerr,anyerr)
- local_pgam=log(xpress(1)/local_K)/log(xrho(1))
- sc = local_K*dens(i,j,k)
-
call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, &
dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, &
Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,&
@@ -391,6 +394,51 @@
cycle
end if
end if
+
+ bdotv=g11(i,j,k)*Bprim(i,j,k,1)*vup(i,j,k,1)+ &
+ g22(i,j,k)*Bprim(i,j,k,2)*vup(i,j,k,2)+ &
+ g33(i,j,k)*Bprim(i,j,k,3)*vup(i,j,k,3)+ &
+ 2.0*(g12(i,j,k)*Bprim(i,j,k,1)*vup(i,j,k,2)+ &
+ g13(i,j,k)*Bprim(i,j,k,1)*vup(i,j,k,3)+ &
+ g23(i,j,k)*Bprim(i,j,k,2)*vup(i,j,k,3))
+
+ magpress = 0.5d0*(b2/w_lorentz(i,j,k)**2+bdotv**2)
+
+ if(rho(i,j,k)*eps(i,j,k)*max_magnetic_to_gas_pressure_ratio.le.magpress) then
+ GRHydro_C2P_failed(i,j,k) = 0
+
+ rho_tmp = rho(i,j,k)
+ press_tmp = press(i,j,k)
+ eps_tmp = eps(i,j,k)
+ velx_tmp = vup(i,j,k,1)
+ vely_tmp = vup(i,j,k,2)
+ velz_tmp = vup(i,j,k,3)
+ w_lorentz_tmp = w_lorentz(i,j,k)
+ Bvecx_tmp = Bprim(i,j,k,1)
+ Bvecy_tmp = Bprim(i,j,k,2)
+ Bvecz_tmp = Bprim(i,j,k,3)
+
+ call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, &
+ dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, &
+ Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,&
+ velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,&
+ Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,&
+ g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), &
+ uxx,uxy,uxz,uyy,uyz,uzz,det, &
+ epsnegative,GRHydro_C2P_failed(i,j,k))
+
+ rho(i,j,k) = rho_tmp
+ press(i,j,k) = press_tmp
+ eps(i,j,k) = eps_tmp
+ vup(i,j,k,1) = velx_tmp
+ vup(i,j,k,2) = vely_tmp
+ vup(i,j,k,3) = velz_tmp
+ w_lorentz(i,j,k) = w_lorentz_tmp
+ Bprim(i,j,k,1) = Bvecx_tmp
+ Bprim(i,j,k,2) = Bvecy_tmp
+ Bprim(i,j,k,3) = Bvecz_tmp
+ cycle
+ end if
end if
else ! if(evolve_temper.eq.0) then
Directory: /trunk/
==================
File [modified]: param.ccl
Delta lines: +5 -0
===================================================================
--- trunk/param.ccl 2013-01-14 14:23:21 UTC (rev 450)
+++ trunk/param.ccl 2013-01-14 14:23:23 UTC (rev 451)
@@ -356,6 +356,11 @@
1.0: :: "Larger values guarantees this sort of rescaling only deep inside the horizon"
} 30.0
+REAL max_magnetic_to_gas_pressure_ratio "consider pressure to be magnetically dominated if magnetic pressure to gas pressure ratio is higher than this" STEERABLE=ALWAYS
+{
+ (0:* :: "any positive value, eg. 100."
+ -1.0 :: "disable"
+} -1.0
boolean c2p_reset_pressure "If the pressure guess is unphysical should we recompute it?"
{
More information about the Commits
mailing list