From 8682216ce83fb74500d4318577e9955b8dd1281e Mon Sep 17 00:00:00 2001 From: cott Date: Mon, 7 Feb 2011 05:36:36 +0000 Subject: * update nuc_eos with some changes in the low-density regime. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/trunk@37 8e189c6b-2ab8-4400-aa02-70a9cfce18b9 --- src/nuc_eos/nuc_eos.F90 | 45 ++++++++++++++++++++++++++++++++------------- 1 file changed, 32 insertions(+), 13 deletions(-) diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90 index b576ce0..64a78a5 100644 --- a/src/nuc_eos/nuc_eos.F90 +++ b/src/nuc_eos/nuc_eos.F90 @@ -160,6 +160,25 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& end subroutine nuc_eos_full +subroutine nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + + implicit none + real*8 xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe + real*8,parameter:: polyK1 = 5.0d14 + real*8,parameter :: polyGamma = 1.33333333333333d0 + integer keytemp + + xprs=polyK1*xrho**(polyGamma) + xenr=xprs/xrho/(polyGamma-1.d0) + xcs2=polyGamma*xprs/xrho + xdpderho=0.0d0 + xdpdrhoe=0.0d0 + + xenr=max(8.0d18,xenr) + xenr=min(3.0d20,xenr) + +end subroutine nuc_poly_eos + subroutine nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) implicit none @@ -212,9 +231,9 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& call CCTK_WARN (0, "nuc_eos: rho > rhomax") endif - if(xrho.lt.eos_rhomin*1.2d0) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) - xent = 4.0d0 + if(xrho.lt.eos_rhomin*1.126d0) then + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + xent = 0.0d0 return endif @@ -234,8 +253,8 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& endif if(xtemp.lt.eos_tempmin) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) - xent = 4.0d0 + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + xent = 0.0d0 return endif endif @@ -335,8 +354,8 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& call CCTK_WARN (0, "nuc_eos: rho > rhomax") endif - if(xrho.lt.eos_rhomin*1.2d0) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + if(xrho.lt.eos_rhomin*1.126d0) then + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) return endif @@ -354,9 +373,9 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& if(xtemp.gt.eos_tempmax) then call CCTK_WARN (0, "nuc_eos: temp > tempmax") endif - + if(xtemp.lt.eos_tempmin) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) return endif endif @@ -414,8 +433,8 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& call CCTK_WARN (0, "nuc_eos: rho > rhomax") endif - if(xrho.lt.eos_rhomin*1.2d0) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + if(xrho.lt.eos_rhomin*1.126d0) then + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) return endif @@ -433,9 +452,9 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& if(xtemp.gt.eos_tempmax) then call CCTK_WARN (0, "nuc_eos: temp > tempmax") endif - + if(xtemp.lt.eos_tempmin) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) return endif endif -- cgit v1.2.3