1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
/*@@
@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_DECLARE(CCTK_REAL, rho, )
CCTK_DECLARE(CCTK_REAL, 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_DECLARE(CCTK_REAL, eps, )
CCTK_DECLARE(CCTK_REAL, 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
CCTK_DECLARE(CCTK_REAL, eps, )
EOS_Hybrid_DPressByDEps = (eos_gamma_th - 1.d0) * rho
end function EOS_Hybrid_DPressByDEps
|
