[Commits] [svn:einsteintoolkit] EOS_Omni/trunk/ (Rev. 27)
cott at tapir.caltech.edu
cott at tapir.caltech.edu
Sat Dec 4 18:46:52 CST 2010
User: cott
Date: 2010/12/04 06:46 PM
Modified:
/trunk/
interface.ccl
/trunk/src/
EOS_Omni_SingleVarCalls.F90
/trunk/src/nuc_eos/
bisection.F90, make.code.defn, nuc_eos.F90
Log:
* updates needed to make TOVSolverHot work
File Changes:
Directory: /trunk/src/nuc_eos/
==============================
File [modified]: bisection.F90
Delta lines: +83 -2
===================================================================
--- trunk/src/nuc_eos/bisection.F90 2010-11-26 19:45:37 UTC (rev 26)
+++ trunk/src/nuc_eos/bisection.F90 2010-12-05 00:46:52 UTC (rev 27)
@@ -7,7 +7,7 @@
real*8 lr,lt0,y,eps0,lt
integer keyerrt
- integer keybisect
+ integer keybisect !this doesn't do anything...
! 1 -> operate on energy
! 2 -> operate on entropy
@@ -19,7 +19,7 @@
integer bcount,i,itmax,maxbcount
- real*8 bivar(*)
+ real*8 bivar
bcount = 0
maxbcount = 80
@@ -80,3 +80,84 @@
end subroutine bisection
+
+subroutine bisection_rho(lr0,lt,y,lpressin,lr,bivar,keyerrr,keybisect)
+
+ use eosmodule
+
+ implicit none
+
+ real*8 lr,lr0,y,lpressin,lt
+ integer keyerrr
+
+ integer keybisect !this doesn't do anything
+ ! 3 -> operate on pressure
+
+ !temporary vars
+ real*8 lr1,lr2,lrmin,lrmax
+ real*8 f1,f2,fmid,dlr,lrmid
+ real*8 d1,d2,d3,tol
+ real*8 f1a,f2a
+
+ integer bcount,i,itmax,maxbcount
+
+ real*8 bivar
+
+ bcount = 0
+ maxbcount = 80
+
+ tol=1.d-9 ! need to find energy to less than 1 in 10^-9
+ itmax=50
+
+ lrmax=logrho(nrho)
+ lrmin=logrho(1)
+
+ lr = lr0
+ lr1 = dlog10(min(10.0d0**lrmax,1.10d0*(10.0d0**lr0)))
+ lr2 = dlog10(max(10.0d0**lrmin,0.90d0*(10.0d0**lr0)))
+
+ call findthis(lr1,lt,y,f1a,bivar,d1,d2,d3)
+ call findthis(lr2,lt,y,f2a,bivar,d1,d2,d3)
+
+ f1=f1a-lpressin
+ f2=f2a-lpressin
+
+ keyerrr=0
+
+ do while(f1*f2.ge.0.0d0)
+ bcount=bcount+1
+ lr1=dlog10(min(10.0d0**lrmax,1.2d0*(10.0d0**lr1)))
+ lr2=dlog10(max(10.0d0**lrmin,0.8d0*(10.0d0**lr2)))
+ call findthis(lr1,lt,y,f1a,bivar,d1,d2,d3)
+ call findthis(lr2,lt,y,f2a,bivar,d1,d2,d3)
+ f1=f1a-lpressin
+ f2=f2a-lpressin
+ if(bcount.ge.maxbcount) then
+ keyerrr=667
+ return
+ endif
+ enddo
+
+ if(f1.lt.0.0d0) then
+ lr=lr1
+ dlr=dlog10((10.0D0**lr2)-(10.0d0**lr1))
+ else
+ lr=lr2
+ dlr=dlog10((10.0D0**lr1)-(10.0d0**lr2))
+ endif
+
+ do i=1,itmax
+ dlr=dlog10((10.0d0**dlr)*0.5d0)
+ lrmid=dlog10(10.0d0**lr+10.0d0**dlr)
+ call findthis(lrmid,lt,y,f2a,bivar,d1,d2,d3)
+ fmid=f2a-lpressin
+ if(fmid.le.0.0d0) lr=lrmid
+ if(abs(1.0d0-f2a/lpressin).lt.tol) then
+ lr=lrmid
+ return
+ endif
+ enddo
+
+
+
+end subroutine bisection_rho
File [modified]: make.code.defn
Delta lines: +1 -1
===================================================================
--- trunk/src/nuc_eos/make.code.defn 2010-11-26 19:45:37 UTC (rev 26)
+++ trunk/src/nuc_eos/make.code.defn 2010-12-05 00:46:52 UTC (rev 27)
@@ -1,5 +1,5 @@
SRCS = eosmodule.F90 nuc_eos.F90 bisection.F90 \
- findtemp.F90 linterp_many.F90 readtable.F90 \
+ findtemp.F90 findrho.F90 linterp_many.F90 readtable.F90 \
linterp.f
SUBDIRS =
File [modified]: nuc_eos.F90
Delta lines: +67 -36
===================================================================
--- trunk/src/nuc_eos/nuc_eos.F90 2010-11-26 19:45:37 UTC (rev 26)
+++ trunk/src/nuc_eos/nuc_eos.F90 2010-12-05 00:46:52 UTC (rev 27)
@@ -24,7 +24,7 @@
use eosmodule
implicit none
- real*8, intent(in) :: xrho,xye
+ real*8, intent(inout) :: xrho,xye
real*8, intent(inout) :: xtemp,xenr,xent
real*8, intent(out) :: xprs,xcs2,xdedt
real*8, intent(out) :: xdpderho,xdpdrhoe,xxa,xxh,xxn,xxp
@@ -35,10 +35,11 @@
! local variables
- real*8 :: lr,lt,y,xx,xeps,leps,xs
+ real*8 :: lr,lt,y,xx,xeps,leps,xs,xpressure
real*8 :: d1,d2,d3
real*8 :: ff(nvars)
integer :: keyerrt = 0
+ integer :: keyerrr = 0
if(xrho.gt.eos_rhomax) then
stop "nuc_eos: rho > rhomax"
@@ -49,11 +50,13 @@
endif
if(xye.gt.eos_yemax) then
- stop "nuc_eos: ye > yemax"
+ keyerr = 101
+ return
endif
if(xye.lt.eos_yemin) then
- stop "nuc_eos: ye < yemin"
+ keyerr = 102
+ return
endif
if(keytemp.eq.1) then
@@ -87,19 +90,33 @@
xs = xent
call findtemp_entropy(lr,lt,y,xs,keyerrt,rfeps)
xtemp = 10.0d0**lt
+
+ elseif(keytemp.eq.3) then
+ !need to find rho based on xprs
+ xpressure = log10(xprs)
+ call findrho_press(lr,lt,y,xpressure,keyerrr,rfeps)
+ if(keyerrr.ne.0) then
+ write(*,*) "Problem in findrho_press:", keyerr
+ keyerr = keyerrr
+ return
+ endif
+ xrho = 10.0d0**lr
+
endif
! have temperature, proceed:
call findall(lr,lt,y,ff)
- xprs = 10.0d0**ff(1)
+ !unless we want xprs to be constant (keytemp==3), reset xprs
+ if(.not.keytemp.eq.3) then
+ xprs = 10.0d0**ff(1)
+ endif
- xeps = 10.0d0**ff(2) - energy_shift
- if(keytemp.eq.1.or.keytemp.eq.2) then
- xenr = xeps
+ if(.not.keytemp.eq.0) then
+ xenr = 10.0d0**ff(2) - energy_shift
endif
- if(keytemp.eq.1.or.keytemp.eq.0) then
+ if(.not.keytemp.eq.2) then
xent = ff(3)
endif
@@ -167,7 +184,7 @@
use eosmodule
implicit none
- real*8, intent(in) :: xrho,xye
+ real*8, intent(inout) :: xrho,xye
real*8, intent(inout) :: xtemp,xenr,xent
real*8, intent(in) :: rfeps
real*8, intent(out) :: xprs,xmunu,xcs2,xdedt
@@ -176,9 +193,10 @@
integer, intent(out) :: keyerr
! local variables
- real*8 :: lr,lt,y,xx,xeps,leps,xs
+ real*8 :: lr,lt,y,xx,xeps,leps,xs,xpressure
real*8 :: d1,d2,d3,ff(8)
integer :: keyerrt = 0
+ integer :: keyerrr = 0
if(xrho.gt.eos_rhomax) then
stop "nuc_eos: rho > rhomax"
@@ -191,13 +209,13 @@
endif
if(xye.gt.eos_yemax) then
- stop "nuc_eos: ye > yemax"
+ keyerr = 101
+ return
endif
if(xye.lt.eos_yemin) then
keyerr = 102
- write(6,*) "ye: ", xye
- stop "nuc_eos ye < yemin"
+ return
endif
if(keytemp.eq.1) then
@@ -220,7 +238,7 @@
keyerr = 0
- if(keytemp.eq.0) then
+if(keytemp.eq.0) then
!need to find temperature based on xeps
call findtemp(lr,lt,y,leps,keyerrt,rfeps)
if(keyerrt.ne.0) then
@@ -239,18 +257,33 @@
endif
xtemp = 10.0d0**lt
+ elseif(keytemp.eq.3) then
+ !need to find rho based on xprs
+ xpressure = log10(xprs)
+ call findrho_press(lr,lt,y,xpressure,keyerrr,rfeps)
+ if (keyerrr.ne.0) then
+ keyerr = keyerrr
+ write(*,*) "Problem in findrho_press:", keyerr
+ return
+ endif
+ xrho = 10.0d0**lr
endif
- ! have temperature, proceed:
+ ! have rho,temp,ye; proceed:
call findall_short(lr,lt,y,ff)
- xprs = 10.0d0**ff(1)
- xeps = 10.0d0**ff(2) - energy_shift
- if((keytemp.eq.1).or.(keytemp.eq.2)) then
- xenr = xeps
+ !unless we want xprs to be constant (keytemp==3), reset xprs
+ if(.not.keytemp.eq.3) then
+ xprs = 10.0d0**ff(1)
endif
- if(keytemp.eq.1.or.keytemp.eq.0) then
+ !unless we want xenr to be constant (keytemp==0), reset xenr
+ if(.not.keytemp.eq.0) then
+ xenr = 10.0d0**ff(2) - energy_shift
+ endif
+
+ !unless we want xent to be constant (keytemp==2), reset xent
+ if(.not.keytemp.eq.2) then
xent = ff(3)
endif
@@ -296,13 +329,13 @@
endif
if(xye.gt.eos_yemax) then
- stop "nuc_eos: ye > yemax"
+ keyerr = 101
+ return
endif
if(xye.lt.eos_yemin) then
keyerr = 102
- write(6,*) "ye: ", xye
- stop "nuc_eos ye < yemin"
+ return
endif
if(keytemp.eq.1) then
@@ -333,17 +366,16 @@
endif
xtemp = 10.0d0**lt
- elseif(keytemp.eq.2) then
- stop "eos_nuc_press does not support keytemp.eq.2"
+ elseif(keytemp.gt.1) then
+ stop "eos_nuc_press does not support keytemp other than 0 and 1"
endif
! have temperature, proceed:
call findall_press_eps(lr,lt,y,ff)
xprs = 10.0d0**ff(1)
- xeps = 10.0d0**ff(2) - energy_shift
- if((keytemp.eq.1).or.(keytemp.eq.2)) then
- xenr = xeps
+ if(keytemp.eq.1) then
+ xenr = 10.0d0**ff(2) - energy_shift
endif
end subroutine nuc_eos_press_eps
@@ -378,13 +410,13 @@
endif
if(xye.gt.eos_yemax) then
- stop "nuc_eos: ye > yemax"
+ keyerr = 101
+ return
endif
if(xye.lt.eos_yemin) then
keyerr = 102
- write(6,*) "ye: ", xye
- stop "nuc_eos ye < yemin"
+ return
endif
if(keytemp.eq.1) then
@@ -415,8 +447,8 @@
endif
xtemp = 10.0d0**lt
- elseif(keytemp.eq.2) then
- stop "eos_nuc_press does not support keytemp.eq.2"
+ elseif(keytemp.gt.1) then
+ stop "eos_nuc_press does not support keytemp > 1"
endif
! have temperature, proceed:
@@ -424,8 +456,7 @@
xdpdrhoe = 10.0d0**ff(1)
xdpderho = 10.0d0**ff(2)
-
- if((keytemp.eq.1).or.(keytemp.eq.2)) then
+ if(keytemp.eq.1) then
call findthis(lr,lt,y,xeps,alltables(:,:,:,2),d1,d2,d3)
xeps = 10.0d0**xeps - energy_shift
xenr = xeps
Directory: /trunk/src/
======================
File [modified]: EOS_Omni_SingleVarCalls.F90
Delta lines: +53 -17
===================================================================
--- trunk/src/EOS_Omni_SingleVarCalls.F90 2010-11-26 19:45:37 UTC (rev 26)
+++ trunk/src/EOS_Omni_SingleVarCalls.F90 2010-12-05 00:46:52 UTC (rev 27)
@@ -509,38 +509,43 @@
end subroutine EOS_Omni_EOS_eps_from_press
-subroutine EOS_Omni_EOS_RestMassDensityFromEpsPress(eoskey,keytemp,rf_precision,&
- npoints,rho,eps,temp,ye,press,keyerr,anyerr)
+subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt(eoskey,&
+ ikeytemp,rf_precision,&
+ npoints,rho,eps,temp,entropy,ye,press,keyerr,anyerr)
use EOS_Omni_Module
implicit none
DECLARE_CCTK_PARAMETERS
- CCTK_INT, intent(in) :: eoskey,keytemp,npoints
+ CCTK_INT, intent(in) :: eoskey,ikeytemp,npoints
CCTK_INT, intent(out) :: keyerr(npoints)
CCTK_INT, intent(out) :: anyerr
CCTK_REAL, intent(in) :: rf_precision
- CCTK_REAL, intent(in) :: ye(npoints),press(npoints),eps(npoints)
- CCTK_REAL, intent(inout) :: rho(npoints),temp(npoints)
+ CCTK_REAL, intent(in) :: ye(npoints),press(npoints)
+ CCTK_REAL, intent(inout) :: rho(npoints),temp(npoints),eps(npoints)
+ CCTK_REAL, intent(inout) :: entropy(npoints)
! local vars
integer :: i
character(256) :: warnstring
+ integer :: keytemp
+ ! temporary vars for nuc_eos
+ real*8 :: xrho,xye,xtemp,xenr,xent
+ real*8 :: xprs,xmunu,xcs2
+ real*8 :: xdedt,xdpderho,xdpdrhoe
+
+ keytemp = ikeytemp
- if(keytemp.eq.1) then
- anyerr = 1
- keyerr(:) = -1
- else
- anyerr = 0
- keyerr(:) = 0
- endif
+ anyerr = 0
+ keyerr(:) = 0
select case (eoskey)
case (1)
! polytropic EOS
do i=1,npoints
- rho(i) = press(i) / ((poly_gamma - 1.0d0)*eps(i))
+ rho(i) = (press(i) / poly_k)**(1.0/poly_gamma)
+ eps(i) = press(i) / (poly_gamma - 1.0d0) / rho(i)
enddo
case (2)
! gamma-law EOS
@@ -549,15 +554,46 @@
enddo
case (3)
! hybrid EOS
- write(warnstring,*) "EOS_Omni_RestMassDensityFromEpsPress not supported for hybrid EOS"
+ write(warnstring,*) "EOS_Omni_RestMassDensityFromPressEpsTemp not supported for hybrid EOS"
call CCTK_WARN(0,warnstring)
case (4)
! nuc EOS
- write(warnstring,*) "EOS_Omni_RestMassDensityFromEpsPress call not supported for nuc_eos yet"
- call CCTK_WARN(0,warnstring)
+ if(ikeytemp.eq.2) then
+ call CCTK_WARN(0,"This function does not work yet when coming in with entropy")
+ else if(ikeytemp.eq.1) then
+ keytemp = 3
+ else
+ call CCTK_WARN(0,"This function does not work yet when coming in with this keytemp")
+ endif
+
+ do i=1,npoints
+ ! initial guess
+ xprs = press(i) * inv_press_gf
+ xrho = rho(i) * inv_rho_gf
+ xtemp = temp(i)
+ xent = entropy(i)
+ xye = ye(i)
+ xenr = eps(i) * inv_eps_gf
+ call nuc_eos_short(xrho,xtemp,xye,xenr,xprs,&
+ xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,&
+ keytemp,keyerr(i),rf_precision)
+
+ if(keyerr(i).ne.0) then
+ anyerr = 1
+ endif
+
+ if(keytemp.eq.3) then
+ eps(i) = xenr * eps_gf
+ rho(i) = xrho * rho_gf
+ entropy(i) = xent
+ else
+ call CCTK_WARN(0, "This keytemp is not implemented yet")
+ endif
+
+ enddo
case DEFAULT
write(warnstring,*) "eoskey ",eoskey," not implemented!"
call CCTK_WARN(0,warnstring)
end select
-end subroutine EOS_Omni_EOS_RestMassDensityFromEpsPress
+end subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt
Directory: /trunk/
==================
File [modified]: interface.ccl
Delta lines: +4 -3
===================================================================
--- trunk/interface.ccl 2010-11-26 19:45:37 UTC (rev 26)
+++ trunk/interface.ccl 2010-12-05 00:46:52 UTC (rev 27)
@@ -100,19 +100,20 @@
PROVIDES FUNCTION EOS_Omni_EpsFromPress WITH EOS_Omni_EOS_eps_from_press LANGUAGE Fortran
-void FUNCTION EOS_Omni_RestMassDensityFromEpsPress(CCTK_INT IN eoskey, \
+void FUNCTION EOS_Omni_RhoFromPressEpsTempEnt(CCTK_INT IN eoskey, \
CCTK_INT IN havetemp, \
CCTK_REAL IN rf_precision, \
CCTK_INT IN npoints, \
CCTK_REAL OUT ARRAY rho, \
- CCTK_REAL IN ARRAY eps, \
+ CCTK_REAL INOUT ARRAY eps, \
CCTK_REAL INOUT ARRAY temp, \
+ CCTK_REAL INOUT ARRAY ent, \
CCTK_REAL IN ARRAY ye, \
CCTK_REAL IN ARRAY press, \
CCTK_INT OUT ARRAY keyerr, \
CCTK_INT OUT anyerr)
-PROVIDES FUNCTION EOS_Omni_RestMassDensityFromEpsPress WITH EOS_Omni_EOS_RestMassDensityFromEpsPress LANGUAGE Fortran
+PROVIDES FUNCTION EOS_Omni_RhoFromPressEpsTempEnt WITH EOS_Omni_EOS_RhoFromPressEpsTempEnt LANGUAGE Fortran
################################################################################
# short and long and specific calls
More information about the Commits
mailing list