Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-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 eos_HELM_eval
21 : use eos_def
22 : use const_def, only: avo, crad, ln10, arg_not_provided, mp, kerg, dp, qp
23 : use utils_lib, only: is_bad, mesa_error
24 : use math_lib
25 : use eospc_eval
26 : use helm
27 :
28 : implicit none
29 :
30 : logical, parameter :: stop_for_is_bad = .false.
31 : logical, parameter :: dbg = .false.
32 :
33 : contains
34 :
35 0 : subroutine get_helm_for_eosdt( &
36 : handle, dbg, Z, X, abar, zbar, &
37 0 : species, chem_id, net_iso, xa, &
38 : rho, logRho, T, logT, remaining_fraction, &
39 0 : res, d_dlnd, d_dlnT, d_dxa, &
40 : skip, ierr)
41 : use chem_lib, only: composition_info
42 : integer, intent(in) :: handle
43 : logical, intent(in) :: dbg
44 : real(dp), intent(in) :: &
45 : Z, X, abar, zbar, remaining_fraction
46 : integer, intent(in) :: species
47 : integer, pointer :: chem_id(:), net_iso(:)
48 : real(dp), intent(in) :: xa(:)
49 : real(dp), intent(in) :: rho, logRho, T, logT
50 : real(dp), intent(inout), dimension(nv) :: &
51 : res, d_dlnd, d_dlnT
52 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
53 : logical, intent(out) :: skip
54 : integer, intent(out) :: ierr
55 : type (EoS_General_Info), pointer :: rq
56 0 : real(dp), dimension(nv) :: d_dabar, d_dzbar
57 0 : real(dp) :: helm_res(num_helm_results)
58 :
59 0 : real(dp) :: xh, xhe, zz, abar_ci, zbar_ci, z2bar_ci, z53bar_ci, ye_ci, mass_correction, sumx
60 0 : real(dp), dimension(species) :: dabar_dx, dzbar_dx, dmc_dx
61 : integer :: i
62 :
63 0 : rq => eos_handles(handle)
64 : call Get_HELMEOS_Results( &
65 : rq, Z, abar, zbar, Rho, logRho, T, logT, &
66 : res, d_dlnd, d_dlnT, d_dabar, d_dzbar, &
67 0 : helm_res, skip, ierr)
68 0 : if (ierr /= 0) return
69 :
70 : ! need this call to get dabar_dx, dzbar_dx
71 : ! might want to pass these in eventually
72 : call composition_info( &
73 : species, chem_id, xa, xh, xhe, zz, &
74 : abar_ci, zbar_ci, z2bar_ci, z53bar_ci, ye_ci, mass_correction, &
75 0 : sumx, dabar_dx, dzbar_dx, dmc_dx)
76 :
77 : ! zero these for now
78 0 : d_dabar(i_phase:i_latent_ddlnRho) = 0d0
79 0 : d_dzbar(i_phase:i_latent_ddlnRho) = 0d0
80 0 : d_dabar(i_frac:i_frac+num_eos_frac_results-1) = 0d0
81 0 : d_dzbar(i_frac:i_frac+num_eos_frac_results-1) = 0d0
82 :
83 0 : do i=1, species
84 0 : d_dxa(:,i) = d_dabar(:)*dabar_dx(i) + d_dzbar(:)*dzbar_dx(i)
85 : end do
86 :
87 : ! zero phase information
88 0 : res(i_phase:i_latent_ddlnRho) = 0d0
89 0 : d_dlnT(i_phase:i_latent_ddlnRho) = 0d0
90 0 : d_dlnd(i_phase:i_latent_ddlnRho) = 0d0
91 :
92 : ! zero all components
93 0 : res(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
94 0 : d_dlnd(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
95 0 : d_dlnT(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
96 :
97 : ! mark this one
98 0 : res(i_frac_HELM) = 1.0d0
99 :
100 0 : end subroutine get_helm_for_eosdt
101 :
102 :
103 0 : subroutine Get_HELMEOS_Results( &
104 : rq, Z, abar, zbar, Rho, logRho, T, logT, &
105 : res, d_dlnd, d_dlnT, d_dabar, d_dzbar, &
106 : helm_res, off_table, ierr)
107 : type (EoS_General_Info), pointer :: rq
108 : real(dp), intent(in) :: Z, abar, zbar
109 : real(dp), intent(in) :: Rho, logRho, T, logT
110 : real(dp), intent(inout), dimension(nv) :: &
111 : res, d_dlnd, d_dlnT, d_dabar, d_dzbar
112 : real(dp), intent(inout) :: helm_res(num_helm_results)
113 : logical, intent(out) :: off_table
114 : integer, intent(out) :: ierr
115 :
116 : logical, parameter :: clip_to_table_boundaries = .true.
117 :
118 : logical :: include_elec_pos, include_radiation
119 :
120 : include 'formats'
121 :
122 : ierr = 0
123 : off_table = .false.
124 :
125 0 : include_elec_pos = rq% include_elec_pos
126 0 : include_radiation = rq% include_radiation
127 :
128 : call helmeos2( &
129 : T, logT, Rho, logRho, abar, zbar, &
130 : rq% coulomb_temp_cut_HELM, rq% coulomb_den_cut_HELM, &
131 : helm_res, clip_to_table_boundaries, include_radiation, &
132 : include_elec_pos, &
133 0 : off_table, ierr)
134 0 : if (off_table) return
135 0 : if (ierr /= 0) then
136 : if (dbg) then
137 : write(*,*) 'failed in helmeos2'
138 : write(*,1) 'T', T
139 : write(*,1) 'logT', logT
140 : write(*,1) 'Rho', Rho
141 : write(*,1) 'logRho', logRho
142 : write(*,1) 'abar', abar
143 : write(*,1) 'zbar', zbar
144 : write(*,*) 'clip_to_table_boundaries', clip_to_table_boundaries
145 : write(*,*) 'include_radiation', include_radiation
146 : write(*,*) 'include_elec_pos', include_elec_pos
147 : call mesa_error(__FILE__,__LINE__,'Get1_HELMEOS_Results')
148 : end if
149 : !write(*,*) 'failed in helmeos2'
150 : return
151 : end if
152 : call do_convert_helm_results( &
153 : helm_res, Z, abar, zbar, Rho, T, &
154 0 : res, d_dlnd, d_dlnT, d_dabar, d_dzbar, ierr)
155 0 : if (ierr /= 0) then
156 : if (dbg) write(*,*) 'failed in do_convert_helm_results'
157 : return
158 : end if
159 :
160 0 : end subroutine Get_HELMEOS_Results
161 :
162 :
163 0 : subroutine do_convert_helm_results( &
164 : helm_res, Z, abar, zbar, Rho, T, &
165 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dabar_c_TRho, d_dzbar_c_TRho, ierr)
166 : use helm
167 : use chem_def
168 : real(dp), intent(in) :: helm_res(num_helm_results)
169 : real(dp), intent(in) :: Z, abar, zbar, Rho, T
170 : real(dp), intent(inout) :: res(nv)
171 : real(dp), intent(inout) :: d_dlnRho_c_T(nv)
172 : real(dp), intent(inout) :: d_dlnT_c_Rho(nv)
173 : real(dp), intent(inout) :: d_dabar_c_TRho(nv)
174 : real(dp), intent(inout) :: d_dzbar_c_TRho(nv)
175 : integer, intent(out) :: ierr
176 :
177 0 : real(dp) :: mu, P, Pgas, energy, entropy, free_e, dse, dpe, dsp
178 : integer :: j, k, ci
179 :
180 : include 'formats'
181 :
182 0 : ierr = 0
183 :
184 : if (.false. .and. eos_test_partials) then
185 : eos_test_partials_val = helm_res(h_etot)
186 : eos_test_partials_dval_dx = helm_res(h_dea)
187 : write(*,1) 'logRho', log10(Rho)
188 : write(*,1) 'logT', log10(T)
189 : write(*,1) 'Rho', Rho
190 : write(*,1) 'T', T
191 : write(*,1) 'Z', Z
192 : write(*,1) 'abar', abar
193 : write(*,1) 'zbar', zbar
194 : write(*,1) 'etot', helm_res(h_etot)
195 : write(*,1) 'detot_dT', helm_res(h_det)
196 : write(*,1) 'detot_dRho', helm_res(h_ded)
197 : write(*,1) 'detot_dabar', helm_res(h_dea)
198 : write(*,1) 'detot_dzbar', helm_res(h_dez)
199 : call mesa_error(__FILE__,__LINE__,'do_convert_helm_results')
200 : end if
201 :
202 0 : energy = helm_res(h_etot)
203 0 : entropy = helm_res(h_stot)
204 0 : P = helm_res(h_ptot)
205 0 : Pgas = helm_res(h_pgas)
206 :
207 0 : res(i_lnE) = log(energy)
208 0 : res(i_lnS) = log(entropy)
209 0 : res(i_lnPgas) = log(Pgas)
210 :
211 0 : res(i_grad_ad) = helm_res(h_nabad)
212 0 : res(i_chiRho) = helm_res(h_chid)
213 0 : res(i_chiT) = helm_res(h_chit)
214 0 : res(i_Cp) = helm_res(h_cp)
215 0 : res(i_Cv) = helm_res(h_cv)
216 0 : res(i_dE_dRho) = helm_res(h_ded)
217 0 : res(i_dS_dT) = helm_res(h_dst)
218 0 : res(i_dS_dRho) = helm_res(h_dsd)
219 0 : mu = abar/(1 + zbar)
220 0 : res(i_mu) = mu
221 0 : free_e = max(1d-99, helm_res(h_xne))/(avo*Rho) ! assuming complete ionization
222 0 : res(i_lnfree_e) = log(free_e)
223 0 : res(i_gamma1) = helm_res(h_gam1)
224 0 : res(i_gamma3) = helm_res(h_gam3)
225 0 : res(i_eta) = helm_res(h_etaele)
226 :
227 0 : d_dlnRho_c_T(i_lnS) = helm_res(h_dsd)*Rho/entropy
228 0 : d_dlnRho_c_T(i_lnPgas) = helm_res(h_dpgasd)*Rho/Pgas
229 0 : d_dlnRho_c_T(i_lnE) = helm_res(h_ded)*Rho/energy
230 :
231 0 : d_dlnRho_c_T(i_grad_ad) = helm_res(h_dnabdd)*Rho
232 0 : d_dlnRho_c_T(i_chiRho) = helm_res(h_dchiddd)*Rho
233 0 : d_dlnRho_c_T(i_chiT) = helm_res(h_dchitdd)*Rho
234 0 : d_dlnRho_c_T(i_Cp) = helm_res(h_dcpdd)*Rho
235 0 : d_dlnRho_c_T(i_Cv) = helm_res(h_dcvdd)*Rho
236 0 : d_dlnRho_c_T(i_dE_dRho) = helm_res(h_dedd)*Rho
237 0 : d_dlnRho_c_T(i_dS_dT) = helm_res(h_dsdt)*Rho
238 0 : d_dlnRho_c_T(i_dS_dRho) = helm_res(h_dsdd)*Rho
239 0 : d_dlnRho_c_T(i_mu) = 0
240 0 : d_dlnRho_c_T(i_lnfree_e) = helm_res(h_dxned)/(avo*free_e) - 1d0
241 0 : d_dlnRho_c_T(i_gamma1) = helm_res(h_dgam1dd)*Rho
242 0 : d_dlnRho_c_T(i_gamma3) = helm_res(h_dgam3dd)*Rho
243 0 : d_dlnRho_c_T(i_eta) = helm_res(h_detad)*Rho
244 :
245 0 : d_dlnT_c_Rho(i_lnS) = helm_res(h_dst)*T/entropy
246 0 : d_dlnT_c_Rho(i_lnPgas) = helm_res(h_dpgast)*T/Pgas
247 0 : d_dlnT_c_Rho(i_lnE) = helm_res(h_det)*T/energy
248 :
249 0 : d_dlnT_c_Rho(i_grad_ad) = helm_res(h_dnabdt)*T
250 0 : d_dlnT_c_Rho(i_chiRho) = helm_res(h_dchiddt)*T
251 0 : d_dlnT_c_Rho(i_chiT) = helm_res(h_dchitdt)*T
252 0 : d_dlnT_c_Rho(i_Cp) = helm_res(h_dcpdt)*T
253 0 : d_dlnT_c_Rho(i_Cv) = helm_res(h_dcvdt)*T
254 0 : d_dlnT_c_Rho(i_dE_dRho) = helm_res(h_dedt)*T
255 0 : d_dlnT_c_Rho(i_dS_dT) = helm_res(h_dstt)*T
256 0 : d_dlnT_c_Rho(i_dS_dRho) = helm_res(h_dsdt)*T
257 0 : d_dlnT_c_Rho(i_mu) = 0
258 0 : d_dlnT_c_Rho(i_lnfree_e) = (helm_res(h_dxnet)*T/(avo*Rho))/free_e
259 0 : d_dlnT_c_Rho(i_gamma1) = helm_res(h_dgam1dt)*T
260 0 : d_dlnT_c_Rho(i_gamma3) = helm_res(h_dgam3dt)*T
261 0 : d_dlnT_c_Rho(i_eta) = helm_res(h_detat)*T
262 :
263 : d_dlnRho_c_T(i_lnE) = helm_res(h_ded)*Rho/energy
264 : d_dlnT_c_Rho(i_lnE) = helm_res(h_det)*T/energy
265 :
266 : d_dlnRho_c_T(i_lnS) = helm_res(h_dsd)*Rho/entropy
267 : d_dlnT_c_Rho(i_lnS) = helm_res(h_dst)*T/entropy
268 :
269 : d_dlnRho_c_T(i_lnPgas) = helm_res(h_dpgasd)*Rho/Pgas
270 : d_dlnT_c_Rho(i_lnPgas) = helm_res(h_dpgast)*T/Pgas
271 :
272 : ! composition partials need set
273 0 : d_dabar_c_TRho(i_lnS) = helm_res(h_dsa)/entropy
274 0 : d_dabar_c_TRho(i_lnPgas) = helm_res(h_dpgasa)/Pgas
275 0 : d_dabar_c_TRho(i_lnE) = helm_res(h_dea)/energy
276 :
277 0 : d_dabar_c_TRho(i_grad_ad) = helm_res(h_dnabda)
278 0 : d_dabar_c_TRho(i_chiRho) = helm_res(h_dchidda)
279 0 : d_dabar_c_TRho(i_chiT) = helm_res(h_dchitda)
280 0 : d_dabar_c_TRho(i_Cp) = helm_res(h_dcpda)
281 0 : d_dabar_c_TRho(i_Cv) = helm_res(h_dcvda)
282 0 : d_dabar_c_TRho(i_dE_dRho) = helm_res(h_deda)
283 0 : d_dabar_c_TRho(i_dS_dT) = helm_res(h_dsta)
284 0 : d_dabar_c_TRho(i_dS_dRho) = helm_res(h_dsda)
285 0 : d_dabar_c_TRho(i_mu) = 0
286 0 : d_dabar_c_TRho(i_lnfree_e) = (helm_res(h_dxnea)/(avo*Rho))/free_e
287 0 : d_dabar_c_TRho(i_gamma1) = helm_res(h_dgam1da)
288 0 : d_dabar_c_TRho(i_gamma3) = helm_res(h_dgam3da)
289 0 : d_dabar_c_TRho(i_eta) = helm_res(h_detaa)
290 :
291 0 : d_dzbar_c_TRho(i_lnS) = helm_res(h_dsz)/entropy
292 0 : d_dzbar_c_TRho(i_lnPgas) = helm_res(h_dpgasz)/Pgas
293 0 : d_dzbar_c_TRho(i_lnE) = helm_res(h_dez)/energy
294 :
295 0 : d_dzbar_c_TRho(i_grad_ad) = helm_res(h_dnabdz)
296 0 : d_dzbar_c_TRho(i_chiRho) = helm_res(h_dchiddz)
297 0 : d_dzbar_c_TRho(i_chiT) = helm_res(h_dchitdz)
298 0 : d_dzbar_c_TRho(i_Cp) = helm_res(h_dcpdz)
299 0 : d_dzbar_c_TRho(i_Cv) = helm_res(h_dcvdz)
300 0 : d_dzbar_c_TRho(i_dE_dRho) = helm_res(h_dedz)
301 0 : d_dzbar_c_TRho(i_dS_dT) = helm_res(h_dstz)
302 0 : d_dzbar_c_TRho(i_dS_dRho) = helm_res(h_dsdz)
303 0 : d_dzbar_c_TRho(i_mu) = 0
304 0 : d_dzbar_c_TRho(i_lnfree_e) = (helm_res(h_dxnez)/(avo*Rho))/free_e
305 0 : d_dzbar_c_TRho(i_gamma1) = helm_res(h_dgam1dz)
306 0 : d_dzbar_c_TRho(i_gamma3) = helm_res(h_dgam3dz)
307 0 : d_dzbar_c_TRho(i_eta) = helm_res(h_detaz)
308 :
309 :
310 0 : end subroutine do_convert_helm_results
311 :
312 :
313 0 : subroutine Get_HELM_Results( &
314 : abar, zbar, arho, alogrho, atemp, alogtemp, &
315 : coulomb_temp_cut, coulomb_den_cut, &
316 : include_radiation, include_elec_pos, &
317 0 : res, off_table, ierr)
318 0 : use const_def, only: dp
319 : use helm
320 :
321 : type (EoS_General_Info), pointer :: rq
322 : real(dp), intent(in) :: abar, zbar
323 : real(dp), intent(in) :: arho, alogrho
324 : real(dp), intent(in) :: atemp, alogtemp
325 : real(dp), intent(in) :: coulomb_temp_cut, coulomb_den_cut
326 : logical, intent(in) :: include_radiation, include_elec_pos
327 : real(dp), intent(inout) :: res(:) ! (num_helm_results)
328 : logical, intent(out) :: off_table
329 : integer, intent(out) :: ierr ! 0 means AOK.
330 :
331 0 : real(dp) :: Rho, logRho, T, logT, dse, dpe, dsp
332 :
333 : logical, parameter :: clip_to_table_boundaries = .true.
334 :
335 : include 'formats'
336 :
337 0 : ierr = 0
338 0 : off_table = .false.
339 :
340 : !..get temp and rho args
341 0 : T = atemp; logT = alogtemp
342 0 : if (atemp == arg_not_provided .and. alogtemp == arg_not_provided) then
343 0 : ierr = -2; return
344 : end if
345 0 : if (atemp == arg_not_provided) T = exp10(logT)
346 :
347 0 : Rho = arho; logrho = alogrho
348 0 : if (arho == arg_not_provided .and. alogrho == arg_not_provided) then
349 0 : ierr = -3; return
350 : end if
351 0 : if (arho == arg_not_provided) Rho = exp10(logRho)
352 :
353 : call helmeos2(T, logT, Rho, logRho, abar, zbar, &
354 : coulomb_temp_cut, coulomb_den_cut, &
355 : res, clip_to_table_boundaries, include_radiation, &
356 0 : include_elec_pos, off_table, ierr)
357 0 : res(h_valid) = 1
358 :
359 0 : end subroutine Get_HELM_Results
360 :
361 : end module eos_HELM_eval
|
