[Commits] [svn:einsteintoolkit] EOS_Omni/trunk/ (Rev. 103)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Wed Mar 12 22:02:02 CDT 2014
User: rhaas
Date: 2014/03/12 10:02 PM
Modified:
/trunk/
interface.ccl
/trunk/src/
EOS_Omni_MultiVarCalls.F90
Log:
make new EOS call that returns press and cs2 [and eps if needed]
From: Christian D. Ott <cott at tapir.caltech.edu>
File Changes:
Directory: /trunk/src/
======================
File [modified]: EOS_Omni_MultiVarCalls.F90
Delta lines: +177 -8
===================================================================
--- trunk/src/EOS_Omni_MultiVarCalls.F90 2014-03-13 03:01:59 UTC (rev 102)
+++ trunk/src/EOS_Omni_MultiVarCalls.F90 2014-03-13 03:02:01 UTC (rev 103)
@@ -101,7 +101,6 @@
CCTK_REAL, intent(out) :: muhat(npoints)
! local vars
- integer :: i
character(256) :: warnstring
if(eoskey.ne.4) then
@@ -157,10 +156,6 @@
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
@@ -264,9 +259,7 @@
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
+ real*8 :: xdpderho,xdpdrhoe
anyerr = 0
@@ -342,3 +335,179 @@
end select
end subroutine EOS_Omni_EOS_DEpsByDRho_DEpsByDPress
+
+
+subroutine EOS_Omni_EOS_Press_cs2(eoskey,keytemp,rf_precision,npoints,&
+ rho,eps,temp,ye,press,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) :: rf_precision
+ CCTK_REAL, intent(in) :: rho(npoints),ye(npoints)
+ CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints)
+ CCTK_REAL, intent(out) :: press(npoints)
+ CCTK_REAL, intent(out) :: cs2(npoints)
+
+ ! local vars
+ integer :: i
+ character(256) :: warnstring
+ real*8 :: hybrid_local_gamma, hybrid_local_k_cgs, &
+ hybrid_p_poly, hybrid_p_th
+ real*8,parameter :: zero = 0.0d0
+ ! temporary vars for nuc_eos
+ real*8 :: xrho
+ real*8 :: 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, cs2_cold, cs2_th, h, h_cold
+
+ 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
+ cs2(i) = poly_gamma * press(i) / rho(i) / &
+ (1 + eps(i) + press(i)/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
+ press(i) = (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 + press(i) * xdpderho / (rho(i)**2)) / &
+ (1.0d0 + eps(i) + press(i)/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_p_poly = press_gf * hybrid_local_k_cgs * &
+ (rho(i) * inv_rho_gf)**hybrid_local_gamma
+ hybrid_p_th = - press_gf * hybrid_local_k_cgs * (hybrid_gamma_th - 1.d0) / &
+ (hybrid_local_gamma - 1.0d0) * (rho(i) * inv_rho_gf)**hybrid_local_gamma + &
+ (hybrid_gamma_th - 1.0d0) * rho(i) * 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) * rho(i)
+ hybrid_p_th = max(zero, hybrid_p_th)
+ press(i) = hybrid_p_poly + hybrid_p_th
+
+ cs2(i) = (hybrid_local_gamma * hybrid_p_poly + hybrid_gamma_th * hybrid_p_th) / &
+ rho(i) / (1.0d0 + eps(i) + press(i)/rho(i))
+ enddo
+
+ case (4)
+ ! nuc eos
+ if(keytemp.eq.1) then
+ call nuc_eos_m_kt1_press_eps_cs2(npoints,rho,temp,ye,&
+ eps,press,cs2,keyerr,anyerr)
+ else if(keytemp.eq.0) then
+ call nuc_eos_m_kt0_press_cs2(npoints,rho,temp,ye,eps,press,&
+ cs2,rf_precision,keyerr,anyerr)
+ else
+ call CCTK_ERROR("This keytemp is not suppported!")
+ STOP
+ endif
+ 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
+ press(i) = coldeos_low_kappa * rho(i)**coldeos_low_gamma
+ eps(i) = press(i) / (coldeos_low_gamma - 1.0) / rho(i)
+ cs2(i) = coldeos_low_kappa * coldeos_low_gamma * &
+ rho(i)**(coldeos_low_gamma-1.0d0) / &
+ (1.0d0 + eps(i) + press(i))
+ cycle
+ 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
+ press(i) = press_cold + press_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 (6)
+ ! barotropic tabular EOS with gamma law
+ write(warnstring,*) "eoskey ",eoskey," for barotropic EOS not implemented!"
+ call CCTK_ERROR(warnstring)
+ STOP
+ case DEFAULT
+ write(warnstring,*) "eoskey ",eoskey," not implemented!"
+ call CCTK_ERROR(warnstring)
+ STOP
+ end select
+
+end subroutine EOS_Omni_EOS_Press_cs2
+
Directory: /trunk/
==================
File [modified]: interface.ccl
Delta lines: +15 -0
===================================================================
--- trunk/interface.ccl 2014-03-13 03:01:59 UTC (rev 102)
+++ trunk/interface.ccl 2014-03-13 03:02:01 UTC (rev 103)
@@ -22,6 +22,21 @@
PROVIDES FUNCTION EOS_Omni_press WITH EOS_Omni_EOS_Press LANGUAGE Fortran
+void FUNCTION EOS_Omni_press_cs2(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 press, \
+ CCTK_REAL OUT ARRAY cs2, \
+ CCTK_INT OUT ARRAY keyerr, \
+ CCTK_INT OUT anyerr)
+
+PROVIDES FUNCTION EOS_Omni_press_cs2 WITH EOS_Omni_EOS_Press_cs2 LANGUAGE Fortran
+
void FUNCTION EOS_Omni_pressOMP(CCTK_INT IN eoskey, \
CCTK_INT IN havetemp, \
CCTK_REAL IN rf_precision, \
More information about the Commits
mailing list