Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2021 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 eos_lib
21 :
22 : use const_def, only: dp, crad
23 : use math_lib
24 :
25 : implicit none
26 :
27 : contains
28 :
29 15 : subroutine eos_init( &
30 : eosDT_cache_dir, use_cache, info)
31 : use eos_initialize, only : Init_eos
32 : character(*), intent(in) :: eosDT_cache_dir ! blank string means use default
33 : logical, intent(in) :: use_cache
34 : integer, intent(out) :: info ! 0 means AOK.
35 : info = 0
36 : call Init_eos( &
37 15 : eosDT_cache_dir, use_cache, info)
38 : if (info /= 0) return
39 : end subroutine eos_init
40 :
41 5 : subroutine eos_shutdown
42 : use eos_def
43 : use helm_alloc,only:free_helm_table
44 5 : if (associated(eos_ht)) call free_helm_table(eos_ht)
45 5 : call eos_def_shutdown
46 5 : end subroutine eos_shutdown
47 :
48 :
49 : ! after eos_init has finished, you can allocate a "handle"
50 : ! and set control parameter values using an inlist
51 :
52 10 : integer function alloc_eos_handle(ierr) result(handle)
53 : integer, intent(out) :: ierr ! 0 means AOK.
54 : character (len=0) :: inlist
55 10 : handle = alloc_eos_handle_using_inlist(inlist, ierr)
56 5 : end function alloc_eos_handle
57 :
58 19 : integer function alloc_eos_handle_using_inlist(inlist,ierr) result(handle)
59 : use eos_def, only:do_alloc_eos
60 : use eos_ctrls_io, only:read_namelist
61 : character (len=*), intent(in) :: inlist ! empty means just use defaults.
62 : integer, intent(out) :: ierr ! 0 means AOK.
63 : ierr = 0
64 19 : handle = do_alloc_eos(ierr)
65 19 : if (ierr /= 0) return
66 19 : call read_namelist(handle, inlist, ierr)
67 19 : end function alloc_eos_handle_using_inlist
68 :
69 5 : subroutine free_eos_handle(handle)
70 : ! frees the handle and all associated data
71 19 : use eos_def,only:do_free_eos_handle
72 : integer, intent(in) :: handle
73 5 : call do_free_eos_handle(handle)
74 5 : end subroutine free_eos_handle
75 :
76 :
77 3 : subroutine eos_ptr(handle,rq,ierr)
78 5 : use eos_def,only:EoS_General_Info,get_eos_ptr
79 : integer, intent(in) :: handle ! from alloc_eos_handle
80 : type (EoS_General_Info), pointer :: rq
81 : integer, intent(out):: ierr
82 3 : call get_eos_ptr(handle,rq,ierr)
83 3 : end subroutine eos_ptr
84 :
85 :
86 : ! as a convenience
87 :
88 13080 : elemental real(dp) function Radiation_Pressure(T)
89 3 : use const_def, only: crad
90 : real(dp), intent(in) :: T
91 13080 : Radiation_Pressure = crad*T*T*T*T/3d0
92 13080 : end function Radiation_Pressure
93 :
94 0 : elemental real(dp) function Radiation_Energy(T)
95 : use const_def, only: crad
96 : real(dp), intent(in) :: T
97 0 : Radiation_Energy = crad*T*T*T*T
98 0 : end function Radiation_Energy
99 :
100 :
101 : ! eos evaluation
102 :
103 : ! you can call these routines after you've allocated a handle.
104 : ! NOTE: the information referenced via the handle is read-only,
105 : ! so you can do multiple evaluations in parallel using the same handle.
106 :
107 :
108 195688 : subroutine eosDT_get( &
109 195688 : handle, species, chem_id, net_iso, xa, &
110 : Rho, logRho, T, logT, &
111 195688 : res, d_dlnd, d_dlnT, d_dxa, ierr)
112 : use eos_def
113 : use eosDT_eval, only: Get_eosDT_Results
114 : use chem_lib, only: basic_composition_info
115 : integer, intent(in) :: handle ! eos handle; from star, pass s% eos_handle
116 : integer, intent(in) :: species ! number of species
117 : integer, pointer :: chem_id(:) ! maps species to chem id
118 : integer, pointer :: net_iso(:) ! maps chem id to species number
119 : real(dp), intent(in) :: xa(:) ! mass fractions
120 : real(dp), intent(in) :: Rho, logRho ! the density
121 : real(dp), intent(in) :: T, logT ! the temperature
122 : real(dp), intent(inout) :: res(:) ! (num_eos_basic_results)
123 : real(dp), intent(inout) :: d_dlnd(:) ! (num_eos_basic_results)
124 : real(dp), intent(inout) :: d_dlnT(:) ! (num_eos_basic_results)
125 : real(dp), intent(inout) :: d_dxa(:,:) ! (num_eos_d_dxa_results,species)
126 : integer, intent(out) :: ierr ! 0 means AOK.
127 195688 : real(dp), allocatable :: d_dxa_eos(:,:) ! eos internally returns derivs of all quantities
128 : type (EoS_General_Info), pointer :: rq
129 195688 : real(dp) :: X, Y, Z, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
130 195688 : call get_eos_ptr(handle,rq,ierr)
131 195688 : if (ierr /= 0) then
132 0 : write(*,*) 'invalid handle for eos_get -- did you call alloc_eos_handle?'
133 : return
134 : end if
135 : call basic_composition_info( &
136 : species, chem_id, xa, X, Y, Z, &
137 195688 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
138 195688 : allocate(d_dxa_eos(num_eos_basic_results, species))
139 : call Get_eosDT_Results( &
140 : rq, Z, X, abar, zbar, &
141 : species, chem_id, net_iso, xa, &
142 : Rho, logRho, T, logT, &
143 195688 : res, d_dlnd, d_dlnT, d_dxa_eos, ierr)
144 : ! only return 1st two d_dxa results (lnE and lnPgas) to star
145 5015710 : d_dxa(1:num_eos_d_dxa_results,1:species) = d_dxa_eos(1:num_eos_d_dxa_results, 1:species)
146 391376 : end subroutine eosDT_get
147 :
148 :
149 0 : subroutine eosDT_get_component( &
150 : handle, which_eos, &
151 0 : species, chem_id, net_iso, xa, &
152 : Rho, log10Rho, T, log10T, &
153 0 : res, d_dlnRho_const_T, d_dlnT_const_Rho, d_dxa_const_TRho, &
154 : ierr)
155 :
156 195688 : use chem_lib, only: basic_composition_info
157 : use eos_def
158 : use eosDT_eval, only: Test_one_eosDT_component
159 :
160 : ! INPUT
161 :
162 : integer, intent(in) :: handle ! eos handle; from star, pass s% eos_handle
163 : integer, intent(in) :: which_eos ! see eos_def: i_eos_<component>
164 : integer, intent(in) :: species ! number of species
165 : integer, pointer :: chem_id(:) ! maps species to chem id
166 : ! index from 1 to species
167 : ! value is between 1 and num_chem_isos
168 : integer, pointer :: net_iso(:) ! maps chem id to species number
169 : ! index from 1 to num_chem_isos (defined in chem_def)
170 : ! value is 0 if the iso is not in the current net
171 : ! else is value between 1 and number of species in current net
172 : real(dp), intent(in) :: xa(:) ! mass fractions
173 :
174 : real(dp), intent(in) :: Rho, log10Rho ! the density
175 : ! provide both if you have them. else pass one and set the other to arg_not_provided
176 : ! "arg_not_provided" is defined in mesa const_def
177 :
178 : real(dp), intent(in) :: T, log10T ! the temperature
179 : ! provide both if you have them. else pass one and set the other to arg_not_provided
180 :
181 : ! OUTPUT
182 :
183 : real(dp), intent(inout) :: res(:) ! (num_eos_basic_results)
184 : ! partial derivatives of the basic results wrt lnd and lnT
185 :
186 : real(dp), intent(inout) :: d_dlnRho_const_T(:) ! (num_eos_basic_results)
187 : ! d_dlnRho_const_T(i) = d(res(i))/dlnd|T,X where X = composition
188 : real(dp), intent(inout) :: d_dlnT_const_Rho(:) ! (num_eos_basic_results)
189 : ! d_dlnT(i) = d(res(i))/dlnT|Rho,X where X = composition
190 : real(dp), intent(inout) :: d_dxa_const_TRho(:,:) ! (num_eos_d_dxa_results,species)
191 :
192 : integer, intent(out) :: ierr ! 0 means AOK.
193 :
194 0 : real(dp), allocatable :: d_dxa_eos(:,:) ! eos internally returns derivs of all quantities
195 :
196 : type (EoS_General_Info), pointer :: rq
197 :
198 0 : real(dp) :: X, Y, Z, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
199 :
200 0 : call get_eos_ptr(handle,rq,ierr)
201 0 : if (ierr /= 0) then
202 0 : write(*,*) 'invalid handle for eos_get -- did you call alloc_eos_handle?'
203 : return
204 : end if
205 :
206 : call basic_composition_info( &
207 : species, chem_id, xa, X, Y, Z, &
208 0 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
209 :
210 0 : allocate(d_dxa_eos(num_eos_basic_results,species))
211 :
212 : call Test_one_eosDT_component( &
213 : rq, which_eos, Z, X, abar, zbar, &
214 : species, chem_id, net_iso, xa, &
215 : Rho, log10Rho, T, log10T, &
216 0 : res, d_dlnRho_const_T, d_dlnT_const_Rho, d_dxa_eos, ierr)
217 :
218 : ! only return 1st two d_dxa results (lnE and lnPgas)
219 0 : d_dxa_const_TRho(1:num_eos_d_dxa_results,1:species) = d_dxa_eos(1:num_eos_d_dxa_results, 1:species)
220 :
221 0 : end subroutine eosDT_get_component
222 :
223 :
224 0 : subroutine helmeos2_eval( &
225 : T, logT, Rho, logRho, abar, zbar, &
226 0 : coulomb_temp_cut, coulomb_den_cut, helm_res, &
227 : clip_to_table_boundaries, include_radiation, &
228 : include_elec_pos, &
229 : off_table, ierr)
230 0 : use helm
231 : real(dp), intent(in) :: T, logT, Rho, logRho, abar, zbar, &
232 : coulomb_temp_cut, coulomb_den_cut
233 : real(dp), intent(inout) :: helm_res(:) ! (num_helm_results)
234 : logical, intent(in) :: clip_to_table_boundaries, include_radiation, &
235 : include_elec_pos
236 : logical, intent(out) :: off_table
237 : integer, intent(out) :: ierr ! 0 means AOK.
238 : call helmeos2( &
239 : T, logT, Rho, logRho, abar, zbar, coulomb_temp_cut, coulomb_den_cut, &
240 : helm_res, clip_to_table_boundaries, include_radiation, include_elec_pos, &
241 0 : off_table, ierr)
242 0 : end subroutine helmeos2_eval
243 :
244 :
245 : ! the following routine uses gas pressure and temperature as input variables
246 12150 : subroutine eosPT_get( &
247 : handle, &
248 12150 : species, chem_id, net_iso, xa, &
249 : Pgas, log10Pgas, T, log10T, &
250 : Rho, log10Rho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas, &
251 12150 : res, d_dlnRho_const_T, d_dlnT_const_Rho, &
252 12150 : d_dxa_const_TRho, ierr)
253 :
254 0 : use chem_lib, only: basic_composition_info
255 : use eos_def
256 : use eosPT_eval, only: Get_eosPT_Results
257 :
258 : ! INPUT
259 :
260 : integer, intent(in) :: handle ! eos handle; from star, pass s% eos_handle
261 :
262 : integer, intent(in) :: species ! number of species
263 : integer, pointer :: chem_id(:) ! maps species to chem id
264 : ! index from 1 to species
265 : ! value is between 1 and num_chem_isos
266 : integer, pointer :: net_iso(:) ! maps chem id to species number
267 : ! index from 1 to num_chem_isos (defined in chem_def)
268 : ! value is 0 if the iso is not in the current net
269 : ! else is value between 1 and number of species in current net
270 : real(dp), intent(in) :: xa(:) ! mass fractions
271 :
272 : real(dp), intent(in) :: Pgas, log10Pgas ! the gas pressure
273 : ! provide both if you have them. else pass one and set the other to arg_not_provided
274 : ! "arg_not_provided" is defined in mesa const_def
275 :
276 : real(dp), intent(in) :: T, log10T ! the temperature
277 : ! provide both if you have them. else pass one and set the other to arg_not_provided
278 :
279 : ! OUTPUT
280 :
281 : real(dp), intent(out) :: Rho, log10Rho ! density
282 : real(dp), intent(out) :: dlnRho_dlnPgas_const_T
283 : real(dp), intent(out) :: dlnRho_dlnT_const_Pgas
284 :
285 : real(dp), intent(inout) :: res(:) ! (num_eos_basic_results)
286 :
287 : ! partial derivatives of the basic results
288 :
289 : real(dp), intent(inout) :: d_dlnRho_const_T(:) ! (num_eos_basic_results)
290 : ! d_dlnRho_const_T(i) = d(res(i))/dlnd|T,X where X = composition
291 : real(dp), intent(inout) :: d_dlnT_const_Rho(:) ! (num_eos_basic_results)
292 : ! d_dlnT_const_Rho(i) = d(res(i))/dlnT|Rho,X where X = composition
293 : real(dp), intent(inout) :: d_dxa_const_TRho(:,:) ! (num_eos_d_dxa_results, species)
294 : ! d_dxa_const_TRho(i) = d(res(i))/X|T,Rho,X where X = composition
295 :
296 12150 : real(dp), allocatable :: d_dxa_eos(:,:) ! eos internally returns derivs of all quantities
297 :
298 : integer, intent(out) :: ierr ! 0 means AOK.
299 :
300 : type (EoS_General_Info), pointer :: rq
301 :
302 12150 : real(dp) :: X, Y, Z, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
303 :
304 12150 : call get_eos_ptr(handle,rq,ierr)
305 12150 : if (ierr /= 0) then
306 0 : write(*,*) 'invalid handle for eos -- did you call alloc_eos_handle?'
307 : return
308 : end if
309 :
310 : call basic_composition_info( &
311 : species, chem_id, xa, X, Y, Z, &
312 12150 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
313 :
314 12150 : allocate(d_dxa_eos(num_eos_basic_results, species))
315 :
316 : call Get_eosPT_Results( &
317 : rq, Z, X, abar, zbar, &
318 : species, chem_id, net_iso, xa, &
319 : Pgas, log10Pgas, T, log10T, &
320 : Rho, log10Rho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas, &
321 : res, d_dlnRho_const_T, d_dlnT_const_Rho, d_dxa_eos, &
322 12150 : ierr)
323 : ! only return 1st two d_dxa results (lnE and lnPgas) to star
324 516360 : d_dxa_const_TRho(1:num_eos_d_dxa_results,1:species) = d_dxa_eos(1:num_eos_d_dxa_results, 1:species)
325 :
326 24300 : end subroutine eosPT_get
327 :
328 :
329 : ! gamma law eos routines -- ignores radiation (e.g.., P = Pgas only)
330 :
331 :
332 0 : subroutine eos_gamma_DP_get_ET( &
333 : abar, rho, P, gamma, energy, T, ierr)
334 12150 : use const_def, only: avo, kerg
335 : real(dp), intent(in) :: abar, rho, P, gamma
336 : real(dp), intent(out) :: energy, T
337 : integer, intent(out) :: ierr
338 0 : ierr = 0
339 0 : energy = (P/rho)/(gamma - 1d0)
340 0 : T = (gamma - 1d0)*energy*abar/(avo*kerg)
341 0 : end subroutine eos_gamma_DP_get_ET
342 :
343 :
344 0 : subroutine eos_gamma_DE_get_PT( &
345 : abar, rho, energy, gamma, P, T, ierr)
346 : use const_def, only: avo, kerg
347 : real(dp), intent(in) :: abar, rho, energy, gamma
348 : real(dp), intent(out) :: P, T
349 : integer, intent(out) :: ierr
350 0 : ierr = 0
351 0 : P = (gamma - 1d0)*energy*rho
352 0 : T = (gamma - 1d0)*energy*abar/(avo*kerg)
353 0 : end subroutine eos_gamma_DE_get_PT
354 :
355 :
356 0 : subroutine eos_gamma_DT_get_P_energy( &
357 : abar, rho, T, gamma, P, energy, ierr)
358 : use const_def, only: avo, kerg
359 : real(dp), intent(in) :: abar, rho, T, gamma
360 : real(dp), intent(out) :: P, energy
361 : integer, intent(out) :: ierr
362 0 : ierr = 0
363 0 : P = avo*kerg*rho*T/abar
364 0 : energy = (P/rho)/(gamma - 1d0)
365 0 : end subroutine eos_gamma_DT_get_P_energy
366 :
367 :
368 0 : subroutine eos_gamma_PRho_get_T_energy( &
369 : abar, P, rho, gamma, T, energy, ierr)
370 : use const_def, only: avo, kerg
371 : real(dp), intent(in) :: abar, P, rho, gamma
372 : real(dp), intent(out) :: T, energy
373 : integer, intent(out) :: ierr
374 0 : ierr = 0
375 0 : energy = (P/rho)/(gamma - 1d0)
376 0 : T = (gamma - 1d0)*energy*abar/(avo*kerg)
377 0 : end subroutine eos_gamma_PRho_get_T_energy
378 :
379 :
380 0 : subroutine eos_gamma_PT_get_rho_energy( &
381 : abar, P, T, gamma, rho, energy, ierr)
382 : use const_def, only: avo, kerg
383 : real(dp), intent(in) :: abar, P, T, gamma
384 : real(dp), intent(out) :: rho, energy
385 : integer, intent(out) :: ierr
386 0 : ierr = 0
387 0 : rho = (P/T)*abar/(avo*kerg)
388 0 : energy = (P/rho)/(gamma - 1d0)
389 0 : end subroutine eos_gamma_PT_get_rho_energy
390 :
391 :
392 0 : subroutine eos_gamma_DE_get( &
393 : handle, abar, energy, log10E, rho, log10Rho, gamma, &
394 0 : T, log10T, res, d_dlnRho_const_T, d_dlnT_const_Rho, &
395 : dlnT_dlnE_c_Rho, dlnT_dlnd_c_E, &
396 : dlnPgas_dlnE_c_Rho, dlnPgas_dlnd_c_E, ierr)
397 : use eos_def
398 : use eosDE_eval, only: Get_eos_gamma_DE_Results
399 : integer, intent(in) :: handle
400 : real(dp), intent(in) :: abar, energy, log10E, Rho, log10Rho, gamma
401 : real(dp), intent(out) :: T, log10T
402 : real(dp), intent(inout), dimension(:) :: &
403 : res, d_dlnRho_const_T, d_dlnT_const_Rho
404 : real(dp), intent(out) :: &
405 : dlnT_dlnE_c_Rho, dlnT_dlnd_c_E, &
406 : dlnPgas_dlnE_c_Rho, dlnPgas_dlnd_c_E
407 : integer, intent(out) :: ierr
408 : type (EoS_General_Info), pointer :: rq
409 0 : call get_eos_ptr(handle,rq,ierr)
410 0 : if (ierr /= 0) then
411 0 : write(*,*) 'invalid handle for eos -- did you call alloc_eos_handle?'
412 0 : return
413 : end if
414 : ! require positive density and energy
415 0 : if ((rho <= 0) .or. (energy <= 0)) then
416 0 : ierr = -1
417 0 : return
418 : end if
419 : call Get_eos_gamma_DE_Results( &
420 : rq, abar, energy, log10E, rho, log10Rho, gamma, &
421 : T, log10T, res, d_dlnRho_const_T, d_dlnT_const_Rho, &
422 : dlnT_dlnE_c_Rho, dlnT_dlnd_c_E, &
423 0 : dlnPgas_dlnE_c_Rho, dlnPgas_dlnd_c_E, ierr)
424 0 : end subroutine eos_gamma_DE_get
425 :
426 :
427 0 : subroutine eos_gamma_PT_get( &
428 : handle, abar, P, log10P, T, log10T, gamma, &
429 0 : rho, log10Rho, res, d_dlnRho_const_T, d_dlnT_const_Rho, &
430 : ierr)
431 0 : use eos_def
432 : use eosDE_eval, only: Get_eos_gamma_DE_Results
433 : use math_lib
434 : integer, intent(in) :: handle
435 : real(dp), intent(in) :: abar, P, log10P, T, log10T, gamma
436 : real(dp), intent(out) :: rho, log10Rho
437 : real(dp), intent(inout), dimension(:) :: &
438 : res, d_dlnRho_const_T, d_dlnT_const_Rho
439 : integer, intent(out) :: ierr
440 : type (EoS_General_Info), pointer :: rq
441 0 : real(dp) :: energy, temp, log10temp, &
442 0 : dlnT_dlnE_c_Rho, dlnT_dlnd_c_E, &
443 0 : dlnPgas_dlnE_c_Rho, dlnPgas_dlnd_c_E
444 0 : call get_eos_ptr(handle,rq,ierr)
445 0 : if (ierr /= 0) then
446 0 : write(*,*) 'invalid handle for eos -- did you call alloc_eos_handle?'
447 0 : return
448 : end if
449 : ! require positive pressure and temperature
450 0 : if ((P <= 0) .or. (T <= 0)) then
451 0 : ierr = -1
452 0 : return
453 : end if
454 : call eos_gamma_PT_get_rho_energy( &
455 0 : abar, P, T, gamma, rho, energy, ierr)
456 0 : log10Rho = log10(rho)
457 : if (ierr /= 0) return
458 : call eos_gamma_DE_get( &
459 : handle, abar, energy, log10(energy), rho, log10Rho, gamma, &
460 : temp, log10temp, res, d_dlnRho_const_T, d_dlnT_const_Rho, &
461 : dlnT_dlnE_c_Rho, dlnT_dlnd_c_E, &
462 0 : dlnPgas_dlnE_c_Rho, dlnPgas_dlnd_c_E, ierr)
463 0 : end subroutine eos_gamma_PT_get
464 :
465 :
466 0 : subroutine eos_gamma_DT_get( &
467 : handle, abar, rho, log10Rho, T, log10T, gamma, &
468 0 : res, d_dlnRho_const_T, d_dlnT_const_Rho, &
469 : Pgas, Prad, energy, entropy, ierr)
470 0 : use eos_def
471 : use eosDE_eval, only: Get_eos_gamma_DE_Results
472 : use math_lib
473 : integer, intent(in) :: handle
474 : real(dp), intent(in) :: abar, rho, log10Rho, T, log10T, gamma
475 : real(dp), intent(inout), dimension(:) :: &
476 : res, d_dlnRho_const_T, d_dlnT_const_Rho
477 : real(dp), intent(out) :: Pgas, Prad, energy, entropy
478 : integer, intent(out) :: ierr
479 : type (EoS_General_Info), pointer :: rq
480 0 : real(dp) :: P, temp, log10temp, &
481 0 : dlnT_dlnE_c_Rho, dlnT_dlnd_c_E, &
482 0 : dlnPgas_dlnE_c_Rho, dlnPgas_dlnd_c_E
483 0 : call get_eos_ptr(handle,rq,ierr)
484 0 : if (ierr /= 0) then
485 0 : write(*,*) 'invalid handle for eos -- did you call alloc_eos_handle?'
486 0 : return
487 : end if
488 : ! require positive density and temperature
489 0 : if ((rho <= 0) .or. (T <= 0)) then
490 0 : ierr = -1
491 0 : return
492 : end if
493 : call eos_gamma_DT_get_P_energy( &
494 0 : abar, rho, T, gamma, P, energy, ierr)
495 : if (ierr /= 0) return
496 : call eos_gamma_DE_get( &
497 : handle, abar, energy, log10(energy), rho, log10Rho, gamma, &
498 : temp, log10temp, res, d_dlnRho_const_T, d_dlnT_const_Rho, &
499 : dlnT_dlnE_c_Rho, dlnT_dlnd_c_E, &
500 0 : dlnPgas_dlnE_c_Rho, dlnPgas_dlnd_c_E, ierr)
501 0 : Pgas = exp(res(i_lnPgas))
502 0 : Prad = crad*T*T*T*T/3d0
503 0 : energy = exp(res(i_lnE))
504 0 : entropy = exp(res(i_lnS))
505 0 : end subroutine eos_gamma_DT_get
506 :
507 :
508 : ! misc
509 :
510 :
511 112 : subroutine eos_fermi_dirac_integral(dk, eta, theta, fd, fdeta, fdtheta)
512 : !..from Frank Timmes' site, http://www.cococubed.com/code_pages/fermi_dirac.shtml
513 : !..this routine computes the fermi-dirac integrals of
514 : !..index dk, with degeneracy parameter eta and relativity parameter theta.
515 : !..input is dk the real(dp) index of the fermi-dirac function,
516 : !..eta the degeneracy parameter, and theta the relativity parameter.
517 : !..theta = (k * T)/(mass_electron * c^2), k = Boltzmann const.
518 : !..the output is fd is computed by applying three 10-point
519 : !..gauss-legendre and one 10-point gauss-laguerre rules over
520 : !..four appropriate subintervals. the derivative with respect to eta is
521 : !..output in fdeta, and the derivative with respct to theta is in fdtheta.
522 : !..within each subinterval the fd kernel.
523 : !..
524 : !..this routine delivers 14 significant figure accuracy
525 : !..
526 : !..reference: j.m. aparicio, apjs 117, 632 1998
527 0 : use gauss_fermi, only: dfermi
528 :
529 : real(dp), intent(in) :: dk
530 : real(dp), intent(in) :: eta
531 : real(dp), intent(in) :: theta
532 : real(dp), intent(out) :: fd
533 : real(dp), intent(out) :: fdeta
534 : real(dp), intent(out) :: fdtheta
535 112 : call dfermi(dk, eta, theta, fd, fdeta, fdtheta)
536 112 : end subroutine eos_fermi_dirac_integral
537 :
538 :
539 0 : subroutine eos_get_helm_results( &
540 : X, abar, zbar, Rho, log10Rho, T, log10T, &
541 : coulomb_temp_cut, coulomb_den_cut, &
542 : include_radiation, include_elec_pos, &
543 0 : res, ierr)
544 : ! direct call to the helm eos.
545 : ! returns much more info than the standard
546 :
547 112 : use eos_def
548 : use eosDT_eval, only: Get_HELM_Results
549 :
550 : ! INPUT
551 :
552 : real(dp), intent(in) :: X ! the hydrogen mass fraction
553 :
554 : real(dp), intent(in) :: abar
555 : ! mean atomic number (nucleons per nucleus; grams per mole)
556 : real(dp), intent(in) :: zbar ! mean charge per nucleus
557 :
558 : real(dp), intent(in) :: Rho, log10Rho ! the density
559 : ! provide both if you have them.
560 : ! else pass one and set the other to arg_not_provided
561 :
562 : real(dp), intent(in) :: T, log10T ! the temperature
563 : ! provide both if you have them.
564 : ! else pass one and set the other to arg_not_provided
565 :
566 : real(dp), intent(in) :: coulomb_temp_cut, coulomb_den_cut
567 :
568 : logical, intent(in) :: include_radiation, include_elec_pos
569 :
570 : ! OUTPUT
571 :
572 : real(dp), intent(inout) :: res(:) ! (num_helm_results)
573 : ! array to hold the results
574 : integer, intent(out) :: ierr ! 0 means AOK.
575 :
576 : logical :: off_table
577 :
578 : call Get_HELM_Results( &
579 : abar, zbar, Rho, log10Rho, T, log10T, &
580 : coulomb_temp_cut, coulomb_den_cut, &
581 : include_radiation, include_elec_pos, &
582 0 : res, off_table, ierr)
583 :
584 0 : end subroutine eos_get_helm_results
585 :
586 :
587 : !subroutine eos_convert_helm_results( &
588 : ! helm_res, Z, X, abar, zbar, Rho, T, res, &
589 : ! d_dlnRho_const_T, d_dlnT_const_Rho, ierr)
590 0 : subroutine eos_convert_helm_results( &
591 0 : helm_res, Z, X, abar, zbar, Rho, T, basic_flag, res, &
592 0 : d_dlnRho_const_T, d_dlnT_const_Rho, &
593 0 : d_dabar_const_TRho, d_dzbar_const_TRho, ierr)
594 0 : use eos_def
595 : use eos_helm_eval, only: do_convert_helm_results
596 : real(dp), intent(in) :: helm_res(:) ! (num_helm_results)
597 : real(dp), intent(in) :: Z, X, abar, zbar, Rho, T
598 : logical, intent(in) :: basic_flag ! if true, then only want basic results
599 : real(dp), intent(inout) :: res(:) ! (num_eos_basic_results)
600 : real(dp), intent(inout) :: d_dlnRho_const_T(:) ! (num_eos_basic_results)
601 : real(dp), intent(inout) :: d_dlnT_const_Rho(:) ! (num_eos_basic_results)
602 : real(dp), intent(inout) :: d_dabar_const_TRho(:) ! (num_eos_basic_results)
603 : real(dp), intent(inout) :: d_dzbar_const_TRho(:) ! (num_eos_basic_results)
604 : !real(dp), intent(inout), dimension(:) :: d2_dlnd2, d2_dlnd_dlnT, d2_dlnT2
605 : integer, intent(out) :: ierr
606 0 : d_dabar_const_TRho = 0
607 0 : d_dzbar_const_TRho = 0
608 : call do_convert_helm_results( &
609 : helm_res, Z, abar, zbar, Rho, T, &
610 : res, d_dlnRho_const_T, d_dlnT_const_Rho, &
611 : d_dabar_const_TRho, d_dzbar_const_TRho, &
612 0 : ierr)
613 0 : end subroutine eos_convert_helm_results
614 :
615 :
616 : ! eosDT search routines. these use num_lib safe_root to find T or Rho.
617 :
618 1476 : subroutine eosDT_get_T( &
619 : handle, &
620 1476 : species, chem_id, net_iso, xa, &
621 : logRho, which_other, other_value, &
622 : logT_tol, other_tol, max_iter, logT_guess, &
623 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
624 1476 : logT_result, res, d_dlnRho_const_T, d_dlnT_const_Rho, &
625 1476 : d_dxa_const_TRho, eos_calls, ierr)
626 :
627 : ! finds log10 T given values for density and 'other', and initial guess for temperature.
628 : ! does up to max_iter attempts until logT changes by less than tol.
629 :
630 : ! 'other' can be any of the basic result variables for the eos
631 : ! specify 'which_other' by means of the definitions in eos_def (e.g., i_lnE)
632 :
633 0 : use chem_lib, only: basic_composition_info
634 : use eos_def
635 : use eosDT_eval, only : get_T
636 :
637 : integer, intent(in) :: handle ! eos handle; from star, pass s% eos_handle
638 :
639 : integer, intent(in) :: species ! number of species
640 : integer, pointer :: chem_id(:) ! maps species to chem id
641 : ! index from 1 to species
642 : ! value is between 1 and num_chem_isos
643 : integer, pointer :: net_iso(:) ! maps chem id to species number
644 : ! index from 1 to num_chem_isos (defined in chem_def)
645 : ! value is 0 if the iso is not in the current net
646 : ! else is value between 1 and number of species in current net
647 : real(dp), intent(in) :: xa(:) ! mass fractions
648 :
649 : real(dp), intent(in) :: logRho ! log10 of density
650 : integer, intent(in) :: which_other ! from eos_def. e.g., i_lnE
651 : real(dp), intent(in) :: other_value ! desired value for the other variable
652 : real(dp), intent(in) :: other_tol
653 :
654 : real(dp), intent(in) :: logT_tol
655 : integer, intent(in) :: max_iter ! max number of iterations
656 :
657 : real(dp), intent(in) :: logT_guess ! log10 of temperature
658 : real(dp), intent(in) :: logT_bnd1, logT_bnd2 ! bounds for logT
659 : ! if don't know bounds, just set to arg_not_provided (defined in const_def)
660 : real(dp), intent(in) :: other_at_bnd1, other_at_bnd2 ! values at bounds
661 : ! if don't know these values, just set to arg_not_provided (defined in const_def)
662 :
663 : real(dp), intent(inout) :: logT_result ! log10 of temperature
664 : real(dp), intent(inout) :: res(:) ! (num_eos_basic_results)
665 : real(dp), intent(inout) :: d_dlnRho_const_T(:) ! (num_eos_basic_results)
666 : real(dp), intent(inout) :: d_dlnT_const_Rho(:) ! (num_eos_basic_results)
667 : real(dp), intent(inout) :: d_dxa_const_TRho(:,:) ! (num_eos_d_dxa_results, species)
668 : real(dp), allocatable :: d_dxa_eos(:,:) ! eos internally returns derivs of all quantities
669 :
670 : integer, intent(out) :: eos_calls
671 : integer, intent(out) :: ierr ! 0 means AOK.
672 :
673 : ! compute composition info
674 : real(dp) :: Y, Z, X, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
675 :
676 : call basic_composition_info( &
677 : species, chem_id, xa, X, Y, Z, &
678 1476 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
679 :
680 1476 : allocate(d_dxa_eos(num_eos_basic_results, species))
681 :
682 : call get_T( &
683 : handle, Z, X, abar, zbar, &
684 : species, chem_id, net_iso, xa, &
685 : logRho, which_other, other_value, &
686 : logT_tol, other_tol, max_iter, logT_guess, &
687 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
688 : logT_result, res, d_dlnRho_const_T, d_dlnT_const_Rho, &
689 1476 : d_dxa_eos, eos_calls, ierr)
690 : ! only return 1st two d_dxa results (lnE and lnPgas) to star
691 36888 : d_dxa_const_TRho(1:num_eos_d_dxa_results,1:species) = d_dxa_eos(1:num_eos_d_dxa_results, 1:species)
692 :
693 1476 : deallocate(d_dxa_eos)
694 :
695 1476 : end subroutine eosDT_get_T
696 :
697 :
698 2 : subroutine eosDT_get_Rho( &
699 : handle, &
700 2 : species, chem_id, net_iso, xa, &
701 : logT, which_other, other_value, &
702 : logRho_tol, other_tol, max_iter, logRho_guess, &
703 : logRho_bnd1, logRho_bnd2, other_at_bnd1, other_at_bnd2, &
704 2 : logRho_result, res, d_dlnRho_const_T, d_dlnT_const_Rho, &
705 2 : d_dxa_const_TRho, eos_calls, ierr)
706 :
707 : ! finds log10 Rho given values for temperature and 'other', and initial guess for density.
708 : ! does up to max_iter attempts until logRho changes by less than tol.
709 :
710 : ! 'other' can be any of the basic result variables for the eos
711 : ! specify 'which_other' by means of the definitions in eos_def (e.g., i_lnE)
712 :
713 1476 : use chem_lib, only: basic_composition_info
714 : use eos_def
715 : use eosDT_eval, only : get_Rho
716 :
717 : integer, intent(in) :: handle ! eos handle; from star, pass s% eos_handle
718 :
719 : integer, intent(in) :: species ! number of species
720 : integer, pointer :: chem_id(:) ! maps species to chem id
721 : ! index from 1 to species
722 : ! value is between 1 and num_chem_isos
723 : integer, pointer :: net_iso(:) ! maps chem id to species number
724 : ! index from 1 to num_chem_isos (defined in chem_def)
725 : ! value is 0 if the iso is not in the current net
726 : ! else is value between 1 and number of species in current net
727 : real(dp), intent(in) :: xa(:) ! mass fractions
728 :
729 : real(dp), intent(in) :: logT ! log10 of temperature
730 :
731 : integer, intent(in) :: which_other ! from eos_def. e.g., i_lnE
732 : real(dp), intent(in) :: other_value ! desired value for the other variable
733 : real(dp), intent(in) :: other_tol
734 :
735 : real(dp), intent(in) :: logRho_tol
736 :
737 : integer, intent(in) :: max_iter ! max number of Newton iterations
738 :
739 : real(dp), intent(in) :: logRho_guess ! log10 of density
740 : real(dp), intent(in) :: logRho_bnd1, logRho_bnd2 ! bounds for logRho
741 : ! if don't know bounds, just set to arg_not_provided (defined in const_def)
742 : real(dp), intent(in) :: other_at_bnd1, other_at_bnd2 ! values at bounds
743 : ! if don't know these values, just set to arg_not_provided (defined in const_def)
744 :
745 : real(dp), intent(out) :: logRho_result ! log10 of density
746 :
747 : real(dp), intent(inout) :: res(:) ! (num_eos_basic_results)
748 : real(dp), intent(inout) :: d_dlnRho_const_T(:) ! (num_eos_basic_results)
749 : real(dp), intent(inout) :: d_dlnT_const_Rho(:) ! (num_eos_basic_results)
750 : real(dp), intent(inout) :: d_dxa_const_TRho(:,:) ! (num_eos_d_dxa_results, species)
751 : real(dp), allocatable :: d_dxa_eos(:,:) ! eos internally returns derivs of all quantities
752 :
753 : integer, intent(out) :: eos_calls
754 : integer, intent(out) :: ierr ! 0 means AOK.
755 :
756 : ! compute composition info
757 : real(dp) :: Y, Z, X, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
758 :
759 : call basic_composition_info( &
760 : species, chem_id, xa, X, Y, Z, &
761 2 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
762 :
763 2 : allocate(d_dxa_eos(num_eos_basic_results, species))
764 :
765 : call get_Rho( &
766 : handle, Z, X, abar, zbar, &
767 : species, chem_id, net_iso, xa, &
768 : logT, which_other, other_value, &
769 : logRho_tol, other_tol, max_iter, logRho_guess, &
770 : logRho_bnd1, logRho_bnd2, other_at_bnd1, other_at_bnd2, &
771 : logRho_result, res, d_dlnRho_const_T, d_dlnT_const_Rho, &
772 2 : d_dxa_eos, eos_calls, ierr)
773 : ! only return 1st two d_dxa results (lnE and lnPgas) to star
774 44 : d_dxa_const_TRho(1:num_eos_d_dxa_results,1:species) = d_dxa_eos(1:num_eos_d_dxa_results, 1:species)
775 :
776 2 : deallocate(d_dxa_eos)
777 :
778 2 : end subroutine eosDT_get_Rho
779 :
780 :
781 : ! eosPT search routines. these use num_lib safe_root to find T or Pgas.
782 :
783 0 : subroutine eosPT_get_T( &
784 : handle, &
785 0 : species, chem_id, net_iso, xa, &
786 : logPgas, which_other, other_value, &
787 : logT_tol, other_tol, max_iter, logT_guess, &
788 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
789 : logT_result, Rho, log10Rho, &
790 : dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas, &
791 0 : res, d_dlnRho_const_T, d_dlnT_const_Rho, &
792 0 : d_dxa_const_TRho, &
793 : eos_calls, ierr)
794 :
795 : ! finds log10 T given values for gas pressure and 'other',
796 : ! and initial guess for temperature.
797 : ! does up to max_iter attempts until logT changes by less than tol.
798 :
799 : ! 'other' can be any of the basic result variables for the eos
800 : ! specify 'which_other' by means of the definitions in eos_def (e.g., i_lnE)
801 :
802 2 : use chem_lib, only: basic_composition_info
803 : use eos_def
804 : use eosPT_eval, only : get_T
805 :
806 : integer, intent(in) :: handle ! eos handle; from star, pass s% eos_handle
807 :
808 : integer, intent(in) :: species ! number of species
809 : integer, pointer :: chem_id(:) ! maps species to chem id
810 : ! index from 1 to species
811 : ! value is between 1 and num_chem_isos
812 : integer, pointer :: net_iso(:) ! maps chem id to species number
813 : ! index from 1 to num_chem_isos (defined in chem_def)
814 : ! value is 0 if the iso is not in the current net
815 : ! else is value between 1 and number of species in current net
816 : real(dp), intent(in) :: xa(:) ! mass fractions
817 :
818 : real(dp), intent(in) :: logPgas ! log10 of gas pressure
819 : integer, intent(in) :: which_other ! from eos_def. e.g., i_lnE
820 : real(dp), intent(in) :: other_value ! desired value for the other variable
821 : real(dp), intent(in) :: other_tol
822 :
823 : real(dp), intent(in) :: logT_tol
824 : integer, intent(in) :: max_iter ! max number of iterations
825 :
826 : real(dp), intent(in) :: logT_guess ! log10 of temperature
827 : real(dp), intent(in) :: logT_bnd1, logT_bnd2 ! bounds for logT
828 : ! if don't know bounds, just set to arg_not_provided (defined in const_def)
829 : real(dp), intent(in) :: other_at_bnd1, other_at_bnd2 ! values at bounds
830 : ! if don't know these values, just set to arg_not_provided (defined in const_def)
831 :
832 : real(dp), intent(out) :: logT_result ! log10 of temperature
833 : real(dp), intent(out) :: Rho, log10Rho ! density
834 : real(dp), intent(out) :: dlnRho_dlnPgas_const_T
835 : real(dp), intent(out) :: dlnRho_dlnT_const_Pgas
836 :
837 : real(dp), intent(inout) :: res(:) ! (num_eos_basic_results)
838 : real(dp), intent(inout) :: d_dlnRho_const_T(:) ! (num_eos_basic_results)
839 : real(dp), intent(inout) :: d_dlnT_const_Rho(:) ! (num_eos_basic_results)
840 : real(dp), intent(inout) :: d_dxa_const_TRho(:,:) ! (num_eos_d_dxa_results, species)
841 : real(dp), allocatable :: d_dxa_eos(:,:) ! eos internally returns derivs of all quantities
842 :
843 : integer, intent(out) :: eos_calls
844 : integer, intent(out) :: ierr ! 0 means AOK.
845 :
846 : ! compute composition info
847 : real(dp) :: Y, Z, X, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
848 :
849 : call basic_composition_info( &
850 : species, chem_id, xa, X, Y, Z, &
851 0 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
852 :
853 0 : allocate(d_dxa_eos(num_eos_basic_results, species))
854 :
855 : call get_T( &
856 : handle, Z, X, abar, zbar, &
857 : species, chem_id, net_iso, xa, &
858 : logPgas, which_other, other_value, &
859 : logT_tol, other_tol, max_iter, logT_guess, &
860 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
861 : logT_result, Rho, log10Rho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas, &
862 : res, d_dlnRho_const_T, d_dlnT_const_Rho, d_dxa_eos, &
863 0 : eos_calls, ierr)
864 : ! only return 1st two d_dxa results (lnE and lnPgas) to star
865 0 : d_dxa_const_TRho(1:num_eos_d_dxa_results,1:species) = d_dxa_eos(1:num_eos_d_dxa_results, 1:species)
866 :
867 0 : deallocate(d_dxa_eos)
868 :
869 0 : end subroutine eosPT_get_T
870 :
871 :
872 0 : subroutine num_eos_files_loaded( &
873 : num_DT, num_FreeEOS)
874 0 : use eos_def
875 : integer, intent(out) :: &
876 : num_DT, num_FreeEOS
877 0 : num_DT = count(eosDT_XZ_loaded)
878 0 : num_FreeEOS = count(FreeEOS_XZ_loaded)
879 0 : end subroutine num_eos_files_loaded
880 :
881 :
882 0 : subroutine eos_get_control_namelist(handle, name, val, ierr)
883 0 : use eos_def
884 : use eos_ctrls_io, only: get_eos_controls
885 : integer, intent(in) :: handle ! eos handle; from star, pass s% eos_handle
886 : character(len=*),intent(in) :: name
887 : character(len=*),intent(out) :: val
888 : integer, intent(out) :: ierr
889 : type (EoS_General_Info), pointer :: rq
890 : ierr = 0
891 0 : call get_eos_ptr(handle,rq,ierr)
892 0 : if(ierr/=0) return
893 0 : call get_eos_controls(rq, name, val, ierr)
894 0 : end subroutine eos_get_control_namelist
895 :
896 0 : subroutine eos_set_control_namelist(handle, name, val, ierr)
897 0 : use eos_def
898 : use eos_ctrls_io, only: set_eos_controls
899 : integer, intent(in) :: handle ! eos handle; from star, pass s% eos_handle
900 : character(len=*),intent(in) :: name
901 : character(len=*),intent(in) :: val
902 : integer, intent(out) :: ierr
903 : type (EoS_General_Info), pointer :: rq
904 : ierr = 0
905 0 : call get_eos_ptr(handle,rq,ierr)
906 0 : if(ierr/=0) return
907 0 : call set_eos_controls(rq, name, val, ierr)
908 0 : end subroutine eos_set_control_namelist
909 :
910 : end module eos_lib
|
