* updates needed to make TOVSolverHot work

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/trunk@27 8e189c6b-2ab8-4400-aa02-70a9cfce18b9

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