* EOS_Omni thorn:

    Initial implementation providing a polytrope and a gamma-law EOS for testing.



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

1 files changed, 333 insertions, 0 deletions

diff --git a/src/EOS_Omni_SingleVarCalls.F90 b/src/EOS_Omni_SingleVarCalls.F90
new file mode 100644
index 0000000..aa0f79d
--- /dev/null
+++ b/src/EOS_Omni_SingleVarCalls.F90
@@ -0,0 +1,333 @@
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+
+! eoskey:
+! 1 --- polytropic EOS
+! 2 --- gamma-law EOS
+! 3 --- hybrid EOS
+! 4 --- finite-T microphysical NSE EOS
+
+subroutine EOS_Omni_EOS_Press(eoskey,keytemp,npoints,&
+                              rho,eps,temp,ye,press,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)    :: rho(npoints),ye(npoints)
+  CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints)
+  CCTK_REAL, intent(out)   :: press(npoints)
+
+  ! local vars
+  integer        :: i
+  character(256) :: warnstring
+
+  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
+           press(i) = press_gf * poly_k_cgs * &
+                (rho(i)*inv_rho_gf)**poly_gamma
+        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
+           press(i) = (gl_gamma - 1.0d0) * rho(i) * eps(i)
+        enddo
+
+     case DEFAULT
+        write(warnstring,*) "eoskey ",eoskey," not implemented!"
+        call CCTK_WARN(0,warnstring)
+     end select
+
+end subroutine EOS_Omni_EOS_Press
+
+
+
+subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,npoints,&
+                              rho,eps,temp,ye,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)    :: rho(npoints),ye(npoints)
+  CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints)
+  CCTK_REAL, intent(out)   :: dpdrhoe(npoints)
+
+  ! local vars
+  integer        :: i
+  character(256) :: warnstring
+
+  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)
+        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)
+        enddo
+
+     case DEFAULT
+        write(warnstring,*) "eoskey ",eoskey," not implemented!"
+        call CCTK_WARN(0,warnstring)
+   end select
+     
+end subroutine EOS_Omni_EOS_DPressByDRho
+
+subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,npoints,&
+                              rho,eps,temp,ye,dpdepsrho,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)    :: rho(npoints),ye(npoints)
+  CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints)
+  CCTK_REAL, intent(out)   :: dpdepsrho(npoints)
+
+  ! local vars
+  integer        :: i
+  character(256) :: warnstring
+
+  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
+           dpdepsrho(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
+           dpdepsrho(i) = (gl_gamma - 1.0d0) * &
+                rho(i)
+        enddo
+
+     case DEFAULT
+        write(warnstring,*) "eoskey ",eoskey," not implemented!"
+        call CCTK_WARN(0,warnstring)
+   end select
+     
+end subroutine EOS_Omni_EOS_DPressByDEps
+
+subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,npoints,&
+                            rho,eps,temp,ye,cs2,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)    :: rho(npoints),ye(npoints)
+  CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints)
+  CCTK_REAL, intent(out)   :: cs2(npoints)
+
+  ! local vars
+  integer        :: i
+  character(256) :: warnstring
+  real*8         :: xpress,xdpdrhoe,xdpderho
+
+  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
+           xpress = press_gf*poly_k_cgs * &
+                (rho(i)*inv_rho_gf)**(poly_gamma)
+           cs2(i) = poly_gamma * xpress / rho(i) / &
+                (1 + eps(i) + xpress/rho(i))
+        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
+           xpress = (gl_gamma-1.0d0)*rho(i)*eps(i)
+           xdpdrhoe = (gl_gamma-1.0d0)*eps(i)
+           xdpderho = (gl_gamma-1.0d0)*rho(i)
+           cs2(i) = (xdpdrhoe + xpress * xdpderho / (rho(i)**2)) / &
+                (1.0d0 + eps(i) + xpress/rho(i))
+        enddo
+
+     case DEFAULT
+        write(warnstring,*) "eoskey ",eoskey," not implemented!"
+        call CCTK_WARN(0,warnstring)
+   end select
+     
+end subroutine EOS_Omni_EOS_cs2
+
+
+subroutine EOS_Omni_EOS_eps_from_press(eoskey,keytemp,npoints,&
+                              rho,eps,temp,ye,press,xeps,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)    :: rho(npoints),ye(npoints),press(npoints)
+  CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints)
+  CCTK_REAL, intent(out)   :: xeps(npoints)
+
+  ! local vars
+  integer        :: i
+  character(256) :: warnstring
+
+  if(keytemp.eq.1) then
+     anyerr = 1
+     keyerr(:) = -1
+  else
+     anyerr    = 0
+     keyerr(:) = 0
+  endif
+
+  select case (eoskey)
+     case (1)
+        ! polytropic EOS
+        do i=1,npoints
+           xeps(i) = press(i) / (poly_gamma - 1.0d0) / rho(i)
+           eps(i)  = xeps(i)
+        enddo
+     case (2)
+        ! gamma-law EOS
+        do i=1,npoints
+           xeps(i) = press(i) / (gl_gamma - 1.0d0) / rho(i)
+           eps(i)  = xeps(i)
+        enddo
+     case DEFAULT
+        write(warnstring,*) "eoskey ",eoskey," not implemented!"
+        call CCTK_WARN(0,warnstring)
+     end select
+
+end subroutine EOS_Omni_EOS_eps_from_press
+
+
+subroutine EOS_Omni_EOS_RestMassDensityFromEpsPress(eoskey,keytemp,npoints,&
+                              rho,eps,temp,ye,press,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)    :: ye(npoints),press(npoints),eps(npoints)
+  CCTK_REAL, intent(inout) :: rho(npoints),temp(npoints)
+
+
+  ! local vars
+  integer        :: i
+  character(256) :: warnstring
+
+  if(keytemp.eq.1) then
+     anyerr = 1
+     keyerr(:) = -1
+  else
+     anyerr    = 0
+     keyerr(:) = 0
+  endif
+
+  select case (eoskey)
+     case (1)
+        ! polytropic EOS
+        do i=1,npoints
+           rho(i) = press(i) / ((poly_gamma - 1.0d0)*eps(i))
+        enddo
+     case (2)
+        ! gamma-law EOS
+        do i=1,npoints
+           rho(i) = press(i) / ((gl_gamma - 1.0d0)*eps(i))
+        enddo
+     case DEFAULT
+        write(warnstring,*) "eoskey ",eoskey," not implemented!"
+        call CCTK_WARN(0,warnstring)
+     end select
+
+end subroutine EOS_Omni_EOS_RestMassDensityFromEpsPress