From d9199f1361e1a22b63e59f11dcf575fc4b6fcfc4 Mon Sep 17 00:00:00 2001 From: cott Date: Fri, 15 Mar 2013 21:50:51 +0000 Subject: * add possibility to stich on a polytrope to the lower end of a tabulated cold EOS git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/trunk@76 8e189c6b-2ab8-4400-aa02-70a9cfce18b9 --- src/EOS_Omni_Module.F90 | 2 ++ src/EOS_Omni_SingleVarCalls.F90 | 28 +++++++++++++++++----------- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/EOS_Omni_Module.F90 b/src/EOS_Omni_Module.F90 index 04dd10d..908614e 100644 --- a/src/EOS_Omni_Module.F90 +++ b/src/EOS_Omni_Module.F90 @@ -33,6 +33,8 @@ module EOS_Omni_Module real*8 :: coldeos_gammath = 0.0d0 real*8 :: coldeos_rhomin = 0.0d0 real*8 :: coldeos_rhomax = 0.0d0 + real*8 :: coldeos_low_kappa = 0.0d0 + real*8 :: coldeos_low_gamma = 0.0d0 real*8 :: coldeos_kappa = 0.0d0 real*8 :: coldeos_thfac = 1.0d0 real*8 :: coldeos_dlrho = 1.0d0 diff --git a/src/EOS_Omni_SingleVarCalls.F90 b/src/EOS_Omni_SingleVarCalls.F90 index 350e622..f789806 100644 --- a/src/EOS_Omni_SingleVarCalls.F90 +++ b/src/EOS_Omni_SingleVarCalls.F90 @@ -122,8 +122,9 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,& ! cold tabular EOS with gamma law do i=1,npoints if(rho(i).lt.coldeos_rhomin) then - keyerr(i) = 104 - anyerr = 1 + press(i) = coldeos_low_kappa * rho(i)**coldeos_low_gamma + eps(i) = press(i) / (coldeos_low_gamma - 1.0) / rho(i) + cycle else if(rho(i).gt.coldeos_rhomax) then keyerr(i) = 103 anyerr = 1 @@ -281,8 +282,9 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,& ! and entropy (the latter, because T=0) do i=1,npoints if(rho(i).lt.coldeos_rhomin) then - keyerr(i) = 104 - anyerr = 1 + dpdrhoe(i) = coldeos_low_kappa * coldeos_low_gamma * & + rho(i)**(coldeos_low_gamma-1.0d0) + cycle else if(rho(i).gt.coldeos_rhomax) then keyerr(i) = 103 anyerr = 1 @@ -421,8 +423,8 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,& ! only the gamma law has non-zero dPdeps do i=1,npoints if(rho(i).lt.coldeos_rhomin) then - keyerr(i) = 104 - anyerr = 1 + dpdepsrho(i) = 0.0d0 + cycle else if(rho(i).gt.coldeos_rhomax) then keyerr(i) = 103 anyerr = 1 @@ -573,13 +575,17 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,& enddo case (5) - ! with the cold eos we have to assume P = P(rho), so - ! by definition dPdrho is at constant internal energy - ! and entropy (the latter, because T=0) + ! with the cold eos we have to assume P = P(rho), so + ! by definition dPdrho is at constant internal energy + ! and entropy (the latter, because T=0) do i=1,npoints if(rho(i).lt.coldeos_rhomin) then - keyerr(i) = 104 - anyerr = 1 + xprs = coldeos_low_kappa * rho(i)**coldeos_low_gamma + eps(i) = xprs / (coldeos_low_gamma - 1.0d0) / rho(i) + cs2(i) = coldeos_low_kappa * coldeos_low_gamma * & + rho(i)**(coldeos_low_gamma-1.0d0) / & + (1.0d0 + eps(i) + xprs) + cycle else if(rho(i).gt.coldeos_rhomax) then keyerr(i) = 103 anyerr = 1 -- cgit v1.2.3