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