diff options
Diffstat (limited to 'src/EOS_Hybrid.F')
| -rw-r--r-- | src/EOS_Hybrid.F | 144 |
1 files changed, 144 insertions, 0 deletions
diff --git a/src/EOS_Hybrid.F b/src/EOS_Hybrid.F new file mode 100644 index 0000000..c74af75 --- /dev/null +++ b/src/EOS_Hybrid.F @@ -0,0 +1,144 @@ + /*@@ + @file EOS_Hybrid.F + @date Wed Mar 20 14:56:35 2002 + @author Ian Hawke + @desc + Routines to calculate the EOS used by + Dimmelmeier Novak Font Ibanez Mueller PRD71 064023 (2005) + in supernova core collapse simulations. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" + + CCTK_REAL function EOS_Hybrid_Pressure(rho, eps) + + USE EOS_Polytrope_Scalars + USE EOS_Hybrid_Scalars + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_REAL rho, eps, local_eos_gamma, p_poly, p_th, local_eos_k_cgs, zero + + zero = 0.d0 + + if (rho > rho_nuc) then + local_eos_gamma = eos_gamma_supernuclear + local_eos_k_cgs = eos_k_supernuclear_cgs + else + local_eos_gamma = eos_gamma + local_eos_k_cgs = eos_k_cgs + end if + + p_poly = p_geom_factor * local_eos_k_cgs * + . (rho * rho_geom_factor_inv)**local_eos_gamma + + p_th = - p_geom_factor * local_eos_k_cgs * (eos_gamma_th - 1.d0) / + . (local_eos_gamma - 1.d0) * (rho * rho_geom_factor_inv)**local_eos_gamma + + . (eos_gamma_th - 1.d0) * rho * eps - + . (eos_gamma_th - 1.d0) * (local_eos_gamma - eos_gamma) / + . (eos_gamma - 1.d0) / (eos_gamma_supernuclear - 1.d0) * + . p_geom_factor * eos_k_cgs * rho_geom_factor_inv**eos_gamma * + . rho_nuc**(eos_gamma - 1.d0) * rho + + p_th = max(zero, p_th) + + EOS_Hybrid_Pressure = p_poly + p_th + + end function EOS_Hybrid_Pressure + +c The specific internal energy isn''t correct yet + + CCTK_REAL function EOS_Hybrid_SpecificIE(rho, press) + + USE EOS_Polytrope_Scalars + USE EOS_Hybrid_Scalars + + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_REAL rho, press + +c$$$ if (rho > rho_nuc) then +c$$$ EOS_Hybrid_SpecificIE = eos_k_cgs / (eos_gamma_2 - 1.d0) * +c$$$ . (rho / rho_nuc) ** eos_gamma_2 * rho_nuc ** eos_gamma_1 / +c$$$ . rho + +c$$$ . eos_k_cgs * (eos_gamma_2 - eos_gamma_1) / +c$$$ . (eos_gamma_2 - 1.d0) / (eos_gamma_1 - 1.d0) * +c$$$ . rho_nuc ** (eos_gamma_1 - 1.d0) +c$$$ else +c$$$ EOS_Hybrid_SpecificIE = eos_k_cgs * (rho ** (eos_gamma_1 - 1.d0)) / +c$$$ . (eos_gamma_1 - 1.d0) +c$$$ end if + + EOS_Hybrid_SpecificIE = 2.34567890d0 + + end function EOS_Hybrid_SpecificIE + +c The rest mass density isn''t correct yet + + CCTK_REAL function EOS_Hybrid_RestMassDens(eps, press) + + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_REAL eps, press + + EOS_Hybrid_RestMassDens = 1.23456789d0 + + end function EOS_Hybrid_RestMassDens + + CCTK_REAL function EOS_Hybrid_DPressByDRho(rho, eps) + + USE EOS_Polytrope_Scalars + USE EOS_Hybrid_Scalars + + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_REAL rho, eps, local_eos_gamma, local_eos_k_cgs, d_p_poly, d_p_th_1, + . d_p_th_2, zero + + zero = 0.d0 + + if (rho > rho_nuc) then + local_eos_gamma = eos_gamma_supernuclear + local_eos_k_cgs = eos_k_supernuclear_cgs + else + local_eos_gamma = eos_gamma + local_eos_k_cgs = eos_k_cgs + end if + + d_p_poly = local_eos_gamma * p_geom_factor * local_eos_k_cgs * + . rho**(local_eos_gamma - 1.d0) * rho_geom_factor_inv**local_eos_gamma + + d_p_th_1 = - local_eos_gamma * p_geom_factor * local_eos_k_cgs * + . (eos_gamma_th - 1.d0) / (local_eos_gamma - 1.d0) * + . rho**(local_eos_gamma - 1.d0) * rho_geom_factor_inv**local_eos_gamma + + d_p_th_2 = (eos_gamma_th - 1.d0) * eps + . - (eos_gamma_th - 1.d0) * (local_eos_gamma - eos_gamma) / + . (eos_gamma - 1.d0) / (eos_gamma_supernuclear - 1.d0) * + . p_geom_factor * eos_k_cgs * rho_geom_factor_inv**eos_gamma * + . rho_nuc**(eos_gamma - 1.d0) + +! d_p_th_1 = max(d_p_th_1, zero) +! d_p_th_2 = max(d_p_th_2, zero) + + EOS_Hybrid_DPressByDRho = d_p_poly + d_p_th_1 + d_p_th_2 + + end function EOS_Hybrid_DPressByDRho + + CCTK_REAL function EOS_Hybrid_DPressByDEps(rho, eps) + + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_REAL rho, eps + + EOS_Hybrid_DPressByDEps = (eos_gamma_th - 1.d0) * rho + + end function EOS_Hybrid_DPressByDEps |