* improve solution for density based on pressure, temperature, Y_e.

  The new implementation is robust, but slow and should only be used
  at initial data (which is generally the case)



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

1 files changed, 61 insertions, 21 deletions

diff --git a/src/EOS_Omni_SingleVarCalls.F90 b/src/EOS_Omni_SingleVarCalls.F90
index 783f2c1..8e3b92c 100644
--- a/src/EOS_Omni_SingleVarCalls.F90
+++ b/src/EOS_Omni_SingleVarCalls.F90
@@ -515,6 +515,11 @@ subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt(eoskey,&
                            ikeytemp,rf_precision,&
                            npoints,rho,eps,temp,entropy,ye,press,keyerr,anyerr)
 
+  ! This routine returns the density and spec. internal energy based on inputs of
+  ! pressure, temperature, and Y_e.
+  ! The current implementation is robust, but very slow; it should be used only
+  ! for initial data where speed does not matter.
+
   use EOS_Omni_Module
   implicit none
   DECLARE_CCTK_PARAMETERS
@@ -537,6 +542,12 @@ subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt(eoskey,&
   real*8           :: xprs,xmunu,xcs2
   real*8           :: xdedt,xdpderho,xdpdrhoe
   
+  ! helper vars
+  logical :: cont
+  integer :: counter
+  real*8  :: rho_guess, rho_guess2
+  real*8  :: press_guess,press_guess2, mydpdrho
+  real*8  :: fac
   keytemp = ikeytemp
 
   anyerr    = 0
@@ -563,36 +574,65 @@ subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt(eoskey,&
         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
+           keytemp = 1
         else
            call CCTK_WARN(0,"This function does not work yet when coming in with this keytemp")
         endif
 
+        keytemp = 1
         do i=1,npoints
-           ! initial guess
+           fac = 1.0d0
+           counter = 0
+           cont = .true.
+           xprs = press(i) * inv_press_gf
+           press_guess = xprs
+           rho_guess = rho(i) * inv_rho_gf
            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
-
+           do while(cont)
+              counter = counter + 1
+              rho_guess2 = rho_guess * 1.0001d0
+              call nuc_eos_short(rho_guess2,xtemp,xye,xenr,press_guess2,& 
+                   xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,&
+                   keytemp,keyerr(i),rf_precision)
+              call nuc_eos_short(rho_guess,xtemp,xye,xenr,press_guess,& 
+                   xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,&
+                   keytemp,keyerr(i),rf_precision)
+
+              mydpdrho = (press_guess2-press_guess)/(rho_guess2-rho_guess)
+              if (mydpdrho.lt.0.0d0) then
+                 !$OMP CRITICAL
+                 write(warnstring,"(A25,1P10E15.6)") "Issue with table, dpdrho.lt.0",&
+                      rho_guess,xtemp,xye
+                 call CCTK_WARN(0,warnstring)
+                 !$OMP END CRITICAL
+              endif
+              
+              if (abs(1.0d0-press_guess/xprs).lt.rf_precision) then
+                 cont = .false.
+                 rho(i) = rho_guess*rho_gf
+                 eps(i) = xenr*eps_gf
+              else
+                 if (fac*(xprs-press_guess)/mydpdrho/rho_guess.gt.0.1d0) then
+                    rho_guess = 0.99d0*rho_guess
+                 else
+                    rho_guess = rho_guess + fac*(xprs-press_guess)/mydpdrho
+                 endif
+              endif
+
+              if (counter.gt.100) then
+                 fac = 0.01d0
+              endif
+
+              if (rho_guess.lt.1.0d3.or.counter.gt.100000) then
+                 keyerr(i) = 473
+                 anyerr = 1
+                 return
+              endif
+           enddo
         enddo
+
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)