diff options
author | cott <cott@8e189c6b-2ab8-4400-aa02-70a9cfce18b9> | 2010-12-05 00:46:52 +0000 |
---|---|---|
committer | cott <cott@8e189c6b-2ab8-4400-aa02-70a9cfce18b9> | 2010-12-05 00:46:52 +0000 |
commit | cb4a0f15f66ea5a14554dc2699a002a480e9e82a (patch) | |
tree | cb7a6231a2c5f602ab2e64f1382d5794c785003c | |
parent | 0329ffcbf0b395c9914d2acfa73dc9e67d048871 (diff) |
* updates needed to make TOVSolverHot work
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/trunk@27 8e189c6b-2ab8-4400-aa02-70a9cfce18b9
-rw-r--r-- | interface.ccl | 7 | ||||
-rw-r--r-- | src/EOS_Omni_SingleVarCalls.F90 | 70 | ||||
-rw-r--r-- | src/nuc_eos/bisection.F90 | 85 | ||||
-rw-r--r-- | src/nuc_eos/make.code.defn | 2 | ||||
-rw-r--r-- | src/nuc_eos/nuc_eos.F90 | 103 |
5 files changed, 208 insertions, 59 deletions
diff --git a/interface.ccl b/interface.ccl index 1903bcd..55d6d58 100644 --- a/interface.ccl +++ b/interface.ccl @@ -100,19 +100,20 @@ void FUNCTION EOS_Omni_EpsFromPress(CCTK_INT IN eoskey, \ 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 diff --git a/src/EOS_Omni_SingleVarCalls.F90 b/src/EOS_Omni_SingleVarCalls.F90 index 72abe97..7f1f825 100644 --- a/src/EOS_Omni_SingleVarCalls.F90 +++ b/src/EOS_Omni_SingleVarCalls.F90 @@ -509,38 +509,43 @@ subroutine EOS_Omni_EOS_eps_from_press(eoskey,keytemp,rf_precision,npoints,& 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 @@ subroutine EOS_Omni_EOS_RestMassDensityFromEpsPress(eoskey,keytemp,rf_precision, 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 diff --git a/src/nuc_eos/bisection.F90 b/src/nuc_eos/bisection.F90 index f12b702..fb0eb43 100644 --- a/src/nuc_eos/bisection.F90 +++ b/src/nuc_eos/bisection.F90 @@ -7,7 +7,7 @@ subroutine bisection(lr,lt0,y,eps0,lt,bivar,keyerrt,keybisect) 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 @@ subroutine bisection(lr,lt0,y,eps0,lt,bivar,keyerrt,keybisect) integer bcount,i,itmax,maxbcount - real*8 bivar(*) + real*8 bivar bcount = 0 maxbcount = 80 @@ -80,3 +80,84 @@ subroutine bisection(lr,lt0,y,eps0,lt,bivar,keyerrt,keybisect) 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 diff --git a/src/nuc_eos/make.code.defn b/src/nuc_eos/make.code.defn index 59c8118..1b5de60 100644 --- a/src/nuc_eos/make.code.defn +++ b/src/nuc_eos/make.code.defn @@ -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 = diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90 index 399d89b..51170d0 100644 --- a/src/nuc_eos/nuc_eos.F90 +++ b/src/nuc_eos/nuc_eos.F90 @@ -24,7 +24,7 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& 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 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& ! 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 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& 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 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& 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 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& 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 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& 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 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& 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 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& 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 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& 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 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& 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 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& 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 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& 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 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& 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 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& 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 |