diff options
Diffstat (limited to 'src/EOS_GeneralHybrid.F90')
| -rwxr-xr-x | src/EOS_GeneralHybrid.F90 | 228 |
1 files changed, 228 insertions, 0 deletions
diff --git a/src/EOS_GeneralHybrid.F90 b/src/EOS_GeneralHybrid.F90 new file mode 100755 index 0000000..e842600 --- /dev/null +++ b/src/EOS_GeneralHybrid.F90 @@ -0,0 +1,228 @@ + /*@@ + @file EOS_GeneralHybrid.F90 + @date + @author Ian Hawke, Christian D. Ott + @desc + Routines to calculate the EOS used by Dimmelmeier et al. + in supernova core collapse simulations. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" + +subroutine EOS_GeneralHybrid_Pressure(nelems, rho, eps, press) + + USE EOS_GP_Scalars + USE EOS_GeneralHybrid_Scalars + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: nelems + CCTK_REAL, dimension(nelems), intent(in) :: rho + CCTK_REAL, dimension(nelems), intent(in) :: eps + CCTK_REAL, dimension(nelems), intent(out) :: press + + CCTK_REAL, allocatable, dimension(:) :: p_poly + CCTK_REAL, allocatable, dimension(:) :: p_th + CCTK_REAL, allocatable, dimension(:) :: local_eos_gamma + CCTK_REAL, allocatable, dimension(:) :: local_eos_k_cgs + + CCTK_REAL zero + + allocate(p_poly(1:nelems)) + allocate(p_th(1:nelems)) + allocate(local_eos_gamma(1:nelems)) + allocate(local_eos_k_cgs(1:nelems)) + + zero = 0.d0 + + where (rho > rho_nuc) + local_eos_gamma = eos_gamma_supernuclear + local_eos_k_cgs = eos_k_supernuclear_cgs + elsewhere + local_eos_gamma = eos_gamma + local_eos_k_cgs = eos_k_cgs + end where + + 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) + + press = p_poly + p_th + + deallocate(p_poly) + deallocate(p_th) + deallocate(local_eos_gamma) + deallocate(local_eos_k_cgs) + + +end subroutine EOS_GeneralHybrid_Pressure + +subroutine EOS_GeneralHybrid_DPDRho(nelems, rho, eps, dpdrho) + + USE EOS_GP_Scalars + USE EOS_GeneralHybrid_Scalars + + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: nelems + CCTK_REAL, dimension(nelems), intent(in) :: rho + CCTK_REAL, dimension(nelems), intent(in) :: eps + CCTK_REAL, dimension(nelems), intent(out) :: dpdrho + + CCTK_REAL, allocatable,dimension(:) :: d_p_poly + CCTK_REAL, allocatable,dimension(:) :: d_p_th_1 + CCTK_REAL, allocatable,dimension(:) :: d_p_th_2 + CCTK_REAL, allocatable,dimension(:) :: local_eos_gamma + CCTK_REAL, allocatable,dimension(:) :: local_eos_k_cgs + + CCTK_REAL zero + + allocate(d_p_poly(1:nelems)) + allocate(d_p_th_1(1:nelems)) + allocate(d_p_th_2(1:nelems)) + allocate(local_eos_gamma(1:nelems)) + allocate(local_eos_k_cgs(1:nelems)) + + + zero = 0.d0 + + where (rho > rho_nuc) + local_eos_gamma = eos_gamma_supernuclear + local_eos_k_cgs = eos_k_supernuclear_cgs + elsewhere + local_eos_gamma = eos_gamma + local_eos_k_cgs = eos_k_cgs + end where + + 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) + + dpdrho = d_p_poly + d_p_th_1 + d_p_th_2 + + deallocate(d_p_poly) + deallocate(d_p_th_1) + deallocate(d_p_th_2) + deallocate(local_eos_gamma) + deallocate(local_eos_k_cgs) + + +end subroutine EOS_GeneralHybrid_DPDRho + + +subroutine EOS_GeneralHybrid_cs2(nelems, rho, eps, cs2) + +! This is an ugly and quick hack. +! Uses way to many operations and memory. + + USE EOS_GP_Scalars + USE EOS_GeneralHybrid_Scalars + + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: nelems + CCTK_REAL, dimension(nelems), intent(in) :: rho + CCTK_REAL, dimension(nelems), intent(in) :: eps + CCTK_REAL, dimension(nelems), intent(out) :: cs2 + CCTK_REAL, allocatable,dimension(:) :: p_poly + CCTK_REAL, allocatable,dimension(:) :: p_th + CCTK_REAL, allocatable,dimension(:) :: pressure + + CCTK_REAL, allocatable,dimension(:) :: local_eos_gamma + CCTK_REAL, allocatable,dimension(:) :: local_eos_k_cgs + + CCTK_REAL zero + + allocate(p_poly(1:nelems)) + allocate(p_th(1:nelems)) + allocate(pressure(1:nelems)) + allocate(local_eos_gamma(1:nelems)) + allocate(local_eos_k_cgs(1:nelems)) + + + zero = 0.d0 + + where (rho > rho_nuc) + local_eos_gamma = eos_gamma_supernuclear + local_eos_k_cgs = eos_k_supernuclear_cgs + elsewhere + local_eos_gamma = eos_gamma + local_eos_k_cgs = eos_k_cgs + end where + +! First calculate the pressure + 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 + + pressure = p_poly + p_th + + p_th = max(zero, p_th) + +! This may look incorrect; It's in Harry's code; in his thesis and I worked it out as well. +! Should be okay... + cs2 = (local_eos_gamma * p_poly + eos_gamma_th * p_th) / & + rho / (1.0d0 + eps + pressure/rho) + + deallocate(pressure) + deallocate(p_poly) + deallocate(p_th) + deallocate(local_eos_gamma) + deallocate(local_eos_k_cgs) + +end subroutine EOS_GeneralHybrid_cs2 + + + +subroutine EOS_GeneralHybrid_DPDIE(nelems, rho, eps, dpdie) + + + USE EOS_GP_Scalars + USE EOS_GeneralHybrid_Scalars + + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: nelems + CCTK_REAL, dimension(nelems), intent(in) :: rho + CCTK_REAL, dimension(nelems), intent(in) :: eps + CCTK_REAL, dimension(nelems), intent(out) :: dpdie + + dpdie = (eos_gamma_th - 1.d0) * rho + +end subroutine EOS_GeneralHybrid_DPDIE + |