Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2022 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 skye
21 : use const_def, only: dp, crad, kerg, mp
22 : use math_lib
23 : use auto_diff
24 : use eos_def
25 :
26 : implicit none
27 :
28 : logical, parameter :: dbg = .false.
29 :
30 :
31 : private
32 : public :: Get_Skye_EOS_Results, Get_Skye_alfa, Get_Skye_alfa_simple, get_Skye_for_eosdt
33 :
34 : contains
35 :
36 286772 : subroutine Get_Skye_alfa( &
37 : rq, logRho, logT, Z, abar, zbar, &
38 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
39 : ierr)
40 : use const_def, only: dp
41 : use eos_blend
42 : type (EoS_General_Info), pointer :: rq
43 : real(dp), intent(in) :: logRho, logT, Z, abar, zbar
44 : real(dp), intent(out) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
45 : integer, intent(out) :: ierr
46 :
47 : logical :: contained
48 : type(auto_diff_real_2var_order1) :: p(2), blend, dist
49 :
50 : ! Blend parameters
51 241252 : real(dp) :: skye_blend_width
52 : integer, parameter :: num_points = 8
53 4583788 : real(dp) :: bounds(8,2)
54 : type (Helm_Table), pointer :: ht
55 :
56 241252 : ierr = 0
57 241252 : ht => eos_ht
58 241252 : skye_blend_width = 0.1d0
59 :
60 : ! Avoid catastrophic loss of precision in HELM tables
61 241252 : bounds(1,1) = ht% logdlo
62 241252 : bounds(1,2) = 8.3d0
63 :
64 : ! Rough ionization temperature from Jermyn+2021 Equation 52 (treating denominator as ~1).
65 : ! We put a lower bound of logT=7.3 to ensure that solar models never use Skye.
66 : ! This is because the blend even in regions that are 99+% ionized produces noticeable
67 : ! kinks in the sound speed profile on a scale testable by the observations.
68 241252 : bounds(2,1) = ht% logdlo
69 241252 : bounds(2,2) = max(7.3d0,log10(1d5 * pow2(zbar))) + skye_blend_width
70 :
71 : ! Rough ionization density from Jermyn+2021 Equation 53, dividing by 3 so we get closer to Dragons.
72 : ! Don't let the density get below rho = 200 g/cc so that solar model stays away from blend into Skye.
73 241252 : bounds(3,1) = max(2.3d0,log10(abar * pow3(zbar))) + skye_blend_width
74 241252 : bounds(3,2) = max(7.3d0,log10(1d5 * pow2(zbar))) + skye_blend_width
75 :
76 : ! HELM low-T bound
77 241252 : bounds(4,1) = max(2.3d0,log10(abar * pow3(zbar))) + skye_blend_width
78 241252 : bounds(4,2) = ht% logtlo
79 :
80 : ! Lower-right of (rho,T) plane
81 241252 : bounds(5,1) = ht% logdhi
82 241252 : bounds(5,2) = ht% logtlo
83 :
84 : ! Upper-right of (rho,T) plane
85 241252 : bounds(6,1) = ht% logdhi
86 241252 : bounds(6,2) = ht% logthi
87 :
88 : ! Avoid catastrophic loss of precision in HELM tables
89 241252 : bounds(7,1) = 3d0 * ht% logthi + log10(abar * mp * crad / (3d0 * kerg * (zbar + 1d0))) - 6d0
90 241252 : bounds(7,2) = ht% logthi
91 :
92 : ! Avoid catastrophic loss of precision in HELM tables
93 241252 : bounds(8,1) = 3d0 * 8.3d0 + log10(abar * mp * crad / (3d0 * kerg * (zbar + 1d0))) - 6d0
94 241252 : bounds(8,2) = 8.3d0
95 :
96 : ! Set up auto_diff point
97 241252 : p(1) = logRho
98 241252 : p(1)%d1val1 = 1d0
99 241252 : p(2) = logT
100 241252 : p(2)%d1val2 = 1d0
101 :
102 241252 : contained = is_contained(num_points, bounds, p)
103 241252 : dist = min_distance_to_polygon(num_points, bounds, p)
104 :
105 241252 : if (contained) then ! Make distance negative for points inside the polygon
106 30586 : dist = -dist
107 : end if
108 :
109 241252 : dist = dist / skye_blend_width
110 241252 : blend = max(dist, 0d0)
111 241252 : blend = min(blend, 1d0)
112 :
113 241252 : alfa = blend%val
114 241252 : d_alfa_dlogRho = blend%d1val1
115 241252 : d_alfa_dlogT = blend%d1val2
116 :
117 241252 : end subroutine Get_Skye_alfa
118 :
119 :
120 0 : subroutine Get_Skye_alfa_simple( &
121 : rq, logRho, logT, Z, abar, zbar, &
122 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
123 : ierr)
124 241252 : use eos_blend
125 : type (EoS_General_Info), pointer :: rq
126 : real(dp), intent(in) :: logRho, logT, Z, abar, zbar
127 : real(dp), intent(out) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
128 : integer, intent(out) :: ierr
129 :
130 : type(auto_diff_real_2var_order1) :: logT_auto, logRho_auto
131 : type(auto_diff_real_2var_order1) :: blend, blend_logT, blend_logRho
132 :
133 : include 'formats'
134 :
135 0 : ierr = 0
136 :
137 : ! logRho is val1
138 0 : logRho_auto% val = logRho
139 0 : logRho_auto% d1val1 = 1d0
140 0 : logRho_auto% d1val2 = 0d0
141 :
142 : ! logT is val2
143 0 : logT_auto% val = logT
144 0 : logT_auto% d1val1 = 0d0
145 0 : logT_auto% d1val2 = 1d0
146 :
147 : ! logT blend
148 0 : if (logT_auto < rq% logT_min_for_any_Skye) then
149 0 : blend_logT = 0d0
150 0 : else if (logT_auto <= rq% logT_min_for_all_Skye) then
151 0 : blend_logT = (logT_auto - rQ% logT_min_for_any_Skye) / (rq% logT_min_for_all_Skye - rq% logT_min_for_any_Skye)
152 0 : else if (logT_auto > rq% logT_min_for_all_Skye) then
153 0 : blend_logT = 1d0
154 : end if
155 :
156 :
157 : ! logRho blend
158 0 : if (logRho_auto < rq% logRho_min_for_any_Skye) then
159 0 : blend_logRho = 0d0
160 0 : else if (logRho_auto <= rq% logRho_min_for_all_Skye) then
161 0 : blend_logRho = (logRho_auto - rQ% logRho_min_for_any_Skye) / (rq% logRho_min_for_all_Skye - rq% logRho_min_for_any_Skye)
162 0 : else if (logRho_auto > rq% logRho_min_for_all_Skye) then
163 0 : blend_logRho = 1d0
164 : end if
165 :
166 : ! combine blends
167 0 : blend = (1d0 - blend_logRho) * (1d0 - blend_logT)
168 :
169 0 : alfa = blend% val
170 0 : d_alfa_dlogRho = blend% d1val1
171 0 : d_alfa_dlogT = blend% d1val2
172 :
173 0 : end subroutine get_Skye_alfa_simple
174 :
175 :
176 45520 : subroutine get_Skye_for_eosdt(handle, dbg, Z, X, abar, zbar, species, chem_id, net_iso, xa, &
177 45520 : rho, logRho, T, logT, remaining_fraction, res, d_dlnd, d_dlnT, d_dxa, skip, ierr)
178 : integer, intent(in) :: handle
179 : logical, intent(in) :: dbg
180 : real(dp), intent(in) :: &
181 : Z, X, abar, zbar, remaining_fraction
182 : integer, intent(in) :: species
183 : integer, pointer :: chem_id(:), net_iso(:)
184 : real(dp), intent(in) :: xa(:)
185 : real(dp), intent(in) :: rho, logRho, T, logT
186 : real(dp), intent(inout), dimension(nv) :: res, d_dlnd, d_dlnT
187 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
188 : logical, intent(out) :: skip
189 : integer, intent(out) :: ierr
190 : type (EoS_General_Info), pointer :: rq
191 :
192 45520 : rq => eos_handles(handle)
193 :
194 : call Get_Skye_EOS_Results(rq, Z, X, abar, zbar, rho, logRho, T, logT, species, chem_id, xa, &
195 45520 : res, d_dlnd, d_dlnT, d_dxa, ierr)
196 45520 : skip = .false.
197 :
198 : ! zero all components
199 364160 : res(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
200 364160 : d_dlnd(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
201 364160 : d_dlnT(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
202 :
203 : ! mark this one
204 45520 : res(i_frac_Skye) = 1.0d0
205 :
206 0 : end subroutine get_Skye_for_eosdt
207 :
208 45520 : subroutine Get_Skye_EOS_Results( &
209 : rq, Z, X, abar, zbar, Rho, logRho, T, logT, &
210 45520 : species, chem_id, xa, res, d_dlnd, d_dlnT, d_dxa, ierr)
211 : type (EoS_General_Info), pointer :: rq
212 : real(dp), intent(in) :: Z, X, abar, zbar
213 : real(dp), intent(in) :: Rho, logRho, T, logT
214 : integer, intent(in) :: species
215 : integer, pointer :: chem_id(:)
216 : real(dp), intent(in) :: xa(:)
217 : integer, intent(out) :: ierr
218 : real(dp), intent(out), dimension(nv) :: res, d_dlnd, d_dlnT
219 : real(dp), intent(out), dimension(nv, species) :: d_dxa
220 :
221 : real(dp) :: logT_ion, logT_neutral
222 :
223 : include 'formats'
224 :
225 : ierr = 0
226 :
227 : call skye_eos( &
228 : T, Rho, X, abar, zbar, &
229 : rq%Skye_min_gamma_for_solid, rq%Skye_max_gamma_for_liquid, &
230 : rq%Skye_solid_mixing_rule, rq%mass_fraction_limit_for_Skye, &
231 : rq%Skye_use_ion_offsets, &
232 : species, chem_id, xa, &
233 45520 : res, d_dlnd, d_dlnT, d_dxa, ierr)
234 :
235 : ! composition derivatives not provided
236 19315690 : d_dxa = 0
237 :
238 : if (ierr /= 0) then
239 : if (dbg) then
240 : write(*,*) 'failed in Get_Skye_EOS_Results'
241 : write(*,1) 'T', T
242 : write(*,1) 'logT', logT
243 : write(*,1) 'Rho', Rho
244 : write(*,1) 'logRho', logRho
245 : write(*,1) 'abar', abar
246 : write(*,1) 'zbar', zbar
247 : write(*,1) 'X', X
248 : call mesa_error(__FILE__,__LINE__,'Get_Skye_EOS_Results')
249 : end if
250 : return
251 : end if
252 :
253 : end subroutine Get_Skye_EOS_Results
254 :
255 :
256 : !>..given a temperature temp [K], density den [g/cm**3], and a composition
257 : !!..this routine returns most of the other
258 : !!..thermodynamic quantities. of prime interest is the pressure [erg/cm**3],
259 : !!..specific thermal energy [erg/gr], the entropy [erg/g/K], along with
260 : !!..their derivatives with respect to temperature, density, abar, and zbar.
261 : !!..other quantities such the normalized chemical potential eta (plus its
262 : !!..derivatives), number density of electrons and positron pair (along
263 : !!..with their derivatives), adiabatic indices, specific heats, and
264 : !!..relativistically correct sound speed are also returned.
265 : !!..
266 : !!..this routine assumes planckian photons, an ideal gas of ions,
267 : !!..and an electron-positron gas with an arbitrary degree of relativity
268 : !!..and degeneracy. interpolation in a table of the helmholtz free energy
269 : !!..is used to return the electron-positron thermodynamic quantities.
270 : !!..all other derivatives are analytic.
271 : !!..
272 : !!..references: cox & giuli chapter 24 ; timmes & swesty apj 1999
273 :
274 : !!..this routine assumes a call to subroutine read_helm_table has
275 : !!..been performed prior to calling this routine.
276 45520 : subroutine skye_eos( &
277 : temp_in, den_in, Xfrac, abar, zbar, &
278 : Skye_min_gamma_for_solid, Skye_max_gamma_for_liquid, &
279 : Skye_solid_mixing_rule, &
280 : mass_fraction_limit, use_ion_offsets, &
281 45520 : species, chem_id, xa, &
282 : res, d_dlnd, d_dlnT, d_dxa, ierr)
283 :
284 : use eos_def
285 : use utils_lib, only: is_bad
286 : use chem_def, only: chem_isos
287 : use ion_offset, only: compute_ion_offset
288 : use skye_ideal
289 : use skye_coulomb
290 : use skye_thermodynamics
291 : use auto_diff
292 :
293 : integer :: j
294 : integer, intent(in) :: species
295 : integer, pointer :: chem_id(:)
296 : real(dp), intent(in) :: xa(:)
297 : real(dp), intent(in) :: temp_in, den_in, mass_fraction_limit, Skye_min_gamma_for_solid, Skye_max_gamma_for_liquid
298 : real(dp), intent(in) :: Xfrac, abar, zbar
299 : logical, intent(in) :: use_ion_offsets
300 : character(len=128), intent(in) :: Skye_solid_mixing_rule
301 : integer, intent(out) :: ierr
302 : real(dp), intent(out), dimension(nv) :: res, d_dlnd, d_dlnT
303 : real(dp), intent(out), dimension(nv, species) :: d_dxa
304 :
305 : integer :: relevant_species, lookup(species)
306 : type(auto_diff_real_2var_order3) :: temp, logtemp, den, logden, din
307 3659590 : real(dp) :: AZION(species), ACMI(species), A(species), select_xa(species), ya(species)
308 : type (Helm_Table), pointer :: ht
309 759230 : real(dp) :: ytot1, ye, norm
310 : type(auto_diff_real_2var_order3) :: etaele, xnefer, phase, latent_ddlnT, latent_ddlnRho
311 : type(auto_diff_real_2var_order3) :: F_ion_gas, F_rad, F_ideal_ion, F_coul
312 : type(auto_diff_real_2var_order3) :: F_ele
313 :
314 45520 : ht => eos_ht
315 :
316 45520 : ierr = 0
317 :
318 45520 : temp = temp_in
319 45520 : temp%d1val1 = 1d0
320 45520 : logtemp = log10(temp)
321 :
322 45520 : den = den_in
323 45520 : den%d1val2 = 1d0
324 45520 : logden = log10(den)
325 :
326 : ! HELM table lookup uses din rather than den
327 45520 : ytot1 = 1.0d0 / abar
328 45520 : ye = ytot1 * zbar
329 45520 : din = ye*den
330 :
331 45520 : F_rad = 0d0
332 45520 : F_ion_gas = 0d0
333 45520 : F_ideal_ion = 0d0
334 45520 : F_coul = 0d0
335 45520 : F_ele = 0d0
336 :
337 : ! Radiation free energy, independent of composition
338 45520 : F_rad = compute_F_rad(temp, den)
339 :
340 : ! Count and pack relevant species for Coulomb corrections. Relevant means mass fraction above limit.
341 45520 : relevant_species = 0
342 45520 : norm = 0d0
343 759230 : do j=1,species
344 759230 : if (xa(j) > mass_fraction_limit) then
345 310914 : relevant_species = relevant_species + 1
346 310914 : AZION(relevant_species) = chem_isos% Z(chem_id(j))
347 310914 : ACMI(relevant_species) = chem_isos% W(chem_id(j))
348 310914 : A(relevant_species) = chem_isos% Z_plus_N(chem_id(j))
349 310914 : select_xa(relevant_species) = xa(j)
350 310914 : norm = norm + xa(j)
351 : end if
352 : end do
353 :
354 : ! Normalize
355 356434 : do j=1,relevant_species
356 356434 : select_xa(j) = select_xa(j) / norm
357 : end do
358 :
359 : ! Compute number fractions
360 : norm = 0d0
361 356434 : do j=1,relevant_species
362 310914 : ya(j) = select_xa(j) / A(j)
363 356434 : norm = norm + ya(j)
364 : end do
365 356434 : do j=1,relevant_species
366 356434 : ya(j) = ya(j) / norm
367 : end do
368 :
369 : ! Ideal ion free energy, only depends on abar
370 45520 : F_ideal_ion = compute_F_ideal_ion(temp, den, abar, relevant_species, ACMI, ya)
371 :
372 45520 : if (use_ion_offsets) then
373 45520 : F_ideal_ion = F_ideal_ion + compute_ion_offset(species, xa, chem_id) ! Offset so ion ground state energy is zero.
374 : end if
375 :
376 : ! Ideal electron-positron thermodynamics (s, e, p)
377 : ! Derivatives are handled by HELM code, so we don't pass *in* any auto_diff types (just get them as return values).
378 : call compute_ideal_ele(temp%val, den%val, din%val, logtemp%val, logden%val, zbar, ytot1, ye, ht, &
379 45520 : F_ele, etaele, xnefer, ierr)
380 :
381 45520 : xnefer = compute_xne(den, ytot1, zbar)
382 :
383 : ! Normalize mass fractions
384 356434 : do j=1,relevant_species
385 356434 : select_xa(j) = select_xa(j) / norm
386 : end do
387 :
388 : ! Compute non-ideal corrections
389 : call nonideal_corrections(relevant_species, ya(1:relevant_species), &
390 : AZION(1:relevant_species), ACMI(1:relevant_species), &
391 : Skye_min_gamma_for_solid, Skye_max_gamma_for_liquid, &
392 : Skye_solid_mixing_rule, den, temp, xnefer, abar, &
393 45520 : F_coul, latent_ddlnT, latent_ddlnRho, phase)
394 :
395 : call pack_for_export(F_ideal_ion, F_coul, F_rad, F_ele, temp, den, xnefer, etaele, abar, zbar, &
396 45520 : phase, latent_ddlnT, latent_ddlnRho, res, d_dlnd, d_dlnT, ierr)
397 45520 : if(ierr/=0) return
398 :
399 45520 : end subroutine skye_eos
400 :
401 :
402 : end module skye
|
