Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2017-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module eospc_eval
21 : use eos_def
22 : use auto_diff
23 : use const_def, only: dp, qe, qp, avo, crad, ln10, arg_not_provided, amu, kerg, one_third, four_thirds_pi
24 : use utils_lib, only: is_bad, mesa_error
25 : use math_lib
26 :
27 : implicit none
28 :
29 : contains
30 :
31 0 : subroutine Get_PC_alfa( &
32 : rq, logRho, logT, Z, abar, zbar, &
33 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
34 : ierr)
35 : use const_def, only: dp
36 : type (EoS_General_Info), pointer :: rq
37 : real(dp), intent(in) :: logRho, logT, Z, abar, zbar
38 : real(dp), intent(out) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
39 : integer, intent(out) :: ierr
40 0 : real(dp) :: logGe0, logGe, logGe_lo, logGe_hi, &
41 0 : A, B, dA_dlnT, dA_dlnRho, dB_dlnT, dB_dlnRho, dlogGe_dlogT, dlogGe_dlogRho, &
42 0 : logT_lo, logT_hi, logRho_lo, logRho_hi
43 :
44 : include 'formats'
45 :
46 0 : ierr = 0
47 :
48 0 : d_alfa_dlogT = 0d0
49 0 : d_alfa_dlogRho = 0d0
50 :
51 0 : logRho_lo = rq% logRho2_PC_limit ! don't use PC for logRho < this
52 0 : logRho_hi = rq% logRho1_PC_limit ! okay for pure PC for logRho > this
53 :
54 0 : if (rq% PC_use_Gamma_limit_instead_of_T) then
55 : !gamma_e = (qe**2)*(four_thirds_pi*avo*Rho*zbar/abar)**one_third/(kerg*T)
56 : !logGe = logGe0 + logRho/3 - logT
57 : ! where Ge0 = (qe**2)*(four_thirds_pi*avo*zbar/abar)**one_third/kerg
58 : logGe0 = log10( &
59 0 : qe*qe*pow(four_thirds_pi*avo*zbar/abar, one_third)/kerg)
60 0 : logGe = logGe0 + logRho/3 - logT
61 0 : logGe_lo = rq% log_Gamma_e_all_HELM ! HELM for logGe <= this
62 0 : logGe_hi = rq% log_Gamma_e_all_PC ! PC for logGe >= this
63 0 : logT_lo = logGe0 + logRho_lo/3 - logGe_lo
64 0 : logT_hi = logGe0 + logRho_hi/3 - logGe_hi
65 0 : if (logRho <= logRho_lo .or. logGe <= logGe_lo) then
66 0 : alfa = 1d0 ! no PC
67 0 : else if (logRho >= logRho_hi .and. logGe >= logGe_hi) then
68 0 : alfa = 0d0 ! pure PC
69 0 : else if (logT >= logT_hi) then ! blend in logGe
70 0 : alfa = (logGe_hi - logGe)/(logGe_hi - logGe_lo)
71 0 : dlogGe_dlogT = -1d0
72 0 : dlogGe_dlogRho = 1d0/3d0
73 0 : d_alfa_dlogT = -dlogGe_dlogT/(logGe_hi - logGe_lo)
74 0 : d_alfa_dlogRho = -dlogGe_dlogRho/(logGe_hi - logGe_lo)
75 : else ! blend in logRho
76 0 : if (logT >= logT_lo) logRho_lo = (logT_lo + logGe_lo - logGe0)*3
77 0 : alfa = (logRho_hi - logRho)/(logRho_hi - logRho_lo)
78 0 : d_alfa_dlogRho = -1d0/(logRho_hi - logRho_lo)
79 : end if
80 : else
81 0 : logT_lo = rq% logT1_PC_limit ! okay for pure PC for logT < this (like logT_all_OPAL)
82 0 : logT_hi = rq% logT2_PC_limit ! don't use PC for logT > this (like logT_all_HELM)
83 0 : if (logRho <= logRho_lo .or. logT >= logT_hi) then
84 0 : alfa = 1d0 ! no PC
85 0 : else if (logRho >= logRho_hi .and. logT <= logT_lo) then
86 0 : alfa = 0d0 ! pure PC
87 0 : else if (logT >= logT_lo) then ! blend in logT
88 0 : alfa = (logT - logT_lo)/(logT_hi - logT_lo)
89 0 : d_alfa_dlogT = 1/(logT_hi - logT_lo)
90 : else ! blend in logRho
91 0 : alfa = (logRho_hi - logRho)/(logRho_hi - logRho_lo)
92 0 : d_alfa_dlogRho = -1d0/(logRho_hi - logRho_lo)
93 : end if
94 : end if
95 :
96 0 : end subroutine Get_PC_alfa
97 :
98 :
99 0 : subroutine get_PC_for_eosdt( &
100 : handle, dbg, Z, X, abar, zbar, &
101 0 : species, chem_id, net_iso, xa, &
102 : rho, logRho, T, logT, remaining_fraction, &
103 0 : res, d_dlnd, d_dlnT, d_dxa, &
104 : skip, ierr)
105 : integer, intent(in) :: handle
106 : logical, intent(in) :: dbg
107 : real(dp), intent(in) :: &
108 : Z, X, abar, zbar, remaining_fraction
109 : integer, intent(in) :: species
110 : integer, pointer :: chem_id(:), net_iso(:)
111 : real(dp), intent(in) :: xa(:)
112 : real(dp), intent(in) :: rho, logRho, T, logT
113 : real(dp), intent(inout), dimension(nv) :: &
114 : res, d_dlnd, d_dlnT
115 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
116 : logical, intent(out) :: skip
117 : integer, intent(out) :: ierr
118 : type (EoS_General_Info), pointer :: rq
119 0 : rq => eos_handles(handle)
120 : call Get_PC_Results(rq, &
121 : Z, X, abar, zbar, &
122 : species, chem_id, net_iso, xa, &
123 : rho, logRho, T, logT, &
124 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
125 0 : skip = .false.
126 :
127 : ! zero latent heat information
128 0 : res(i_latent_ddlnT:i_latent_ddlnRho) = 0d0
129 0 : d_dlnT(i_latent_ddlnT:i_latent_ddlnRho) = 0d0
130 0 : d_dlnd(i_latent_ddlnT:i_latent_ddlnRho) = 0d0
131 :
132 : ! zero all components
133 0 : res(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
134 0 : d_dlnd(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
135 0 : d_dlnT(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
136 :
137 : ! mark this one
138 0 : res(i_frac_PC) = 1.0d0
139 :
140 0 : end subroutine get_PC_for_eosdt
141 :
142 :
143 0 : subroutine Get_PC_Results( &
144 : rq, Z, X, abar, zbar, &
145 0 : species, chem_id, net_iso, xa, &
146 : Rho, logRho, T, logT, &
147 0 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa, ierr)
148 : use pc_eos
149 : use chem_def, only: chem_isos
150 : type (EoS_General_Info), pointer :: rq
151 : real(dp), intent(in) :: Z, X, abar, zbar
152 : integer, intent(in) :: species
153 : integer, pointer :: chem_id(:), net_iso(:)
154 : real(dp), intent(in) :: xa(:)
155 : real(dp), intent(in) :: Rho, logRho, T, logT
156 : real(dp), intent(inout) :: res(:) ! (nv)
157 : real(dp), intent(inout) :: d_dlnRho_c_T(:) ! (nv)
158 : real(dp), intent(inout) :: d_dlnT_c_Rho(:) ! (nv)
159 : real(dp), intent(inout) :: d_dxa(:,:) ! (nv, species)
160 : integer, intent(out) :: ierr
161 :
162 0 : real(dp) :: start_crystal, full_crystal
163 0 : real(dp), dimension(species) :: AY, AZion, ACMI
164 : integer :: i, j
165 :
166 : include 'formats'
167 :
168 0 : ierr = 0
169 0 : AZion(1:species) = chem_isos% Z(chem_id(1:species))
170 0 : ACMI(1:species) = chem_isos% W(chem_id(1:species)) ! this really is atomic weight.
171 0 : do j=1,species
172 0 : if (xa(j) < rq% mass_fraction_limit_for_PC) then
173 0 : AY(j) = 0
174 : else
175 0 : AY(j) = xa(j)/ACMI(j)
176 : end if
177 : end do
178 :
179 0 : start_crystal = rq% PC_Gamma_start_crystal
180 0 : full_crystal = rq% PC_Gamma_full_crystal
181 :
182 :
183 0 : call do1(.false.,Rho,T,start_crystal,full_crystal,res,d_dlnT_c_Rho,d_dlnRho_c_T,ierr)
184 0 : if (ierr /= 0) return
185 :
186 : ! composition derivatives not provided
187 0 : d_dxa = 0
188 :
189 0 : if (is_bad(res(i_lnS))) then
190 0 : ierr = -1
191 0 : write(*,1) 'res(i_lnS), logRho, logT', res(i_lnS), logRho, logT
192 0 : call mesa_error(__FILE__,__LINE__,'Get_PC_Results')
193 : end if
194 :
195 : contains
196 :
197 0 : subroutine do1(show,RHO_real,T_real,start_crystal,full_crystal,res,d_dlnT_c_Rho,d_dlnRho_c_T,ierr)
198 : logical, intent(in) :: show
199 : real(dp), intent(in) :: start_crystal, full_crystal
200 : real(dp), intent(in) :: RHO_real, T_real
201 : real(dp), intent(out) :: res(:), d_dlnT_c_Rho(:), d_dlnRho_c_T(:) ! (nv)
202 : integer, intent(out) :: ierr
203 :
204 : integer :: j, LIQSOL
205 : type(auto_diff_real_2var_order1) :: dsp, dse, dpe, &
206 : DENS,GAMI,CHI,TPT,GAMImean,TEMP, &
207 : PnkT,UNkT,SNk,CVNkt,CV,CHIR,CHIT,Tnk,P,Cp,gamma1,gamma3,grad_ad, N, &
208 : Prad,Pgas,PRADnkT,dE_dRho,dS_dT,dS_dRho,mu,lnfree_e, lnPgas, lnE, lnS, RHO, T
209 : type(auto_diff_real_2var_order1) :: phase
210 0 : real(dp) :: Zmean,CMImean,Z2mean
211 : real(dp), parameter :: UN_T6=0.3157746d0
212 :
213 : include 'formats'
214 :
215 0 : ierr = 0
216 :
217 0 : T = T_real
218 0 : T%d1val1 = 1d0
219 :
220 0 : RHO = RHO_real
221 0 : RHO%d1val2 = 1d0
222 :
223 0 : TEMP=T*1d-6/UN_T6 ! T [au]
224 :
225 0 : if (show) then
226 0 : write(*,1) 'RHO', RHO
227 0 : write(*,1) 'TEMP', TEMP
228 0 : write(*,1) 'start_crystal', start_crystal
229 0 : write(*,1) 'full_crystal', full_crystal
230 0 : write(*,1) 'AZion(1:species)', AZion(1:species)
231 0 : write(*,1) 'ACMI(1:species)', ACMI(1:species)
232 0 : write(*,1) 'AY(1:species)', AY(1:species)
233 0 : write(*,1) 'xa(1:species)', xa(1:species)
234 : end if
235 :
236 : call MELANGE9( &
237 : species,AY,AZion,ACMI,RHO,TEMP, &
238 : start_crystal,full_crystal,PRADnkT, &
239 : DENS,Zmean,CMImean,Z2mean,GAMImean,CHI,TPT,LIQSOL, &
240 0 : PnkT,UNkT,SNk,CVNkt,CHIR,CHIT,ierr)
241 :
242 : ! calculate phase information
243 0 : if (GAMImean<start_crystal) then ! liquid
244 0 : phase = 0d0
245 0 : else if (GAMImean>full_crystal) then ! solid
246 0 : phase = 1d0
247 : else ! blend of liquid and solid
248 0 : phase = (GAMImean - start_crystal)/(full_crystal - start_crystal) ! 1 for solid, 0 for liquid
249 : end if
250 :
251 0 : if (ierr /= 0) then
252 0 : return
253 : write(*,1) 'RHO', RHO
254 : write(*,1) 'T', T
255 : write(*,1) 'logRho', log10(RHO)
256 : write(*,1) 'logT', log10(T)
257 : write(*,*) 'ierr from MELANGE9'
258 : call mesa_error(__FILE__,__LINE__,'debug eos')
259 : end if
260 :
261 0 : if (show) then
262 0 : write(*,1) 'PRADnkT', PRADnkT
263 0 : write(*,1) 'DENS', DENS
264 0 : write(*,1) 'Zmean', Zmean
265 0 : write(*,1) 'CMImean', CMImean
266 0 : write(*,1) 'Z2mean', Z2mean
267 0 : write(*,1) 'GAMI', GAMImean
268 0 : write(*,1) 'CHI', CHI
269 0 : write(*,1) 'TPT', TPT
270 0 : write(*,1) 'PnkT', PnkT
271 0 : write(*,1) 'UNkT', UNkT
272 0 : write(*,1) 'SNk', SNk
273 0 : write(*,1) 'CV', CVNkt
274 0 : write(*,1) 'CHIR', CHIR
275 0 : write(*,1) 'CHIT', CHIT
276 0 : write(*,'(A)')
277 : end if
278 :
279 0 : Tnk=8.31447d7/CMImean*RHO*T ! n_i kT [erg/cc]
280 0 : Pgas = PnkT*Tnk
281 0 : if (rq% include_radiation) then
282 : ! results from MELANGE9 do not include radiation. add it now.
283 0 : CHIT=CHIT*(PnkT/(PnkT+PRADnkT)) + 4.d0*PRADnkT/(PnkT+PRADnkT)
284 0 : CHIR=CHIR*(PnkT/(PnkT+PRADnkT))
285 0 : PnkT=PnkT+PRADnkT
286 0 : UNkT=UNkT+3.d0*PRADnkT
287 0 : SNk=SNk+4.d0*PRADnkT
288 0 : CVNkt=CVNkt+12.d0*PRADnkT
289 : end if
290 0 : P = PnkT*Tnk
291 0 : N = 1/(abar*amu)
292 0 : CV = CVNkt*N*kerg
293 0 : gamma3 = 1d0 + P/rho * chit/(T*CV)
294 0 : gamma1 = chit*(gamma3-1d0) + chir
295 0 : grad_ad = (gamma3-1d0)/gamma1
296 0 : Cp = CV * gamma1/chir
297 0 : dE_dRho = (1d0-chiT)*P/(rho*rho)
298 0 : dS_dT = CV/T
299 0 : dS_dRho = -P*chiT/(rho*rho * T)
300 0 : mu = abar / (1d0 + zbar)
301 0 : lnfree_e = log(zbar/abar) ! for complete ionization
302 0 : lnPgas = log(Pgas)
303 0 : lnE = safe_log(UNkt*N*kerg*T)
304 0 : lnS = safe_log(SNk*N*kerg)
305 :
306 :
307 0 : res(i_lnPgas) = lnPgas % val
308 0 : res(i_lnE) = lnE % val
309 0 : res(i_lnS) = lnS % val
310 0 : res(i_grad_ad) = grad_ad % val
311 0 : res(i_chiRho) = CHIR % val
312 0 : res(i_chiT) = CHIT % val
313 0 : res(i_Cp) = Cp % val
314 0 : res(i_Cv) = CV % val
315 0 : res(i_dE_dRho) = dE_dRho % val
316 0 : res(i_dS_dT) = dS_dT % val
317 0 : res(i_dS_dRho) = dS_dRho % val
318 0 : res(i_mu) = mu % val
319 0 : res(i_lnfree_e) = lnfree_e % val
320 0 : res(i_gamma1) = gamma1 % val
321 0 : res(i_gamma3) = gamma3 % val
322 0 : res(i_eta) = CHI % val
323 0 : res(i_phase) = phase % val
324 :
325 0 : d_dlnT_c_Rho(i_lnPgas) = lnPgas % d1val1 * T % val
326 0 : d_dlnT_c_Rho(i_lnE) = lnE % d1val1 * T % val
327 0 : d_dlnT_c_Rho(i_lnS) = lnS % d1val1 * T % val
328 0 : d_dlnT_c_Rho(i_grad_ad) = grad_ad % d1val1 * T % val
329 0 : d_dlnT_c_Rho(i_chiRho) = CHIR % d1val1 * T % val
330 0 : d_dlnT_c_Rho(i_chiT) = CHIT % d1val1 * T % val
331 0 : d_dlnT_c_Rho(i_Cp) = Cp % d1val1 * T % val
332 0 : d_dlnT_c_Rho(i_Cv) = CV % d1val1 * T % val
333 0 : d_dlnT_c_Rho(i_dE_dRho) = dE_dRho % d1val1 * T % val
334 0 : d_dlnT_c_Rho(i_dS_dT) = dS_dT % d1val1 * T % val
335 0 : d_dlnT_c_Rho(i_dS_dRho) = dS_dRho % d1val1 * T % val
336 0 : d_dlnT_c_Rho(i_mu) = mu % d1val1 * T % val
337 0 : d_dlnT_c_Rho(i_lnfree_e) = lnfree_e % d1val1 * T % val
338 0 : d_dlnT_c_Rho(i_gamma1) = gamma1 % d1val1 * T % val
339 0 : d_dlnT_c_Rho(i_gamma3) = gamma3 % d1val1 * T % val
340 0 : d_dlnT_c_Rho(i_eta) = CHI % d1val1 * T % val
341 0 : d_dlnT_c_Rho(i_phase) = phase % d1val1 * T % val
342 :
343 0 : d_dlnRho_c_T(i_lnPgas) = lnPgas % d1val2 * RHO % val
344 0 : d_dlnRho_c_T(i_lnE) = lnE % d1val2 * RHO % val
345 0 : d_dlnRho_c_T(i_lnS) = lnS % d1val2 * RHO % val
346 0 : d_dlnRho_c_T(i_grad_ad) = grad_ad % d1val2 * RHO % val
347 0 : d_dlnRho_c_T(i_chiRho) = CHIR % d1val2 * RHO % val
348 0 : d_dlnRho_c_T(i_chiT) = CHIT % d1val2 * RHO % val
349 0 : d_dlnRho_c_T(i_Cp) = Cp % d1val2 * RHO % val
350 0 : d_dlnRho_c_T(i_Cv) = CV % d1val2 * RHO % val
351 0 : d_dlnRho_c_T(i_dE_dRho) = dE_dRho % d1val2 * RHO % val
352 0 : d_dlnRho_c_T(i_dS_dT) = dS_dT % d1val2 * RHO % val
353 0 : d_dlnRho_c_T(i_dS_dRho) = dS_dRho % d1val2 * RHO % val
354 0 : d_dlnRho_c_T(i_mu) = mu % d1val2 * RHO % val
355 0 : d_dlnRho_c_T(i_lnfree_e) = lnfree_e % d1val2 * RHO % val
356 0 : d_dlnRho_c_T(i_gamma1) = gamma1 % d1val2 * RHO % val
357 0 : d_dlnRho_c_T(i_gamma3) = gamma3 % d1val2 * RHO % val
358 0 : d_dlnRho_c_T(i_eta) = CHI % d1val2 * RHO % val
359 0 : d_dlnRho_c_T(i_phase) = phase % d1val2 * RHO % val
360 :
361 0 : end subroutine do1
362 :
363 : end subroutine Get_PC_Results
364 :
365 : end module eospc_eval
|
