* add support for cold tabulated EOS (P=P(rho)) that

  can be supplemented with a thermal gamma law



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

4 files changed, 296 insertions, 5 deletions

diff --git a/src/EOS_Omni_Handles.c b/src/EOS_Omni_Handles.c
index b744cac..1ecd349 100644
--- a/src/EOS_Omni_Handles.c
+++ b/src/EOS_Omni_Handles.c
@@ -12,6 +12,8 @@ CCTK_INT EOS_Omni_GetHandle_(CCTK_STRING name)
         return 3;
     if (CCTK_EQUALS(name, "nuc_eos"))
         return 4;
+    if (CCTK_EQUALS(name, "cold_tabulated"))
+        return 5;
     return 0;
 }
 
diff --git a/src/EOS_Omni_Module.F90 b/src/EOS_Omni_Module.F90
index e6ecc14..04dd10d 100644
--- a/src/EOS_Omni_Module.F90
+++ b/src/EOS_Omni_Module.F90
@@ -25,7 +25,21 @@ module EOS_Omni_Module
   real*8,parameter :: cliteinv2 = 1.11265005605362d-21
 
   ! These values are set by EOS_Omni_Startup
-
   real*8 :: hybrid_k2 = 0.0d0
 
+  ! stuff for the cold, tabulated EOS with a gamma law
+  ! set by the reader routine
+  integer :: coldeos_nrho = 0
+  real*8 :: coldeos_gammath = 0.0d0
+  real*8 :: coldeos_rhomin = 0.0d0
+  real*8 :: coldeos_rhomax = 0.0d0
+  real*8 :: coldeos_kappa = 0.0d0
+  real*8 :: coldeos_thfac = 1.0d0
+  real*8 :: coldeos_dlrho = 1.0d0
+  real*8 :: coldeos_dlrhoi = 1.0d0
+  real*8, allocatable :: coldeos_logrho(:)
+  real*8, allocatable :: coldeos_eps(:)
+  real*8, allocatable :: coldeos_gamma(:)
+  real*8, allocatable :: coldeos_cs2(:)
+
 end module EOS_Omni_Module
diff --git a/src/EOS_Omni_SingleVarCalls.F90 b/src/EOS_Omni_SingleVarCalls.F90
index bbff17c..350e622 100644
--- a/src/EOS_Omni_SingleVarCalls.F90
+++ b/src/EOS_Omni_SingleVarCalls.F90
@@ -34,6 +34,11 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
   real*8           :: xrho,xye,xtemp,xenr,xent
   real*8           :: xprs,xmunu,xcs2
   real*8           :: xdedt,xdpderho,xdpdrhoe
+  ! temporary vars for coldeos + gamma law
+  integer :: ir
+  real*8 :: press_cold, press_th
+  real*8 :: eps_cold, eps_th
+  real*8 :: gamma
 
   anyerr    = 0
   keyerr(:) = 0
@@ -89,7 +94,7 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
         enddo
 
      case (4)
-
+        ! nuc eos
         do i=1,npoints
 
            xrho = rho(i) * inv_rho_gf
@@ -113,6 +118,41 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
 
         enddo
 
+     case (5)
+        ! cold tabular EOS with gamma law
+        do i=1,npoints
+           if(rho(i).lt.coldeos_rhomin) then
+              keyerr(i) = 104
+              anyerr = 1
+           else if(rho(i).gt.coldeos_rhomax) then
+              keyerr(i) = 103
+              anyerr = 1
+           else
+              xrho = log10(rho(i))
+              ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
+           endif
+           ir = max(2, min(ir,coldeos_nrho))
+
+           gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / &
+                (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1)
+
+           eps_cold = (coldeos_eps(ir) - coldeos_eps(ir-1)) / &
+                    (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                    (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1)
+
+           if(keytemp.eq.0) then
+              eps_th = max(0.0d0,eps(i) - eps_cold)
+           else if (keytemp.eq.1) then
+              eps_th = 0.0d0
+              eps(i) = eps_cold
+           endif
+
+           press_cold = coldeos_kappa * rho(i)**gamma
+           press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th
+           press(i) = press_cold !+ press_th
+        enddo
+
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)
@@ -147,6 +187,11 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
   real*8           :: xrho,xye,xtemp,xenr,xent
   real*8           :: xprs,xmunu,xcs2
   real*8           :: xdedt,xdpderho,xdpdrhoe
+  ! temporary vars for cold tabulated EOS + gamma law
+  integer :: ir
+  real*8 :: press_cold, press_th
+  real*8 :: eps_cold, eps_th
+  real*8 :: gamma, cs2, h
 
   anyerr    = 0
   keyerr(:) = 0
@@ -205,8 +250,6 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
 
            dpdrhoe(i) = hybrid_dp_poly + max(0.0d0,hybrid_dp_th1 + hybrid_dp_th2)
 
-           
-
         enddo
 
      case (4)
@@ -231,6 +274,54 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
 
            dpdrhoe(i) = xdpdrhoe * press_gf * inv_rho_gf
         enddo
+        
+     case (5)
+        ! with the cold eos we have to assume P = P(rho), so
+        ! by definition dPdrho is at constant internal energy
+        ! and entropy (the latter, because T=0)
+        do i=1,npoints
+           if(rho(i).lt.coldeos_rhomin) then
+              keyerr(i) = 104
+              anyerr = 1
+           else if(rho(i).gt.coldeos_rhomax) then
+              keyerr(i) = 103
+              anyerr = 1
+           else
+              xrho = log10(rho(i))
+              ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
+           endif
+           ir = max(2, min(ir,coldeos_nrho))
+
+           gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / &
+                (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1)
+           
+           ! this is the cold speed of sound squared
+           cs2 = (coldeos_cs2(ir) - coldeos_cs2(ir-1)) / &
+                (coldeos_cs2(ir) - coldeos_cs2(ir-1)) * &
+                (xrho - coldeos_logrho(ir-1)) + coldeos_cs2(ir-1)
+
+           eps_cold = (coldeos_eps(ir) - coldeos_eps(ir-1)) / &
+                    (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                    (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1)
+
+           if(keytemp.eq.0) then
+              eps_th = max(0.0d0,eps(i) - eps_cold)
+           else if (keytemp.eq.1) then
+              eps_th = 0.0d0
+              eps(i) = eps_cold
+           endif
+
+           press_cold = coldeos_kappa * rho(i)**gamma
+           press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th
+           
+           ! this is the cold enthalpy, because it belongs to the
+           ! cold speed of sound squared
+           h = 1.0d0 + eps_cold + press_cold / rho(i)
+           dpdrhoe(i) = cs2*h + press_th / rho(i)
+
+        enddo
+
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)
@@ -260,6 +351,9 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,&
   real*8           :: xrho,xye,xtemp,xenr,xent
   real*8           :: xprs,xmunu,xcs2
   real*8           :: xdedt,xdpderho,xdpdrhoe
+  ! temporary vars for coldeos + gamma law                                                              
+  integer :: ir
+  real*8 :: eps_cold, eps_th
 
   anyerr    = 0
   keyerr(:) = 0
@@ -322,6 +416,36 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,&
 
         enddo
 
+     case (5)
+        ! with the cold eos we have to assume P = P(rho), so
+        ! only the gamma law has non-zero dPdeps
+        do i=1,npoints
+           if(rho(i).lt.coldeos_rhomin) then
+              keyerr(i) = 104
+              anyerr = 1
+           else if(rho(i).gt.coldeos_rhomax) then
+              keyerr(i) = 103
+              anyerr = 1
+           else
+              xrho = log10(rho(i))
+              ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
+           endif
+           ir = max(2, min(ir,coldeos_nrho))
+
+           eps_cold = (coldeos_eps(ir) - coldeos_eps(ir-1)) / &
+                    (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                    (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1)
+
+           if(keytemp.eq.0) then
+              eps_th = max(0.0d0,eps(i) - eps_cold)
+           else if (keytemp.eq.1) then
+              eps_th = 0.0d0
+              eps(i) = eps_cold
+           endif
+
+           dpdepsrho(i) = coldeos_thfac * (coldeos_gammath - 1.0d0)*rho(i)
+
+        enddo
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)
@@ -355,6 +479,11 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
   real*8           :: xrho,xye,xtemp,xenr,xent
   real*8           :: xprs,xmunu,xcs2
   real*8           :: xdedt
+  ! temporary vars for cold tabulated EOS + gamma law                                                   
+  integer :: ir
+  real*8 :: press_cold, press_th
+  real*8 :: eps_cold, eps_th
+  real*8 :: gamma, cs2_cold, cs2_th, h, h_cold
 
   anyerr    = 0
   keyerr(:) = 0
@@ -443,8 +572,56 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
                 (1.0d0 + eps(i) + xprs * press_gf / rho(i))
 
         enddo
+     case (5)
+        ! with the cold eos we have to assume P = P(rho), so                                            
+        ! by definition dPdrho is at constant internal energy                                           
+        ! and entropy (the latter, because T=0)                                                         
+        do i=1,npoints
+           if(rho(i).lt.coldeos_rhomin) then
+              keyerr(i) = 104
+              anyerr = 1
+           else if(rho(i).gt.coldeos_rhomax) then
+              keyerr(i) = 103
+              anyerr = 1
+                 else
+              xrho = log10(rho(i))
+              ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
+           endif
+           ir = max(2, min(ir,coldeos_nrho))
+
+           gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / &
+                (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1)
+
+           ! this is the cold speed of sound squared                                                    
+           cs2_cold = (coldeos_cs2(ir) - coldeos_cs2(ir-1)) / &
+                (coldeos_cs2(ir) - coldeos_cs2(ir-1)) * &
+                (xrho - coldeos_logrho(ir-1)) + coldeos_cs2(ir-1)
+
+           eps_cold = (coldeos_eps(ir) - coldeos_eps(ir-1)) / &
+                    (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                    (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1)
+
+           if(keytemp.eq.0) then
+              eps_th = max(0.0d0,eps(i) - eps_cold)
+           else if (keytemp.eq.1) then
+              eps_th = 0.0d0
+              eps(i) = eps_cold
+           endif
 
+           press_cold = coldeos_kappa * rho(i)**gamma
+           press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th
 
+           xdpdrhoe = coldeos_thfac*(coldeos_gammath - 1.0d0)*eps_th
+           xdpderho = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)
+           cs2_th = (xdpdrhoe + press_th * xdpderho / (rho(i)**2))
+
+           h = 1.0d0 + eps(i) + (press_cold+press_th) / rho(i)
+           h_cold = 1.0d0 + eps_cold + press_cold / rho(i)
+           
+           cs2(i) = (cs2_cold * h_cold + cs2_th) / h
+        enddo
+           
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)
@@ -541,6 +718,9 @@ subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt(eoskey,&
   real*8           :: xrho,xye,xtemp,xenr,xent
   real*8           :: xprs,xmunu,xcs2
   real*8           :: xdedt,xdpderho,xdpdrhoe
+  ! temporary vars for cold tabulated EOS
+  integer :: ir
+  real*8  :: gamma
   
   ! helper vars
   logical :: cont
@@ -632,6 +812,101 @@ subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt(eoskey,&
               endif
            enddo
         enddo
+        
+     case (5)
+        ! cold tabulate EOS
+        ! deal only with case in which thermal pressure is zero
+        if(keytemp.ne.1) then
+           call CCTK_WARN(0,"finding rho(press) for tabulated cold EOS only possible if keytemp=1")
+        endif
+
+        do i=1,npoints
+           fac = 1.0d0
+           counter = 0
+           cont = .true.
+           xprs = press(i) 
+           press_guess = xprs
+           rho_guess = rho(i) 
+           xprs = press(i) 
+           do while(cont)
+              counter = counter + 1
+              rho_guess2 = rho_guess * 1.0001d0
+
+              ! press 2
+              if(rho_guess2.lt.coldeos_rhomin) then
+                 keyerr(i) = 104
+                 anyerr = 1
+              else if(rho_guess2.gt.coldeos_rhomax) then
+                 keyerr(i) = 103
+                 anyerr = 1
+              else
+                 xrho = log10(rho_guess2)
+                 ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
+              endif
+              ir = max(2, min(ir,coldeos_nrho))
+           
+              gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / &
+                (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1)
+
+              press_guess2 = coldeos_kappa * rho_guess2**gamma
+              ! end press 2
+              ! press 1
+              if(rho_guess.lt.coldeos_rhomin) then
+                 keyerr(i) = 104
+                 anyerr = 1
+              else if(rho_guess.gt.coldeos_rhomax) then
+                 keyerr(i) = 103
+                 anyerr = 1
+              else
+                 xrho = log10(rho_guess)
+                 ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
+              endif
+              ir = max(2, min(ir,coldeos_nrho))
+           
+              gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / &
+                (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1)
+
+              press_guess = coldeos_kappa * rho_guess**gamma
+              ! end press 1
+
+              ! derivative
+              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,rho_guess2
+                 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
+                 eps(i) = (coldeos_eps(ir) - coldeos_eps(ir-1)) / &
+                    (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                    (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1)
+              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.0d-15.or.counter.gt.100000) then
+                 keyerr(i) = 473
+                 anyerr = 1
+                 return
+              endif
+           enddo
+
+        enddo
 
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
diff --git a/src/make.code.defn b/src/make.code.defn
index c071792..f23abc5 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -3,7 +3,7 @@
 # Source files in this directory
 SRCS = EOS_Omni_Module.F90 EOS_Omni_Startup.F90 EOS_Omni_SingleVarCalls.F90 \
        EOS_Omni_SingleVarCalls_harm.F90 EOS_Omni_Handles.c \
-       EOS_Omni_MultiVarCalls.F90
+       EOS_Omni_MultiVarCalls.F90 EOS_Omni_ColdEOSReadTable.F90
 
 # Subdirectories containing source files
 SUBDIRS = nuc_eos