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_thermodynamics
21 : use math_lib
22 : use auto_diff
23 :
24 : implicit none
25 :
26 : private
27 :
28 : public :: compute_derived_quantities, pack_for_export, thermodynamics_from_free_energy
29 :
30 : contains
31 :
32 : !> Computes the entropy, internal energy, and pressure of a system
33 : !! given its free energy, temperature, and density.
34 : !!
35 : !! @param F Free energy (erg/g)
36 : !! @param temp Temperature (K)
37 : !! @param den Density (g/cm^3)
38 : !! @param s Entropy (erg/g/K)
39 : !! @param e Internal energy (erg/g)
40 : !! @param p Pressure (erg/cm^3)
41 91040 : subroutine thermodynamics_from_free_energy(F, temp, den, s, e, p)
42 : type(auto_diff_real_2var_order3), intent(in) :: F
43 : type(auto_diff_real_2var_order3), intent(in) :: temp, den
44 : type(auto_diff_real_2var_order3), intent(out) :: s, e, p
45 :
46 : ! Entropy
47 : ! s = - F/dT
48 45520 : s = -differentiate_1(F)
49 :
50 : ! Pgas
51 : ! p = -dF/dV = -(1/V)dF/dlnV = (1/V)dF/dlnRho = Rho dF/dlnRho = Rho^2 dF/dRho
52 45520 : p = pow2(den) * differentiate_2(F)
53 :
54 : ! Energy
55 : ! E = F + T * S
56 45520 : e = F + temp * s
57 :
58 45520 : end subroutine thermodynamics_from_free_energy
59 :
60 : !> Computes various thermodynamic derivatives and related quantities given
61 : !! the entropy, internal energy, pressure, temperature, and density of a system.
62 : !!
63 : !! @param temp Temperature (K)
64 : !! @param den Density (g/cm^3)
65 : !! @param s Entropy (erg/g/K)
66 : !! @param e Internal energy (erg/g)
67 : !! @param p Pressure (erg/cm^3)
68 : !! @param cv Specific heat at constant volume (erg/g/K)
69 : !! @param cp Specific heat at constant pressure (erg/g/K)
70 : !! @param chit Thermal susceptibility (dlnP/dlnT at constant Rho)
71 : !! @param chid Volume susceptibility (dlnP/dlnRho at constant T)
72 : !! @param gam1 First adiabatic index
73 : !! @param gam2 Second adiabatic index
74 : !! @param gam3 Third adiabatic index
75 : !! @param nabad Adiabatic gradient (dlnT/dlnP at constant entropy)
76 : !! @param cs Sound speed (cm/s)
77 45520 : subroutine compute_derived_quantities(temp, dens, s, e, p, cv, cp, chit, chid, gam1, gam2, gam3, nabad, cs)
78 : use const_def, only: clight
79 : type(auto_diff_real_2var_order3), intent(in) :: temp, dens, s, e, p
80 : type(auto_diff_real_2var_order3), intent(out) :: cv, cp, chit, chid, gam1, gam2, gam3, nabad, cs
81 :
82 : ! Susceptibilities
83 45520 : chit = differentiate_1(p) * temp / p
84 45520 : chid = differentiate_2(p) * dens / p
85 :
86 : ! Specific heat at constant volume
87 45520 : cv = differentiate_1(e)
88 :
89 : ! Adiabatic indices
90 45520 : gam3 = 1d0 + (p / dens) * chit / (temp * cv)
91 45520 : gam1 = chid + (gam3 - 1d0) * chit
92 45520 : nabad = (gam3 - 1d0) / gam1
93 45520 : gam2 = 1d0 - nabad
94 :
95 : ! Specific heat at constant pressure
96 45520 : cp = cv * gam1 / chid
97 :
98 : ! Sound speed
99 45520 : cs = clight * sqrt(gam1 / (1d0 + (dens / p) * (e + clight*clight)))
100 45520 : end subroutine compute_derived_quantities
101 :
102 : !> Computes thermodynamic quantities from Skye and packs them into the EOS return vectors.
103 : !!
104 : !! @param F_ideal_ion Ideal ion free energy (erg/g)
105 : !! @param F_coul Coulomb corrections to the free energy (erg/g)
106 : !! @param F_rad Radiation free energy (erg/g)
107 : !! @param F_ele Ideal electron-positron free energy (erg/g)
108 : !! @param temp Temperature (K)
109 : !! @param den Density (g/cm^3)
110 : !! @param xnefer Electron density (1/cm^3)
111 : !! @param etaele The electron chemical potential in units of kT.
112 : !! @param abar The mean atomic mass number.
113 : !! @param zbar The mean atomic charge number.
114 : !! @param abar The mean atomic mass number.
115 : !! @param phase The blended phase. 0 for liquid, 1 for solid, smoothly interpolates in between.
116 : !! @param latent_ddlnT The latent heat of the smoothed phase transition in lnT (T dS/dlnT)
117 : !! @param latent_ddlnRho The latent heat of the smoothed phase transition in lnRho (T dS/dlnRho)
118 : !! @param include_radiation True to include radiation effects, false otherwise.
119 : !! @param res The EOS return vector.
120 : !! @param d_dlnRho The derivative of the EOS return vector with respect to lnRho.
121 : !! @param d_dlnT The derivative of the EOS return vector with respect to lnT.
122 45520 : subroutine pack_for_export(F_ideal_ion, F_coul, F_rad, F_ele, temp, dens, xnefer, etaele, abar, zbar, &
123 : phase, latent_ddlnT, latent_ddlnRho, res, d_dlnRho, d_dlnT, ierr)
124 : use const_def, only: dp, crad, avo
125 : use eos_def
126 : type(auto_diff_real_2var_order3), intent(in) :: F_ideal_ion, F_coul, F_rad, F_ele, temp, dens, xnefer, etaele
127 : type(auto_diff_real_2var_order3), intent(in) :: phase, latent_ddlnT, latent_ddlnRho
128 : real(dp), intent(in) :: abar, zbar
129 :
130 : integer, intent(out) :: ierr
131 :
132 : ! Intermediates
133 : type(auto_diff_real_2var_order3) :: srad, erad, prad, sgas, egas, pgas, p, e, s
134 : type(auto_diff_real_2var_order3) :: F_gas, lnS, lnE, lnPgas, mu, lnfree_e
135 : type(auto_diff_real_2var_order3) :: cv, cp, chit, chid, gam1, gam2, gam3, nabad, cs
136 :
137 : ! Outputs
138 : real(dp), intent(inout) :: res(nv)
139 : real(dp), intent(inout) :: d_dlnRho(nv)
140 : real(dp), intent(inout) :: d_dlnT(nv)
141 :
142 45520 : ierr = 0
143 :
144 : ! Form the electron-positron-ion gas plus Coulomb corrections
145 45520 : F_gas = F_ideal_ion + F_coul + F_ele
146 :
147 : ! Compute base thermodynamic quantities
148 45520 : call thermodynamics_from_free_energy(F_gas, temp, dens, sgas, egas, pgas)
149 :
150 : ! Write the radiation terms explicitly to avoid having rho^2/rho^2 term that
151 : ! auto_diff gives for prad when calling thermodynamics_from_free_energy on F_rad.
152 : ! This avoids some subtractions for quantities that should be 0 like chid.
153 45520 : prad = crad * pow4(temp) / 3d0
154 45520 : erad = crad * pow4(temp) / dens
155 45520 : srad = 4d0 * crad * pow3(temp) / (3d0 * dens)
156 :
157 45520 : p = prad + pgas
158 45520 : e = erad + egas
159 45520 : s = srad + sgas
160 :
161 45520 : if(s<0 .or. e<0 .or. pgas<0) then
162 0 : ierr = -1
163 : return
164 : end if
165 :
166 45520 : lnS = log(s)
167 45520 : lnE = log(e)
168 45520 : lnPgas = log(pgas)
169 :
170 : ! assuming complete ionization
171 45520 : mu = abar / (1d0 + zbar)
172 45520 : lnfree_e = log(max(1d-99, xnefer)/(avo*dens))
173 :
174 :
175 45520 : call compute_derived_quantities(temp, dens, s, e, p, cv, cp, chit, chid, gam1, gam2, gam3, nabad, cs)
176 :
177 45520 : res(i_lnS) = lnS%val
178 45520 : res(i_lnE) = lnE%val
179 45520 : res(i_lnPgas) = lnPgas%val
180 45520 : res(i_mu) = mu%val
181 45520 : res(i_grad_ad) = nabad%val
182 45520 : res(i_chiRho) = chid%val
183 45520 : res(i_chiT) = chit%val
184 45520 : res(i_Cp) = cp%val
185 45520 : res(i_Cv) = cv%val
186 45520 : res(i_dE_dRho) = e%d1val2
187 45520 : res(i_dS_dT) = s%d1val1
188 45520 : res(i_dS_dRho) = s%d1val2
189 45520 : res(i_lnfree_e) = lnfree_e%val
190 45520 : res(i_gamma1) = gam1%val
191 45520 : res(i_gamma3) = gam3%val
192 45520 : res(i_eta) = etaele%val
193 45520 : res(i_phase) = phase%val
194 45520 : res(i_latent_ddlnT) = latent_ddlnT%val
195 45520 : res(i_latent_ddlnRho) = latent_ddlnRho%val
196 :
197 45520 : d_dlnT(i_lnS) = lnS%d1val1 * temp%val
198 45520 : d_dlnT(i_lnE) = lnE%d1val1 * temp%val
199 45520 : d_dlnT(i_lnPgas) = lnPgas%d1val1 * temp%val
200 45520 : d_dlnT(i_mu) = mu%d1val1 * temp%val
201 45520 : d_dlnT(i_grad_ad) = nabad%d1val1 * temp%val
202 45520 : d_dlnT(i_chiRho) = chid%d1val1 * temp%val
203 45520 : d_dlnT(i_chiT) = chit%d1val1 * temp%val
204 45520 : d_dlnT(i_Cp) = cp%d1val1 * temp%val
205 45520 : d_dlnT(i_Cv) = cv%d1val1 * temp%val
206 45520 : d_dlnT(i_dE_dRho) = e%d1val1_d1val2 * temp%val
207 45520 : d_dlnT(i_dS_dT) = s%d2val1 * temp%val
208 45520 : d_dlnT(i_dS_dRho) = s%d1val1_d1val2 * temp%val
209 45520 : d_dlnT(i_lnfree_e) = lnfree_e%d1val1 * temp%val
210 45520 : d_dlnT(i_gamma1) = gam1%d1val1 * temp%val
211 45520 : d_dlnT(i_gamma3) = gam3%d1val1 * temp%val
212 45520 : d_dlnT(i_eta) = etaele%d1val1 * temp%val
213 45520 : d_dlnT(i_phase) = phase%d1val1 * temp%val
214 45520 : d_dlnT(i_latent_ddlnT) = latent_ddlnT%d1val1 * temp%val
215 45520 : d_dlnT(i_latent_ddlnRho) = latent_ddlnRho%d1val1 * temp%val
216 :
217 45520 : d_dlnRho(i_lnS) = lnS%d1val2 * dens%val
218 45520 : d_dlnRho(i_lnE) = lnE%d1val2 * dens%val
219 45520 : d_dlnRho(i_lnPgas) = lnPgas%d1val2 * dens%val
220 45520 : d_dlnRho(i_mu) = mu%d1val2 * dens%val
221 45520 : d_dlnRho(i_grad_ad) = nabad%d1val2 * dens%val
222 45520 : d_dlnRho(i_chiRho) = chid%d1val2 * dens%val
223 45520 : d_dlnRho(i_chiT) = chit%d1val2 * dens%val
224 45520 : d_dlnRho(i_Cp) = cp%d1val2 * dens%val
225 45520 : d_dlnRho(i_Cv) = cv%d1val2 * dens%val
226 45520 : d_dlnRho(i_dE_dRho) = e%d2val2 * dens%val
227 45520 : d_dlnRho(i_dS_dT) = s%d1val1_d1val2 * dens%val
228 45520 : d_dlnRho(i_dS_dRho) = s%d2val2 * dens%val
229 45520 : d_dlnRho(i_lnfree_e) = lnfree_e%d1val2 * dens%val
230 45520 : d_dlnRho(i_gamma1) = gam1%d1val2 * dens%val
231 45520 : d_dlnRho(i_gamma3) = gam3%d1val2 * dens%val
232 45520 : d_dlnRho(i_eta) = etaele%d1val2 * dens%val
233 45520 : d_dlnRho(i_phase) = phase%d1val2 * dens%val
234 45520 : d_dlnRho(i_latent_ddlnT) = latent_ddlnT%d1val2 * dens%val
235 45520 : d_dlnRho(i_latent_ddlnRho) = latent_ddlnRho%d1val2 * dens%val
236 :
237 45520 : end subroutine pack_for_export
238 :
239 :
240 :
241 : end module skye_thermodynamics
|
