* add a new routine that computes dpde|rho and dpdrho|e

  at once.



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

4 files changed, 244 insertions, 2 deletions

diff --git a/interface.ccl b/interface.ccl
index 90e6151..1903bcd 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -52,6 +52,22 @@ void FUNCTION EOS_Omni_DPressByDRho(CCTK_INT IN eoskey,         \
 
 PROVIDES FUNCTION EOS_Omni_DPressByDRho WITH EOS_Omni_EOS_DPressByDRho LANGUAGE Fortran
 
+void FUNCTION EOS_Omni_dpderho_dpdrhoe(CCTK_INT IN eoskey,         \
+			     CCTK_INT IN havetemp,              \
+			     CCTK_REAL IN rf_precision,         \
+                             CCTK_INT IN npoints,               \
+			     CCTK_REAL IN ARRAY rho,            \
+			     CCTK_REAL INOUT ARRAY eps,         \
+			     CCTK_REAL INOUT ARRAY temp,        \
+			     CCTK_REAL IN ARRAY ye,             \
+			     CCTK_REAL OUT ARRAY dpderho,       \
+			     CCTK_REAL OUT ARRAY dpdrhoe,       \
+			     CCTK_INT OUT ARRAY keyerr,         \
+                             CCTK_INT OUT anyerr)
+
+PROVIDES FUNCTION EOS_Omni_dpderho_dpdrhoe WITH EOS_Omni_EOS_dpderho_dpdrhoe LANGUAGE Fortran
+
+
 
 void FUNCTION EOS_Omni_cs2(CCTK_INT IN eoskey,                  \
 			     CCTK_INT IN havetemp,              \
diff --git a/src/EOS_Omni_MultiVarCalls.F90 b/src/EOS_Omni_MultiVarCalls.F90
index ec35cfe..25be0a2 100644
--- a/src/EOS_Omni_MultiVarCalls.F90
+++ b/src/EOS_Omni_MultiVarCalls.F90
@@ -82,3 +82,123 @@ subroutine EOS_Omni_EOS_short(eoskey,keytemp,rf_precision,npoints,&
   enddo
 
 end subroutine EOS_Omni_EOS_short
+
+
+subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
+     rho,eps,temp,ye,dpderho,dpdrhoe,keyerr,anyerr)
+
+  use EOS_Omni_Module
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT, intent(in)     :: eoskey,keytemp,npoints
+  CCTK_INT, intent(out)    :: keyerr(npoints)
+  CCTK_INT, intent(out)    :: anyerr
+  CCTK_REAL, intent(in)    :: rf_precision
+  CCTK_REAL, intent(in)    :: rho(npoints),ye(npoints)
+  CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints)
+  CCTK_REAL, intent(out)   :: dpderho(npoints)
+  CCTK_REAL, intent(out)   :: dpdrhoe(npoints)
+
+  ! local vars
+  integer          :: i
+  character(256)   :: warnstring
+  real*8           :: hybrid_local_gamma
+  real*8           :: hybrid_local_k_cgs
+  real*8           :: hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2
+  ! temporary vars for nuc_eos
+  real*8           :: xrho,xye,xtemp,xenr,xent
+  real*8           :: xprs,xmunu,xcs2
+  real*8           :: xdedt,xdpderho,xdpdrhoe
+
+
+  anyerr    = 0
+  keyerr(:) = 0
+
+  select case (eoskey)
+  case (1)
+     ! polytropic EOS                                                                                      
+        if(keytemp.eq.1) then
+           do i=1,npoints
+              eps(i) = press_gf * poly_k_cgs * &
+                   (rho(i)*inv_rho_gf)**(poly_gamma) / &
+                   (poly_gamma - 1.0d0) / rho(i)
+           enddo
+        endif
+        do i=1,npoints
+           dpdrhoe(i) = press_gf * poly_k_cgs *  &
+                poly_gamma * inv_rho_gf *        &
+                (rho(i)*inv_rho_gf) ** (poly_gamma - 1.d0)
+           dpderho(i) = 0.0d0
+        enddo
+  case (2)
+     ! gamma-law EOS                                                                                       
+        if(keytemp.eq.1) then
+           do i=1,npoints
+              eps(i) = press_gf * gl_k_cgs * &
+                   (rho(i)*inv_rho_gf)**(gl_gamma) / &
+                   (gl_gamma - 1.0d0) / rho(i)
+           enddo
+        endif
+        do i=1,npoints
+           dpdrhoe(i) = (gl_gamma-1.0d0) * &
+                eps(i)
+           dpderho(i) = (gl_gamma - 1.0d0) * &
+                rho(i)
+        enddo
+  case (3)
+     ! hybrid EOS                                                                                          
+        do i=1,npoints
+           if(rho(i).gt.hybrid_rho_nuc) then
+              hybrid_local_gamma = hybrid_gamma2
+              hybrid_local_k_cgs = hybrid_k2_cgs
+           else
+              hybrid_local_gamma = hybrid_gamma1
+              hybrid_local_k_cgs = hybrid_k1_cgs
+           endif
+           hybrid_dp_poly = hybrid_local_gamma * press_gf *                 &
+                hybrid_local_k_cgs * rho(i)**(hybrid_local_gamma - 1.0d0) * &
+                inv_rho_gf**hybrid_local_gamma
+
+           hybrid_dp_th1 = - hybrid_local_gamma * press_gf * hybrid_local_k_cgs * &
+                (hybrid_gamma_th - 1.d0) / (hybrid_local_gamma - 1.d0) *      &
+                rho(i)**(hybrid_local_gamma - 1.d0) * inv_rho_gf**hybrid_local_gamma
+
+           hybrid_dp_th2  = (hybrid_gamma_th - 1.d0) * eps(i)                        &
+                - (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) /  &
+                (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) *                    &
+                press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 *               &
+                hybrid_rho_nuc**(hybrid_gamma1 - 1.d0)
+
+           dpdrhoe(i) = hybrid_dp_poly + hybrid_dp_th1 + hybrid_dp_th2
+           dpderho(i) = (hybrid_gamma_th - 1.0d0) * rho(i)
+        enddo
+    case (4)
+        do i=1,npoints
+           xrho = rho(i) * inv_rho_gf
+           xtemp = temp(i)
+           xye = ye(i)
+           xenr = eps(i) * inv_eps_gf
+           call nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,&
+                xdpderho,&
+                keytemp,keyerr(i),rf_precision)
+
+           if(keyerr(i).ne.0) then
+              anyerr = 1
+           endif
+
+           if(keytemp.eq.1) then
+              eps(i) = xenr * eps_gf
+           else
+              temp(i) = xtemp
+           endif
+
+           dpdrhoe(i) = xdpdrhoe * press_gf * inv_rho_gf
+           dpderho(i) = xdpderho * press_gf * inv_eps_gf
+        enddo
+   case DEFAULT
+      write(warnstring,*) "eoskey ",eoskey," not implemented!"
+      call CCTK_WARN(0,warnstring)
+   end select
+
+ end subroutine EOS_Omni_EOS_dpderho_dpdrhoe
diff --git a/src/nuc_eos/findtemp.F90 b/src/nuc_eos/findtemp.F90
index d208a5f..646994d 100644
--- a/src/nuc_eos/findtemp.F90
+++ b/src/nuc_eos/findtemp.F90
@@ -16,10 +16,10 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps)
   integer :: rl = 0
   integer itmax,i,keyerrt
   integer ii,jj,kk
-  
+
   keyerrt=0
 
-  tol=rfeps ! need to find energy to less than 1 in 10^-10
+  tol=rfeps ! need to find energy to less than 1 in 10^10
   itmax=20 ! use at most 20 iterations, then bomb
 
   lt=lt0
@@ -95,12 +95,14 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps)
           keyerrt=0
           goto 12
        else
+#if 0
           ! total failure
           write(*,*) "EOS: Did not converge in findtemp!"
           write(*,*) "rl,logrho,logtemp0,ye,lt,eps,eps0,abs(eps-eps0)/eps0"
           write(*,"(i4,i4,1P10E19.10)") i,rl,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0
           write(*,*) "Tried calling bisection... didn't help... :-/"
           write(*,*) "Bisection error: ",keyerrt
+#endif
        endif
     endif
     
diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90
index e84ac82..399d89b 100644
--- a/src/nuc_eos/nuc_eos.F90
+++ b/src/nuc_eos/nuc_eos.F90
@@ -348,6 +348,91 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,&
 
 end subroutine nuc_eos_press_eps
 
+subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,&
+     xdpderho,&
+     keytemp,keyerr,rfeps)
+
+  use eosmodule
+  implicit none
+
+  real*8, intent(in)    :: xrho,xye
+  real*8, intent(inout) :: xtemp,xenr
+  real*8, intent(in)    :: rfeps
+  real*8, intent(out)   :: xdpdrhoe,xdpderho
+  integer, intent(in)   :: keytemp
+  integer, intent(out)  :: keyerr
+
+  ! local variables
+  real*8 :: xcs2,xprs
+  real*8 :: lr,lt,y,xx,xeps,leps,xs
+  real*8 :: d1,d2,d3,ff(8)
+  integer :: keyerrt = 0
+
+  if(xrho.gt.eos_rhomax) then
+     stop "nuc_eos: rho > rhomax"
+  endif
+
+  if(xrho.lt.eos_rhomin*1.2d0) then
+     call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
+     return
+  endif
+
+  if(xye.gt.eos_yemax) then
+     stop "nuc_eos: ye > yemax"
+  endif
+
+  if(xye.lt.eos_yemin) then
+     keyerr = 102
+     write(6,*) "ye: ", xye
+     stop "nuc_eos ye < yemin"
+  endif
+
+  if(keytemp.eq.1) then
+     if(xtemp.gt.eos_tempmax) then
+        stop "nuc_eos: temp > tempmax"
+     endif
+     
+     if(xtemp.lt.eos_tempmin) then
+        call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
+        return
+     endif
+  endif
+
+  lr = log10(xrho)
+  lt = log10(xtemp)
+  y = xye
+  xeps = xenr + energy_shift
+  leps = log10(max(xeps,1.0d0))
+
+  keyerr = 0
+
+  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
+        keyerr = keyerrt
+        return
+     endif
+     xtemp = 10.0d0**lt
+
+  elseif(keytemp.eq.2) then
+     stop "eos_nuc_press does not support keytemp.eq.2"
+  endif
+
+  ! have temperature, proceed:
+  call findall_dpdr_dpde(lr,lt,y,ff)
+  xdpdrhoe = 10.0d0**ff(1)
+  xdpderho = 10.0d0**ff(2)
+
+
+  if((keytemp.eq.1).or.(keytemp.eq.2)) then
+     call findthis(lr,lt,y,xeps,alltables(:,:,:,2),d1,d2,d3)
+     xeps = 10.0d0**xeps - energy_shift
+     xenr = xeps
+  endif
+
+end subroutine nuc_eos_dpdr_dpde
+
 
 subroutine findthis(lr,lt,y,value,array,d1,d2,d3)
 
@@ -426,3 +511,22 @@ subroutine findall_press_eps(lr,lt,y,ff)
   ff(:) = ffx(:,1)
 
 end subroutine findall_press_eps
+
+subroutine findall_dpdr_dpde(lr,lt,y,ff)
+
+  use eosmodule
+  implicit none
+
+  real*8 ffx(8,1)
+  real*8 ff(8)
+  real*8 lr,lt,y
+  integer i
+  integer :: nvarsx = 2
+
+
+! Ewald's interpolator           
+  call intp3d_many(lr,lt,y,ffx,1,alltables(:,:,:,7:8), &
+       nrho,ntemp,nye,nvarsx,logrho,logtemp,ye)
+  ff(:) = ffx(:,1)
+
+end subroutine findall_dpdr_dpde