User: cott
Date: 2010/07/02 04:57 PM
Modified:
/trunk/src/
GRHydro_Con2Prim.F90, GRHydro_Eigenproblem.F90, GRHydro_Eigenproblem_Marquina.F90, GRHydro_EoSChangeGamma.F90, GRHydro_Marquina.F90, GRHydro_Prim2Con.F90
Log:
* Changes to incorporate the EOS_Omni interace.
GRHydro has now 3 EOS interfaces. Rejoice!
To get the new one compiled in, use -DUSE_EOS_OMNI in
CPPFLAGS and FPPFLAGS.
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Con2Prim.F90
Delta lines: +233 -16
===================================================================
--- trunk/src/GRHydro_Con2Prim.F90 2010-07-02 21:51:08 UTC (rev 131)
+++ trunk/src/GRHydro_Con2Prim.F90 2010-07-02 21:57:07 UTC (rev 132)
@@ -150,6 +150,21 @@
CCTK_REAL :: local_min_tracer
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: poly_eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ poly_eoskey = GRHydro_poly_eoskey
+! end EOS Omni vars
+#endif
+
call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
type2_bits = -1
@@ -163,9 +178,17 @@
else
local_min_tracer = 0d0
end if
-
+
+#if USE_EOS_OMNI
+ ! this is a poly call
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ GRHydro_rho_min,xeps,xtemp,xye,pmin,keyerr,anyerr)
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
+#else
pmin = EOS_Pressure(GRHydro_polytrope_handle, GRHydro_rho_min, 1.0d0)
epsmin = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, GRHydro_rho_min, pmin)
+#endif
!$OMP PARALLEL DO PRIVATE(i,j,itracer,&
!$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, psi4pt, epsnegative)
@@ -221,8 +244,16 @@
scon(i,j,k,:) = 0.d0
vel(i,j,k,:) = 0.d0
w_lorentz(i,j,k) = 1.d0
+#if USE_EOS_OMNI
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
+#else
press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k))
eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k))
+#endif
! w_lorentz=1, so the expression for tau reduces to:
tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k)
@@ -235,9 +266,9 @@
scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vel(i,j,k,1),vel(i,j,k,2), &
vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
- z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
+ z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, &
+ GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
-
if (epsnegative) then
#if 1
@@ -267,8 +298,16 @@
scon(i,j,k,:) = 0.d0
vel(i,j,k,:) = 0.d0
w_lorentz(i,j,k) = 1.d0
- press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k))
- eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k))
+#if USE_EOS_OMNI
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
+#else
+ press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k))
+ eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k))
+#endif
! w_lorentz=1, so the expression for tau reduces to:
tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k)
@@ -315,12 +354,12 @@
subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
velz, epsilon, press, w_lorentz, uxx, uxy, uxz, uyy, &
- uyz, uzz, det, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, GRHydro_reflevel, GRHydro_C2P_failed)
+ uyz, uzz, det, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, &
+ GRHydro_reflevel, GRHydro_C2P_failed)
implicit none
DECLARE_CCTK_PARAMETERS
-
CCTK_REAL dens, sx, sy, sz, tau, rho, velx, vely, velz, epsilon, &
press, uxx, uxy, uxz, uyy, uyz, uzz, det, w_lorentz, x, &
y, z, r, GRHydro_rho_min
@@ -336,6 +375,22 @@
#undef _EOS_BASE_INC_
#endif
#include "EOS_Base.inc"
+
+
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ eoskey = GRHydro_eoskey
+! end EOS Omni vars
+#endif
+
!!$ Undensitize the variables
@@ -348,9 +403,15 @@
2.*usx*usz*uxz + 2.*usy*usz*uyz
!!$ Set initial guess for pressure:
-
+#if USE_EOS_OMNI
+ call EOS_Omni_press(eoskey,keytemp,n,&
+ rho,epsilon,xtemp,xye,xpress,keyerr,anyerr)
+ pold = max(1.d-10,xpress)
+#else
pold = max(1.d-10,EOS_Pressure(handle, rho, epsilon))
+#endif
+
!!$ Check that the variables have a chance of being physical
if( (utau + pold + udens)**2 - s2 .le. 0.0d0) then
@@ -386,8 +447,16 @@
!!$ Calculate the function
+#if USE_EOS_OMNI
+ call EOS_Omni_press(eoskey,keytemp,n,&
+ rho,epsilon,xtemp,xye,xpress,keyerr,anyerr)
+
+ f = pold - xpress
+#else
f = pold - EOS_Pressure(handle, rho, epsilon)
+#endif
+
!!$Find the root
count = 0
@@ -430,8 +499,18 @@
goto 51
end if
+#if USE_EOS_OMNI
+ call EOS_Omni_DPressByDRho(eoskey,keytemp,n,&
+ rho,epsilon,xtemp,xye,dpressbydrho,keyerr,anyerr)
+
+ call EOS_Omni_DPressByDEps(eoskey,keytemp,n,&
+ rho,epsilon,xtemp,xye,dpressbydeps,keyerr,anyerr)
+#else
dpressbydrho = EOS_DPressByDRho(handle,rho,epsilon)
dpressbydeps = EOS_DPressByDEps(handle,rho,epsilon)
+#endif
+
+
temp1 = (utau+udens+pnew)**2 - s2
drhobydpress = udens * s2 / (sqrt(temp1)*(udens+utau+pnew)**2)
depsbydpress = pnew * s2 / (rho * (udens + utau + pnew) * temp1)
@@ -449,7 +528,14 @@
epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - &
udens)/udens
+#if USE_EOS_OMNI
+ call EOS_Omni_press(eoskey,keytemp,n,&
+ rho,epsilon,xtemp,xye,xpress,keyerr,anyerr)
+
+ f = pnew - xpress
+#else
f = pnew - EOS_Pressure(handle, rho, epsilon)
+#endif
!check whether rho is NaN
if (rho .ne. rho) then
@@ -466,8 +552,17 @@
do i=1,GRHydro_polish
+#if USE_EOS_OMNI
+ call EOS_Omni_DPressByDRho(eoskey,keytemp,n,&
+ rho,epsilon,xtemp,xye,dpressbydrho,keyerr,anyerr)
+
+ call EOS_Omni_DPressByDEps(eoskey,keytemp,n,&
+ rho,epsilon,xtemp,xye,dpressbydeps,keyerr,anyerr)
+#else
dpressbydrho = EOS_DPressByDRho(handle,rho,epsilon)
dpressbydeps = EOS_DPressByDEps(handle,rho,epsilon)
+#endif
+
temp1 = (utau+udens+pnew)**2 - s2
drhobydpress = udens * s2 / (sqrt(temp1)*(udens+utau+pnew)**2)
depsbydpress = pnew * s2 / (rho * (udens + utau + pnew) * temp1)
@@ -484,8 +579,16 @@
epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - &
udens)/udens
+#if USE_EOS_OMNI
+ call EOS_Omni_press(eoskey,keytemp,n, &
+ rho,epsilon,xtemp,xye,xpress,keyerr,anyerr)
+
+ f = pold - xpress
+#else
f = pold - EOS_Pressure(handle, rho, epsilon)
+#endif
+
enddo
!!$ BEGIN: Check for NaN value (2nd check)
@@ -575,14 +678,42 @@
CCTK_INT :: type2_bits
CCTK_REAL :: local_min_tracer
-
+
#ifdef _EOS_BASE_INC_
#undef _EOS_BASE_INC_
#endif
#include "EOS_Base.inc"
+
+#if USE_EOS_OMNI
+! begin EOS omni
+ CCTK_INT :: eoskey = 0
+ CCTK_INT :: poly_eoskey = 0
+ CCTK_INT :: keyerr(1) = 0
+ CCTK_INT :: anyerr = 0
+ CCTK_INT :: keytemp = 0
+ CCTK_INT :: n = 1
+ CCTK_REAL :: xye = 0.0d0
+ CCTK_REAL :: xeps = 0.0d0
+ CCTK_REAL :: xtemp = 0.0d0
+ eoskey = GRHydro_eoskey
+ poly_eoskey = GRHydro_poly_eoskey
+! end EOS omni
+#endif
+
+
+#if USE_EOS_OMNI
+ ! this is a poly call
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ GRHydro_rho_min,epsmin,xtemp,xye,pmin,epsmin,keyerr,anyerr)
+#else
pmin=EOS_Pressure(GRHydro_polytrope_handle, GRHydro_rho_min, 1.0d0)
epsmin = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, GRHydro_rho_min, pmin)
+#endif
+
call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
@@ -807,6 +938,7 @@
local_min_tracer = 0d0
end if
+
!!$ do k = GRHydro_stencil + 1, nz - GRHydro_stencil
!!$ do j = GRHydro_stencil + 1, ny - GRHydro_stencil
!!$ do i = GRHydro_stencil + 1, nx - GRHydro_stencil
@@ -907,6 +1039,22 @@
#endif
#include "EOS_Base.inc"
+
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: poly_eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ poly_eoskey = GRHydro_poly_eoskey
+! end EOS Omni vars
+#endif
+
!!$ Undensitize the variables
sqrtdet = sqrt(det)
@@ -919,10 +1067,21 @@
2.d0*usx*usz*uxz + 2.d0*usy*usz*uyz
!!$ Set initial guess for rho:
+ rhoold = max(GRHydro_rho_min,rho)
- rhoold = max(GRHydro_rho_min,rho)
+#if USE_EOS_OMNI
+! call EOS_Omni_press(poly_eoskey,keytemp,n,&
+! rhoold,xeps,xtemp,xye,xpress,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ rhoold,xeps,xtemp,xye,press,xeps,keyerr,anyerr)
+
+ enthalpy = 1.0d0 + xeps + press / rhoold
+#else
enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhoold, press) + &
EOS_Pressure(handle, rhoold, epsilon) / rhoold
+#endif
+
w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 ))
!!$ BEGIN: Check for NaN value (1st check)
@@ -934,8 +1093,12 @@
endif
!!$ END: Check for NaN value (1st check)
+#if USE_EOS_OMNI
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ rhoold,epsilon,xtemp,xye,press,keyerr,anyerr)
+#else
press = EOS_Pressure(handle, rhoold, epsilon)
-
+#endif
!!$ Calculate the function
f = rhoold * w_lorentz - udens
@@ -968,8 +1131,16 @@
rhonew = GRHydro_rho_min
udens = rhonew
dens = sqrtdet * rhonew
+#if USE_EOS_OMNI
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ rhonew,1.0d0,xtemp,xye,press,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
+#else
press = EOS_Pressure(handle, rhonew, 1.d0)
epsilon = EOS_SpecificIntEnergy(handle, rhonew, press)
+#endif
sx = 0.d0
sy = 0.d0
sz = 0.d0
@@ -985,7 +1156,13 @@
!!$ Ian has the feeling that this is an error in general. But for the
!!$ 2D_Polytrope it gives the right answers.
+#if USE_EOS_OMNI
+ call EOS_Omni_DPressByDrho(poly_eoskey,keytemp,n,&
+ rhonew,epsilon,xtemp,xye,denthalpy,keyerr,anyerr)
+ denthalpy = denthalpy / rhonew
+#else
denthalpy = EOS_DPressByDrho(handle, rhonew, epsilon) / rhonew
+#endif
df = w_lorentz - rhonew * s2 * denthalpy / &
(w_lorentz*(udens**2)*(enthalpy**3))
@@ -1005,12 +1182,22 @@
end if
!!$ Recalculate primitive variables and function
-
- enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhonew, press) + &
+
+#if USE_EOS_OMNI
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ rhonew,epsilon,xtemp,xye,press,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
+
+ enthalpy = 1.0d0 + epsilon + press / rhonew
+#else
+ enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhonew, press) + &
EOS_Pressure(handle, rhonew, epsilon) / rhonew
+ press = EOS_Pressure(handle, rhonew, epsilon)
+#endif
+
w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 ))
-
- press = EOS_Pressure(handle, rhonew, epsilon)
tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens
f = rhonew * w_lorentz - udens
@@ -1022,10 +1209,21 @@
if (rhonew .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then
rhonew = GRHydro_rho_min
udens = rhonew
- ! before 2009/10/07 dens was not reset and was used with its negative value below; this has since been changed, but the change produces significant changes in the results
+ ! before 2009/10/07 dens was not reset and was used with its negative
+ ! value below; this has since been changed, but the change produces
+ ! significant changes in the results
dens = sqrtdet * rhonew
+
+#if USE_EOS_OMNI
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ rhonew,1.0d0,xtemp,xye,press,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
+#else
press = EOS_Pressure(handle, rhonew, 1.d0)
epsilon = EOS_SpecificIntEnergy(handle, rhonew, press)
+#endif
! w_lorentz=1, so the expression for utau reduces to:
tau = rhonew + rhonew*epsilon - udens
sx = 0.d0
@@ -1040,9 +1238,19 @@
50 rho = rhonew
+#if USE_EOS_OMNI
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ rhonew,epsilon,xtemp,xye,press,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
+
+ enthalpy = 1.0d0 + epsilon + press / rhonew
+#else
enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhonew, press) + &
EOS_Pressure(handle, rhonew, epsilon) / rhonew
w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 ))
+#endif
!!$ BEGIN: Check for NaN value (2nd check)
if (w_lorentz .ne. w_lorentz) then
@@ -1053,8 +1261,16 @@
endif
!!$ END: Check for NaN value (2nd check)
+#if USE_EOS_OMNI
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ rhonew,epsilon,xtemp,xye,press,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
+#else
press = EOS_Pressure(handle, rhonew, epsilon)
epsilon = EOS_SpecificIntEnergy(handle, rhonew, press)
+#endif
tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens
if (epsilon .lt. 0.0d0) then
@@ -1116,6 +1332,7 @@
@@*/
+
subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS)
implicit none
File [modified]: GRHydro_Eigenproblem.F90
Delta lines: +92 -0
===================================================================
--- trunk/src/GRHydro_Eigenproblem.F90 2010-07-02 21:51:08 UTC (rev 131)
+++ trunk/src/GRHydro_Eigenproblem.F90 2010-07-02 21:57:07 UTC (rev 132)
@@ -63,16 +63,49 @@
#endif
#include "EOS_Base.inc"
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ eoskey = GRHydro_eoskey
+! end EOS Omni vars
+#endif
+
one = 1.0d0
two = 2.0d0
!!$ Set required fluid quantities
+
+
+#if USE_EOS_OMNI
+! call EOS_Omni_cs2(eoskey,keytemp,n,&
+! rho,eps,xtemp,xye,cs2,keyerr,anyerr)
+ call EOS_Omni_press(eoskey,keytemp,n,&
+ rho,eps,xtemp,xye,press,keyerr,anyerr)
+
+ call EOS_Omni_DPressByDEps(eoskey,keytemp,n,&
+ rho,eps,xtemp,xye,dpdeps,keyerr,anyerr)
+
+ call EOS_Omni_DPressByDRho(eoskey,keytemp,n,&
+ rho,eps,xtemp,xye,dpdrho,keyerr,anyerr)
+
+ cs2 = (dpdrho + press * dpdeps / (rho**2))/ &
+ (1.0d0 + eps + press/rho)
+#else
press = EOS_Pressure(handle,rho,eps)
dpdrho = EOS_DPressByDRho(handle,rho,eps)
dpdeps = EOS_DPressByDEps(handle,rho,eps)
cs2 = (dpdrho + press * dpdeps / (rho**2))/ &
(1.0d0 + eps + press/rho)
+#endif
vlowx = gxx*velx + gxy*vely + gxz*velz
vlowy = gxy*velx + gyy*vely + gyz*velz
@@ -100,6 +133,7 @@
lam(4) = lam3
lam(5) = lamp
+
end subroutine eigenvalues
/*@@
@@ -157,16 +191,47 @@
#endif
#include "EOS_Base.inc"
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ eoskey = GRHydro_eoskey
+! end EOS Omni vars
+#endif
+
one = 1.0d0
two = 2.0d0
!!$ Set required fluid quantities
+#if USE_EOS_OMNI
+ call EOS_Omni_press(eoskey,keytemp,n,&
+ rho,eps,xtemp,xye,press,keyerr,anyerr)
+
+ call EOS_Omni_DPressByDEps(eoskey,keytemp,n,&
+ rho,eps,xtemp,xye,dpdeps,keyerr,anyerr)
+
+ call EOS_Omni_DPressByDRho(eoskey,keytemp,n,&
+ rho,eps,xtemp,xye,dpdrho,keyerr,anyerr)
+
+ cs2 = (dpdrho + press * dpdeps / (rho**2))/ &
+ (1.0d0 + eps + press/rho)
+! call EOS_Omni_cs2(eoskey,keytemp,n,&
+! rho,eps,xtemp,xye,cs2,keyerr,anyerr)
+#else
press = EOS_Pressure(handle,rho,eps)
dpdrho = EOS_DPressByDRho(handle,rho,eps)
dpdeps = EOS_DPressByDEps(handle,rho,eps)
cs2 = (dpdrho + press * dpdeps / (rho**2))/ &
(1.0d0 + eps + press/rho)
+#endif
! if (cs2<0) cs2=0 ! this does not modify the roe crashing problem with shocktube
enthalpy = one + eps + press / rho
@@ -345,6 +410,7 @@
rflux(4) = abs(lam1) * du(4) + enthalpy * w * (vlowz * vxa)
rflux(5) = abs(lam1) * du(5) + enthalpy * w * (velx * vxb + vxa) - vxa
+
else
!!$Form Jacobian matrix in characteristic form from right eigenvectors.
@@ -471,6 +537,10 @@
roeflux4 = rflux(4)
roeflux5 = rflux(5)
+ if(roeflux1.ne.roeflux1) then
+ stop "bad herberd"
+ endif
+
end subroutine eigenproblem
/*@@
@@ -522,16 +592,38 @@
#endif
#include "EOS_Base.inc"
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ eoskey = GRHydro_eoskey
+! end EOS Omni vars
+#endif
+
one = 1.0d0
two = 2.0d0
!!$ Set required fluid quantities
+#if USE_EOS_OMNI
+ call EOS_Omni_press(eoskey,keytemp,n,&
+ rho,eps,xtemp,xye,press,keyerr,anyerr)
+ call EOS_Omni_cs2(eoskey,keytemp,n,&
+ rho,eps,xtemp,xye,cs2,keyerr,anyerr)
+#else
press = EOS_Pressure(handle,rho,eps)
dpdrho = EOS_DPressByDRho(handle,rho,eps)
dpdeps = EOS_DPressByDEps(handle,rho,eps)
cs2 = (dpdrho + press * dpdeps / (rho**2))/ &
(1.0d0 + eps + press/rho)
+#endif
! if (cs2<0) cs2=0 ! this does not modify the roe crashing problem with shocktube
enthalpy = one + eps + press / rho
File [modified]: GRHydro_Eigenproblem_Marquina.F90
Delta lines: +44 -5
===================================================================
--- trunk/src/GRHydro_Eigenproblem_Marquina.F90 2010-07-02 21:51:08 UTC (rev 131)
+++ trunk/src/GRHydro_Eigenproblem_Marquina.F90 2010-07-02 21:57:07 UTC (rev 132)
@@ -99,19 +99,46 @@
#endif
#include "EOS_Base.inc"
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ eoskey = GRHydro_eoskey
+! end EOS Omni vars
+#endif
+
one = 1.0d0
two = 2.0d0
!!$ LEFT
!!$ Set required fluid quantities
-
+ invrhol = one / rhol
+#if USE_EOS_OMNI
+ call EOS_Omni_press(eoskey,keytemp,n,&
+ rhol,epsl,xtemp,xye,pressl,keyerr,anyerr)
+! call EOS_Omni_cs2(eoskey,keytemp,n,&
+! rhol,epsl,xtemp,xye,cs2l,keyerr,anyerr)
+ call EOS_Omni_DPressByDEps(eoskey,keytemp,n,&
+ rhol,epsl,xtemp,xye,dpdepsl,keyerr,anyerr)
+ call EOS_Omni_DPressByDRho(eoskey,keytemp,n,&
+ rhol,epsl,xtemp,xye,dpdrhol,keyerr,anyerr)
+ cs2l = (dpdrhol + pressl * dpdepsl * invrhol**2)/ &
+ (1.0d0 + epsl + pressl*invrhol)
+#else
pressl = EOS_Pressure(handle,rhol,epsl)
dpdrhol = EOS_DPressByDRho(handle,rhol,epsl)
dpdepsl = EOS_DPressByDEps(handle,rhol,epsl)
- invrhol = one / rhol
cs2l = (dpdrhol + pressl * dpdepsl * invrhol**2)/ &
(1.0d0 + epsl + pressl*invrhol)
+#endif
enthalpyl = one + epsl + pressl * invrhol
@@ -161,14 +188,25 @@
!!$ RIGHT
!!$ Set required fluid quantities
-
+ invrhor = one / rhor
+#if USE_EOS_OMNI
+ call EOS_Omni_press(eoskey,keytemp,n,&
+ rhor,epsr,xtemp,xye,pressr,keyerr,anyerr)
+! call EOS_Omni_cs2(eoskey,keytemp,n,&
+! rhor,epsr,xtemp,xye,cs2r,keyerr,anyerr)
+ call EOS_Omni_DPressByDEps(eoskey,keytemp,n,&
+ rhor,epsr,xtemp,xye,dpdepsr,keyerr,anyerr)
+ call EOS_Omni_DPressByDRho(eoskey,keytemp,n,&
+ rhor,epsr,xtemp,xye,dpdrhor,keyerr,anyerr)
+ cs2r = (dpdrhor + pressr * dpdepsr * invrhor**2)/ &
+ (1.0d0 + epsr + pressr*invrhor)
+#else
pressr = EOS_Pressure(handle,rhor,epsr)
dpdrhor = EOS_DPressByDRho(handle,rhor,epsr)
dpdepsr = EOS_DPressByDEps(handle,rhor,epsr)
- invrhor = one / rhor
cs2r = (dpdrhor + pressr * dpdepsr * invrhor**2)/ &
(1.0d0 + epsr + pressr*invrhor)
-
+#endif
enthalpyr = one + epsr + pressr * invrhor
vlowxr = gxx*velxr + gxy*velyr + gxz*velzr
@@ -247,6 +285,7 @@
kappal = dpdepsl / (dpdepsl - rhol * cs2l)
+
!!$ Right eigenvector # 1
reivec1l(1) = kappal / (enthalpyl * wl)
File [modified]: GRHydro_EoSChangeGamma.F90
Delta lines: +90 -3
===================================================================
--- trunk/src/GRHydro_EoSChangeGamma.F90 2010-07-02 21:51:08 UTC (rev 131)
+++ trunk/src/GRHydro_EoSChangeGamma.F90 2010-07-02 21:57:07 UTC (rev 132)
@@ -38,7 +38,11 @@
subroutine GRHydro_EoSChangeGamma(CCTK_ARGUMENTS)
+#if USE_EOS_OMNI
+ USE EOS_Omni_Module, only: press_gf, inv_rho_gf, poly_k_cgs, rho_gf
+#else
USE EOS_Polytrope_Scalars
+#endif
implicit none
@@ -53,16 +57,39 @@
#include "EOS_Base.inc"
!!$ Set up the fluid constants
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: poly_eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ poly_eoskey = GRHydro_poly_eoskey
+! end EOS Omni vars
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ 1.0d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ 1.0d0,1.0d0,xtemp,xye,xpress,xeps,keyerr,anyerr)
+ local_Gamma = 1.0d0 + xpress/xeps
+ press = press_gf * poly_k_cgs * &
+ (rho * inv_rho_gf)**local_Gamma
+ eps = press / (rho * (local_Gamma - 1.d0))
+#else
local_Gamma = 1.d0 + EOS_Pressure(GRHydro_polytrope_handle, 1.d0, 1.d0) / &
EOS_SpecificIntEnergy(GRHydro_polytrope_handle, 1.d0, 1.d0)
-
-!!$ Change the pressure and specific internal energy
-
press = p_geom_factor * eos_k_cgs * &
(rho * rho_geom_factor_inv)**local_Gamma
eps = press / (rho * (local_Gamma - 1.d0))
+#endif
+!!$ Change the pressure and specific internal energy
+
!!$ Get the conserved variables. Hardwired polytrope EoS!!!
!!$ Note that this call also sets pressure and eps
@@ -125,11 +152,33 @@
#include "EOS_Base.inc"
!!$ Set up the fluid constants
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: poly_eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ poly_eoskey = GRHydro_poly_eoskey
+! end EOS Omni vars
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ 1.0d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ 1.0d0,1.0d0,xtemp,xye,xpress,xeps,keyerr,anyerr)
+ local_Gamma = 1.0d0 + xpress/xeps
+ local_K = xpress
+#else
local_Gamma = 1.d0 + EOS_Pressure(GRHydro_polytrope_handle, 1.d0, 1.d0) / &
EOS_SpecificIntEnergy(GRHydro_polytrope_handle, 1.d0, 1.d0)
local_K = EOS_Pressure(GRHydro_polytrope_handle, 1.d0, 1.d0)
+#endif
if (abs(local_Gamma - 2.d0) < 1.d-10) then
@@ -184,7 +233,11 @@
subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS)
+#if USE_EOS_OMNI
+ USE EOS_Omni_Module, only: press_gf, inv_rho_gf, poly_k_cgs, rho_gf
+#else
USE EOS_Polytrope_Scalars
+#endif
implicit none
@@ -203,7 +256,39 @@
#include "EOS_Base.inc"
!!$ Set up the fluid constants
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: poly_eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ poly_eoskey = GRHydro_poly_eoskey
+! end EOS Omni vars
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ 1.0d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ 1.0d0,1.0d0,xtemp,xye,xpress,xeps,keyerr,anyerr)
+ local_Gamma = 1.0d0 + xpress/xeps
+ local_K = xpress
+
+ eos_k_initial_cgs = initial_k * rho_gf**initial_Gamma / press_gf
+
+ press = (local_Gamma - 1.d0) / (initial_Gamma - 1.0d0 ) * press_gf * eos_k_initial_cgs * &
+ (rho * rho_gf) ** initial_Gamma
+
+ eps = press_gf * eos_k_initial_cgs * &
+ (rho * inv_rho_gf) ** initial_Gamma / &
+ (rho * (initial_Gamma - 1.0d0))
+
+#else
+
call CCTK_INFO("Pulling the rug...by changing K and Gamma")
local_Gamma = 1.d0 + EOS_Pressure(GRHydro_polytrope_handle, 1.d0, 1.d0) / &
@@ -220,6 +305,8 @@
(rho * rho_geom_factor_inv) ** initial_Gamma / &
(rho * (initial_Gamma - 1.0d0))
+#endif
+
do k = 1, cctk_lsh(3)
do j = 1, cctk_lsh(2)
do i = 1, cctk_lsh(1)
File [modified]: GRHydro_Marquina.F90
Delta lines: +1 -7
===================================================================
--- trunk/src/GRHydro_Marquina.F90 2010-07-02 21:51:08 UTC (rev 131)
+++ trunk/src/GRHydro_Marquina.F90 2010-07-02 21:57:07 UTC (rev 132)
@@ -37,10 +37,6 @@
implicit none
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
@@ -220,7 +216,6 @@
endif
!!$ END: Check for NaN value (1st check)
-!!$ pressp = EOS_Pressure(GRHydro_eos_handle,primp(1),primp(5))
!!$right state
@@ -244,7 +239,6 @@
endif
!!$ END: Check for NaN value (2nd check)
-!!$ pressm_i = EOS_Pressure(GRHydro_eos_handle,primm_i(1),primm_i(5))
!!$eigenvalues and right eigenvectors
@@ -353,7 +347,7 @@
syflux(i,j,k) = f_marquina(3)
szflux(i,j,k) = f_marquina(4)
tauflux(i,j,k) = f_marquina(5)
-
+
enddo
enddo
enddo
File [modified]: GRHydro_Prim2Con.F90
Delta lines: +46 -3
===================================================================
--- trunk/src/GRHydro_Prim2Con.F90 2010-07-02 21:51:08 UTC (rev 131)
+++ trunk/src/GRHydro_Prim2Con.F90 2010-07-02 21:57:07 UTC (rev 132)
@@ -132,6 +132,21 @@
#include "EOS_Base.inc"
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ eoskey = GRHydro_eoskey
+! end EOS Omni vars
+#endif
+
w = 1.d0 / sqrt(1.d0 - (gxx*dvelx*dvelx + gyy*dvely*dvely + gzz &
*dvelz*dvelz + 2*gxy*dvelx*dvely + 2*gxz*dvelx *dvelz + 2*gyz&
*dvely*dvelz))
@@ -145,7 +160,12 @@
endif
!!$ END: Check for NaN value
+#if USE_EOS_OMNI
+ call EOS_Omni_press(eoskey,keytemp,n,&
+ drho,deps,xtemp,xye,dpress,keyerr,anyerr)
+#else
dpress = EOS_Pressure(handle, drho, deps)
+#endif
vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
@@ -339,6 +359,21 @@
#undef _EOS_BASE_INC_
#endif
#include "EOS_Base.inc"
+
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: poly_eoskey = 0
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+ poly_eoskey = GRHydro_poly_eoskey
+! end EOS Omni vars
+#endif
w_tmp = gxx*dvelx*dvelx + gyy*dvely*dvely + gzz *dvelz*dvelz + &
2*gxy*dvelx*dvely + 2*gxz*dvelx*dvelz + 2*gyz*dvely*dvelz
@@ -357,11 +392,19 @@
else
w = 1.d0 / sqrt(1.d0 - w_tmp)
endif
-
+
+#if USE_EOS_OMNI
+ call EOS_Omni_press(poly_eoskey,keytemp,n,&
+ drho,xeps,xtemp,xye,dpress,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,&
+ drho,xeps,xtemp,xye,dpress,deps,keyerr,anyerr)
+#else
if (handle .ge. 0) then
- dpress = EOS_Pressure(handle, drho, deps)
- deps = EOS_SpecificIntEnergy(handle, drho, dpress)
+ dpress = EOS_Pressure(handle, drho, deps)
+ deps = EOS_SpecificIntEnergy(handle, drho, dpress)
end if
+#endif
vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
