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 : !..here is the tabular helmholtz free energy eos:
21 : !..
22 : !..routine helmeos computes the pressure, energy and entropy via tables
23 :
24 : module helm
25 : use const_def, only: dp
26 : use math_lib
27 :
28 : implicit none
29 :
30 : logical, parameter :: dbg = .false.
31 : !logical, parameter :: dbg = .true.
32 :
33 : private :: dbg
34 :
35 : contains
36 :
37 0 : subroutine helmeos2( &
38 : T, logT, Rho, logRho, abar_in, zbar_in, &
39 : coulomb_temp_cut, coulomb_den_cut, helm_res, &
40 : clip_to_table_boundaries, include_radiation, &
41 : include_elec_pos, off_table, ierr)
42 : use eos_def
43 : use const_def, only: pi, avo
44 : use utils_lib, only: is_bad
45 : real(dp), intent(in) :: T, logT, Rho, logRho
46 : real(dp), intent(in) :: abar_in, zbar_in
47 : real(dp), intent(in) :: coulomb_temp_cut, coulomb_den_cut
48 : real(dp), intent(inout) :: helm_res(num_helm_results)
49 : logical, intent(in) :: &
50 : clip_to_table_boundaries, include_radiation, &
51 : include_elec_pos
52 : logical, intent(out) :: off_table
53 : integer, intent(out) :: ierr
54 :
55 : logical :: skip_elec_pos
56 :
57 : include 'formats'
58 :
59 : ierr = 0
60 : off_table = .false.
61 :
62 0 : skip_elec_pos = .not. include_elec_pos
63 : call helmeos2aux( &
64 : T, logT, Rho, logRho, abar_in, zbar_in, &
65 : coulomb_temp_cut, coulomb_den_cut, helm_res, &
66 0 : clip_to_table_boundaries, include_radiation, (.not. include_elec_pos), off_table, ierr)
67 0 : if (off_table) return
68 :
69 0 : end subroutine helmeos2
70 :
71 :
72 0 : subroutine helmeos2aux( &
73 : temp_in, logtemp_in, den_in, logden_in, abar_in, zbar_in, &
74 : coulomb_temp_cut, coulomb_den_cut, helm_res, &
75 : clip_to_table_boundaries, include_radiation, must_skip_elec_pos, off_table, ierr)
76 :
77 0 : use eos_def
78 : use helm_polynomials
79 : use const_def, asol => crad
80 : use utils_lib, only: is_bad
81 :
82 :
83 : real(dp), intent(in) :: temp_in, logtemp_in, den_in, logden_in
84 : real(dp), intent(in) :: abar_in, zbar_in, coulomb_temp_cut, coulomb_den_cut
85 : real(dp), intent(inout) :: helm_res(num_helm_results)
86 : logical, intent(in) :: clip_to_table_boundaries, include_radiation, must_skip_elec_pos
87 : logical, intent(out) :: off_table
88 : integer, intent(out) :: ierr
89 :
90 0 : real(dp) :: h ! = planck_h
91 : type (Helm_Table), pointer :: ht
92 :
93 :
94 : !..declare local variables
95 : include 'helm_declare_local_variables.dek'
96 :
97 :
98 : !..given a temperature temp [K], density den [g/cm**3], and a composition
99 : !..characterized by abar and zbar, this routine returns most of the other
100 : !..thermodynamic quantities. of prime interest is the pressure [erg/cm**3],
101 : !..specific thermal energy [erg/gr], the entropy [erg/g/K], along with
102 : !..their derivatives with respect to temperature, density, abar, and zbar.
103 : !..other quantities such the normalized chemical potential eta (plus its
104 : !..derivatives), number density of electrons and positron pair (along
105 : !..with their derivatives), adiabatic indices, specific heats, and
106 : !..relativistically correct sound speed are also returned.
107 : !..
108 : !..this routine assumes planckian photons, an ideal gas of ions,
109 : !..and an electron-positron gas with an arbitrary degree of relativity
110 : !..and degeneracy. interpolation in a table of the helmholtz free energy
111 : !..is used to return the electron-positron thermodynamic quantities.
112 : !..all other derivatives are analytic.
113 : !..
114 : !..references: cox & giuli chapter 24 ; timmes & swesty apj 1999
115 :
116 : !..this routine assumes a call to subroutine read_helm_table has
117 : !..been performed prior to calling this routine.
118 :
119 :
120 : !..declare
121 :
122 0 : real(dp) :: abar, zbar, temp, logtemp, den, logden
123 : logical :: skip_elec_pos
124 :
125 : !..for the interpolations
126 : integer :: iat, jat
127 0 : real(dp) :: xt, xd, mxt, mxd, fi(36), &
128 0 : din, dindd, dinda, dindz, dindda, dinddz, dindaa, &
129 0 : dindaz, dindzz, dinddaa, dinddaz, &
130 : w0t, w1t, w2t, w0mt, w1mt, w2mt, &
131 : w0d, w1d, w2d, w0md, w1md, w2md, &
132 0 : dpepdd_in, dpepddd_in, dpepddt_in
133 :
134 : ! real(dp) psi0, dpsi0, ddpsi0, dddpsi0, &
135 : ! psi1, dpsi1, ddpsi1, dddpsi1, &
136 : ! psi2, dpsi2, ddpsi2, dddpsi2, &
137 : ! h5
138 :
139 : ! real(dp) xpsi0, xdpsi0, xddpsi0, &
140 : ! xpsi1, xdpsi1, xddpsi1, h3
141 :
142 0 : real(dp) :: si0t, si1t, si2t, si0mt, si1mt, si2mt, &
143 0 : si0d, si1d, si2d, si0md, si1md, si2md, &
144 0 : dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt, &
145 0 : dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md, &
146 0 : ddsi0t, ddsi1t, ddsi2t, ddsi0mt, ddsi1mt, ddsi2mt, &
147 0 : ddsi0d, ddsi1d, ddsi2d, ddsi0md, ddsi1md, ddsi2md, &
148 0 : dddsi0t, dddsi1t, dddsi2t, &
149 0 : dddsi0mt, dddsi1mt, dddsi2mt, &
150 0 : dddsi0d, dddsi1d, dddsi2d, &
151 0 : dddsi0md, dddsi1md, dddsi2md
152 :
153 0 : real(dp) :: free, df_d, df_t, df_dd, df_tt, df_dt, &
154 0 : df_ttt, df_dtt, df_ddt, df_ddd
155 :
156 0 : ht => eos_ht
157 :
158 0 : ierr = 0
159 0 : off_table = .false.
160 :
161 0 : h = planck_h
162 0 : third = 1.0d0/3.0d0
163 0 : sioncon = (2.0d0 * pi * amu * kerg)/(h*h)
164 0 : sifac = 8.6322745944370191d-45
165 0 : kergavo = kerg * avo
166 0 : asoli3 = asol/3.0d0
167 0 : clight2 = clight*clight
168 0 : eostol = 1.0d-13
169 0 : fpmin = 1.0d-14
170 : !..note: sifac = h**3/(2.0d0*pi*amu)**1.5d0
171 0 : forth = 4.0d0/3.0d0
172 0 : fiveth = 5.0d0/3.0d0
173 0 : teninth = 10.0d0/9.0d0
174 0 : esqu = qe*qe
175 0 : forthpi = forth * pi
176 :
177 0 : abar = abar_in
178 0 : zbar = zbar_in
179 0 : temp = temp_in
180 0 : logtemp = logtemp_in
181 0 : den = den_in
182 0 : logden = logden_in
183 :
184 : !..for very low T, convert all H to H2. adjust abar and zbar accordingly.
185 :
186 : ! NOTE: table lookup uses din rather than den
187 0 : ytot1 = 1.0d0/abar
188 0 : ye = ytot1 * zbar
189 0 : din = ye*den
190 :
191 0 : skip_elec_pos = must_skip_elec_pos
192 0 : if (.not. skip_elec_pos) then ! see if need to set it true
193 0 : if (temp < ht% templo) then
194 0 : ierr = 1
195 0 : return
196 : end if
197 :
198 0 : if (din < ht% denlo) then
199 0 : ierr = 1
200 0 : return
201 : end if
202 :
203 0 : if (temp > ht% temphi) then
204 0 : ierr = 1
205 0 : return
206 : end if
207 :
208 0 : if (din > ht% denhi) then
209 0 : ierr = 1
210 0 : return
211 : end if
212 : end if
213 :
214 :
215 0 : if (abar < 0d0) then
216 0 : ierr = 1
217 0 : return
218 : end if
219 :
220 : !..very neutron rich compositions may need to be bounded,
221 : !..avoid that extrema for now in order to increase efficiency
222 : ! ye = max(1.0d-16, ye)
223 :
224 :
225 : !..initialize local variables
226 : include 'helm_initialize_local_variables.dek'
227 : if (ierr /= 0) then
228 : if (dbg) write(*,*) 'failed in helm_initialize_local_variables'
229 : return
230 : end if
231 :
232 : !..radiation section:
233 0 : if (include_radiation) then
234 : include 'helm_radiation.dek'
235 : if (ierr /= 0) then
236 : if (dbg) write(*,*) 'failed in helm_radiation'
237 : return
238 : end if
239 : end if
240 :
241 : !..ion section:
242 : include 'helm_ideal_ions.dek'
243 : if (ierr /= 0) then
244 : if (dbg) write(*,*) 'failed in helm_ideal_ions'
245 : return
246 : end if
247 :
248 : !..electron-positron section:
249 0 : if (.not. skip_elec_pos) then
250 : include 'helm_electron_positron.dek'
251 : if (ierr /= 0) then
252 : if (dbg) write(*,*) 'failed in helm_electron_positron'
253 : return
254 : end if
255 : end if
256 :
257 : !..coulomb section:
258 0 : if ((ht% with_coulomb_corrections) .and. (.not. skip_elec_pos)) then
259 : include 'helm_coulomb2.dek'
260 : if (ierr /= 0) then
261 : if (dbg) write(*,*) 'failed in helm_coulomb2'
262 : return
263 : end if
264 : end if
265 :
266 : !..sum the gas and total (gas + radiation) components
267 : include 'helm_sum_totals.dek'
268 : if (ierr /= 0) then
269 : if (dbg) write(*,*) 'failed in helm_sum_totals'
270 : return
271 : end if
272 :
273 : ! error out for very low pgas,egas,sgas (previously just set hard floor at this value)
274 0 : if(pgas < 1d-20 .or. egas < 1d-20 .or. sgas < 1d-20) then
275 0 : ierr = 1
276 : if (dbg) write(*,*) 'failed in helm, sums too small'
277 0 : return
278 : end if
279 :
280 : !..compute the derivative quantities (cv, gamma1 ...etc)
281 : include 'helm_gammas.dek'
282 : if (ierr /= 0) then
283 : if (dbg) write(*,*) 'failed in helm_gammas'
284 : return
285 : end if
286 :
287 : !..maxwell relations; each is zero if the consistency is perfect
288 : !..if you don't need this, save three divides and comment this out
289 0 : x = den * den
290 0 : dse = temp*dentrdt/denerdt - 1.0d0
291 0 : dpe = (denerdd*x + temp*dpresdt)/pres - 1.0d0
292 0 : dsp = -dentrdd*x/dpresdt - 1.0d0
293 :
294 : !..store results
295 : include 'helm_store_results.dek'
296 0 : helm_res(h_crp:h_valid) = 0
297 :
298 : !..debugging printout
299 : if (.false.) then
300 : include 'helm_print_results.dek'
301 : end if
302 :
303 0 : end subroutine helmeos2aux
304 :
305 0 : subroutine show_h5(fi, ia, ja, w0t, w1t, w2t, w0mt, w1mt, w2mt, w0d, w1d, w2d, w0md, w1md, w2md)
306 : integer, intent(in) :: ia, ja
307 : real(dp), intent(in) :: fi(36), w0t, w1t, w2t, w0mt, w1mt, w2mt, w0d, w1d, w2d, w0md, w1md, w2md
308 0 : write(*,'(99a15)') 'w0t', 'w1t', 'w2t', 'w0mt', 'w1mt', 'w2mt', 'w0d', 'w1d', 'w2d', 'w0md', 'w1md', 'w2md'
309 0 : write(*,'(99e15.6)') w0t, w1t, w2t, w0mt, w1mt, w2mt, w0d, w1d, w2d, w0md, w1md, w2md
310 0 : write(*,'(a30,99e26.16)') 'fi(1)*w0d*w0t', fi(1)*w0d*w0t, fi(1),w0d,w0t
311 0 : write(*,'(a30,99e26.16)') 'fi(2)*w0md*w0t', fi(2)*w0md*w0t, fi(2),w0md,w0t
312 0 : write(*,'(a30,99e26.16)') 'fi(3)*w0d*w0mt', fi(3)*w0d*w0mt, fi(3),w0d,w0mt
313 0 : write(*,'(a30,99e26.16)') 'fi(4)*w0md*w0mt', fi(4)*w0md*w0mt, fi(4),w0md,w0mt
314 0 : write(*,'(a30,99e26.16)') '1 + 2 + 3 + 4', fi(1)*w0d*w0t + fi(2)*w0md*w0t + fi(3)*w0d*w0mt + fi(4)*w0md*w0mt
315 0 : write(*,'(A)')
316 0 : write(*,'(a30,99e26.16)') 'fi(5)*w0d*w1t', fi(5)*w0d*w1t, fi(5),w0d,w1t
317 0 : write(*,'(a30,99e26.16)') 'fi(6)*w0md*w1t', fi(6)*w0md*w1t, fi(6),w0md,w1t
318 0 : write(*,'(a30,99e26.16)') 'fi(7)*w0d*w1mt', fi(7)*w0d*w1mt, fi(7),w0d,w1mt
319 0 : write(*,'(a30,99e26.16)') 'fi(8)*w0md*w1mt', fi(8)*w0md*w1mt, fi(8),w0md,w1mt
320 0 : write(*,'(a30,99e26.16)') 'fi(9)*w0d*w2t', fi(9)*w0d*w2t, fi(9),w0d,w2t
321 0 : write(*,'(a30,99e26.16)') 'fi(10)*w0md*w2t', fi(10)*w0md*w2t, fi(10),w0md,w2t
322 0 : write(*,'(a30,99e26.16)') 'fi(11)*w0d*w2mt', fi(11)*w0d*w2mt, fi(11),w0d,w2mt
323 0 : write(*,'(a30,99e26.16)') 'fi(12)*w0md*w2mt', fi(12)*w0md*w2mt, fi(12),w0md,w2mt
324 0 : write(*,'(a30,99e26.16)') 'fi(13)*w1d*w0t', fi(13)*w1d*w0t, fi(13),w1d,w0t
325 0 : write(*,'(a30,99e26.16)') 'fi(14)*w1md*w0t', fi(14)*w1md*w0t, fi(14),w1md,w0t
326 0 : write(*,'(a30,99e26.16)') 'fi(15)*w1d*w0mt', fi(15)*w1d*w0mt, fi(15),w1d,w0mt
327 0 : write(*,'(a30,99e26.16)') 'fi(16)*w1md*w0mt', fi(16)*w1md*w0mt, fi(16),w1md,w0mt
328 0 : write(*,'(a30,99e26.16)') 'fi(17)*w2d*w0t', fi(17)*w2d*w0t, fi(17),w2d,w0t
329 0 : write(*,'(a30,99e26.16)') 'fi(18)*w2md*w0t', fi(18)*w2md*w0t, fi(18),w2md,w0t
330 0 : write(*,'(a30,99e26.16)') 'fi(19)*w2d*w0mt', fi(19)*w2d*w0mt, fi(19),w2d,w0mt
331 0 : write(*,'(a30,99e26.16)') 'fi(20)*w2md*w0mt', fi(20)*w2md*w0mt, fi(20),w2md,w0mt
332 0 : write(*,'(a30,99e26.16)') 'fi(21)*w1d*w1t', fi(21)*w1d*w1t, fi(21),w1d,w1t
333 0 : write(*,'(a30,99e26.16)') 'fi(22)*w1md*w1t', fi(22)*w1md*w1t, fi(22),w1md,w1t
334 0 : write(*,'(a30,99e26.16)') 'fi(23)*w1d*w1mt', fi(23)*w1d*w1mt, fi(23),w1d,w1mt
335 0 : write(*,'(a30,99e26.16)') 'fi(24)*w1md*w1mt', fi(24)*w1md*w1mt, fi(24),w1md,w1mt
336 0 : write(*,'(a30,99e26.16)') 'fi(25)*w2d*w1t', fi(25)*w2d*w1t, fi(25),w2d,w1t
337 0 : write(*,'(a30,99e26.16)') 'fi(26)*w2md*w1t', fi(26)*w2md*w1t, fi(26),w2md,w1t
338 0 : write(*,'(a30,99e26.16)') 'fi(27)*w2d*w1mt', fi(27)*w2d*w1mt, fi(27),w2d,w1mt
339 0 : write(*,'(a30,99e26.16)') 'fi(28)*w2md*w1mt', fi(28)*w2md*w1mt, fi(28),w2md,w1mt
340 0 : write(*,'(a30,99e26.16)') 'fi(29)*w1d*w2t', fi(29)*w1d*w2t, fi(29),w1d,w2t
341 0 : write(*,'(a30,99e26.16)') 'fi(30)*w1md*w2t', fi(30)*w1md*w2t, fi(30),w1md,w2t
342 0 : write(*,'(a30,99e26.16)') 'fi(31)*w1d*w2mt', fi(31)*w1d*w2mt, fi(31),w1d,w2mt
343 0 : write(*,'(a30,99e26.16)') 'fi(32)*w1md*w2mt', fi(32)*w1md*w2mt, fi(32),w1md,w2mt
344 0 : write(*,'(a30,99e26.16)') 'fi(33)*w2d*w2t', fi(33)*w2d*w2t, fi(33),w2d,w2t
345 0 : write(*,'(a30,99e26.16)') 'fi(34)*w2md*w2t', fi(34)*w2md*w2t, fi(34),w2md,w2t
346 0 : write(*,'(a30,99e26.16)') 'fi(35)*w2d*w2mt', fi(35)*w2d*w2mt, fi(35),w2d,w2mt
347 0 : write(*,'(a30,99e26.16)') 'fi(36)*w2md*w2mt', fi(36)*w2md*w2mt, fi(36),w2md,w2mt
348 0 : end subroutine show_h5
349 :
350 : end module helm
|
