diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/GRHydro_Con2Prim.F90 | 249 | ||||
| -rw-r--r-- | src/GRHydro_Eigenproblem.F90 | 92 | ||||
| -rw-r--r-- | src/GRHydro_Eigenproblem_Marquina.F90 | 49 | ||||
| -rw-r--r-- | src/GRHydro_EoSChangeGamma.F90 | 93 | ||||
| -rw-r--r-- | src/GRHydro_Marquina.F90 | 8 | ||||
| -rw-r--r-- | src/GRHydro_Prim2Con.F90 | 49 |
6 files changed, 506 insertions, 34 deletions
diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90 index 8e13501..c38d412 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -150,6 +150,21 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) 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 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) 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 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) 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,8 +266,8 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) 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 @@ -267,8 +298,16 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) 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 @@ a single point. 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 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & #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,8 +403,14 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & 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 @@ -386,7 +447,15 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & !!$ 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 @@ -430,8 +499,18 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & 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 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & 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 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & 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,7 +579,15 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & 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 @@ -575,14 +678,42 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS) 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 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) 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 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & #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 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & 2.d0*usx*usz*uxz + 2.d0*usy*usz*uyz !!$ Set initial guess for 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 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & 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 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & 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 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & !!$ 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 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & 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 - w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 )) + press = EOS_Pressure(handle, rhonew, epsilon) +#endif - press = EOS_Pressure(handle, rhonew, epsilon) + w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 )) tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens f = rhonew * w_lorentz - udens @@ -1022,10 +1209,21 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & 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 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & 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 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & 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 @@ end subroutine Con2Prim_ptPolytype @@*/ + subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS) implicit none diff --git a/src/GRHydro_Eigenproblem.F90 b/src/GRHydro_Eigenproblem.F90 index 194ba72..d98e2d8 100644 --- a/src/GRHydro_Eigenproblem.F90 +++ b/src/GRHydro_Eigenproblem.F90 @@ -63,16 +63,49 @@ subroutine eigenvalues(handle,rho,velx,vely,velz,eps,w_lorentz,& #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 @@ subroutine eigenvalues(handle,rho,velx,vely,velz,eps,w_lorentz,& lam(4) = lam3 lam(5) = lamp + end subroutine eigenvalues /*@@ @@ -157,16 +191,47 @@ subroutine eigenproblem(handle,rho,velx,vely,velz,eps,& #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 @@ subroutine eigenproblem(handle,rho,velx,vely,velz,eps,& 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 @@ subroutine eigenproblem(handle,rho,velx,vely,velz,eps,& roeflux4 = rflux(4) roeflux5 = rflux(5) + if(roeflux1.ne.roeflux1) then + stop "bad herberd" + endif + end subroutine eigenproblem /*@@ @@ -522,16 +592,38 @@ subroutine eigenproblem_leftright(handle,rho,velx,vely,velz,eps,& #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 diff --git a/src/GRHydro_Eigenproblem_Marquina.F90 b/src/GRHydro_Eigenproblem_Marquina.F90 index 443329e..806f952 100644 --- a/src/GRHydro_Eigenproblem_Marquina.F90 +++ b/src/GRHydro_Eigenproblem_Marquina.F90 @@ -99,19 +99,46 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,& #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 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,& !!$ 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 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,& kappal = dpdepsl / (dpdepsl - rhol * cs2l) + !!$ Right eigenvector # 1 reivec1l(1) = kappal / (enthalpyl * wl) diff --git a/src/GRHydro_EoSChangeGamma.F90 b/src/GRHydro_EoSChangeGamma.F90 index 5c854a6..f486c79 100644 --- a/src/GRHydro_EoSChangeGamma.F90 +++ b/src/GRHydro_EoSChangeGamma.F90 @@ -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,15 +57,38 @@ subroutine GRHydro_EoSChangeGamma(CCTK_ARGUMENTS) #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 @@ subroutine GRHydro_EoSChangeK(CCTK_ARGUMENTS) #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 @@ end subroutine GRHydro_EoSChangeK 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,6 +256,38 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS) #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") @@ -220,6 +305,8 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS) (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) diff --git a/src/GRHydro_Marquina.F90 b/src/GRHydro_Marquina.F90 index a4196cf..144049f 100644 --- a/src/GRHydro_Marquina.F90 +++ b/src/GRHydro_Marquina.F90 @@ -37,10 +37,6 @@ subroutine GRHydro_Marquina(CCTK_ARGUMENTS) 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 @@ subroutine GRHydro_Marquina(CCTK_ARGUMENTS) endif !!$ END: Check for NaN value (1st check) -!!$ pressp = EOS_Pressure(GRHydro_eos_handle,primp(1),primp(5)) !!$right state @@ -244,7 +239,6 @@ subroutine GRHydro_Marquina(CCTK_ARGUMENTS) 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 @@ subroutine GRHydro_Marquina(CCTK_ARGUMENTS) syflux(i,j,k) = f_marquina(3) szflux(i,j,k) = f_marquina(4) tauflux(i,j,k) = f_marquina(5) - + enddo enddo enddo diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90 index 9710a60..06ed464 100644 --- a/src/GRHydro_Prim2Con.F90 +++ b/src/GRHydro_Prim2Con.F90 @@ -132,6 +132,21 @@ subroutine prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & #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 @@ subroutine prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & 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 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, & #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 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, & 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 |