[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 132)

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