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_support
21 :
22 : use const_def, only: dp, ln10, arg_not_provided
23 : use star_private_def
24 : use utils_lib, only : is_bad, mesa_error
25 :
26 : implicit none
27 :
28 : private
29 : public :: get_eos
30 : public :: solve_eos_given_DE
31 : public :: solve_eos_given_DEgas
32 : public :: solve_eos_given_DP
33 : public :: solve_eos_given_DS
34 : public :: solve_eos_given_PT
35 : public :: solve_eos_given_PgasT
36 : public :: solve_eos_given_PgasT_auto
37 :
38 : integer, parameter :: MAX_ITER_FOR_SOLVE = 100
39 :
40 : contains
41 :
42 : ! Get eos results data given density & temperature
43 :
44 192472 : subroutine get_eos( &
45 192472 : s, k, xa, &
46 : Rho, logRho, T, logT, &
47 : res, dres_dlnRho, dres_dlnT, &
48 192472 : dres_dxa, ierr)
49 :
50 : use eos_lib, only: eosDT_get
51 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results, num_helm_results, i_lnE
52 :
53 : type (star_info), pointer :: s
54 : integer, intent(in) :: k ! 0 means not being called for a particular cell
55 : real(dp), intent(in) :: xa(:), Rho, logRho, T, logT
56 : real(dp), dimension(num_eos_basic_results), intent(out) :: &
57 : res, dres_dlnRho, dres_dlnT
58 : real(dp), intent(out) :: dres_dxa(num_eos_d_dxa_results,s% species)
59 : integer, intent(out) :: ierr
60 :
61 : integer :: j
62 :
63 : include 'formats'
64 :
65 192472 : ierr = 0
66 :
67 192472 : if (s% doing_timing) &
68 0 : s% timing_num_get_eos_calls = s% timing_num_get_eos_calls + 1
69 :
70 192472 : if(logRho < -25) then
71 : ! Provide some hard lower limit on what we would even try to evalue the eos at
72 : ! Going to low causes FPE's when we try to evaluate certain derivatives that need (rho**power)
73 0 : s% retry_message = 'eos evaluated at too low a density'
74 0 : ierr = -1
75 0 : return
76 : end if
77 :
78 : call eosDT_get( &
79 : s% eos_handle, s% species, s% chem_id, s% net_iso, xa, &
80 : Rho, logRho, T, logT, &
81 192472 : res, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
82 :
83 192472 : if (ierr /= 0) then
84 0 : s% retry_message = 'get_eos failed'
85 0 : if (s% report_ierr) then
86 0 : !$OMP critical (get_eos_critical)
87 0 : write(*,*) 'get_eos ierr', ierr
88 0 : write(*,2) 'k', k
89 0 : do j=1,s% species
90 0 : write(*,2) 'xa(j) ' // trim(s% nameofequ(j+s% nvar_hydro)), j, xa(j)
91 : end do
92 0 : write(*,1) 'log10Rho', logRho
93 0 : write(*,1) 'log10T', logT
94 0 : if (s% stop_for_bad_nums .and. &
95 0 : is_bad(logRho+logT)) call mesa_error(__FILE__,__LINE__,'do_eos_for_cell')
96 : !$OMP end critical (get_eos_critical)
97 : end if
98 0 : return
99 : end if
100 :
101 192472 : end subroutine get_eos
102 :
103 :
104 : ! Solve for temperature & eos results data given density & energy
105 :
106 1472 : subroutine solve_eos_given_DE( &
107 1472 : s, k, xa, &
108 : logRho, logE, logT_guess, logT_tol, logE_tol, &
109 1472 : logT, res, dres_dlnRho, dres_dlnT, dres_dxa, &
110 : ierr)
111 :
112 192472 : use eos_def
113 : use eos_lib, only: eosDT_get_T
114 :
115 : type (star_info), pointer :: s
116 : integer, intent(in) :: k ! 0 indicates not for a particular cell.
117 : real(dp), intent(in) :: &
118 : xa(:), logRho, logE, &
119 : logT_guess, logT_tol, logE_tol
120 : real(dp), intent(out) :: logT
121 : real(dp), dimension(num_eos_basic_results), intent(out) :: &
122 : res, dres_dlnRho, dres_dlnT
123 : real(dp), intent(out) :: dres_dxa(num_eos_d_dxa_results,s% species)
124 : integer, intent(out) :: ierr
125 :
126 : integer :: eos_calls
127 :
128 : include 'formats'
129 :
130 : ierr = 0
131 :
132 : call eosDT_get_T( &
133 : s% eos_handle, &
134 : s% species, s% chem_id, s% net_iso, xa, &
135 : logRho, i_lnE, logE*ln10, &
136 : logT_tol, logE_tol*ln10, MAX_ITER_FOR_SOLVE, logT_guess, &
137 : arg_not_provided, arg_not_provided, arg_not_provided, arg_not_provided, &
138 : logT, res, dres_dlnRho, dres_dlnT, dres_dxa, &
139 1472 : eos_calls, ierr)
140 :
141 1472 : if (s% doing_timing) s% timing_num_solve_eos_calls = s% timing_num_solve_eos_calls + eos_calls
142 :
143 1472 : end subroutine solve_eos_given_DE
144 :
145 :
146 : ! Solve for temperature & eos results data given density & gas energy
147 :
148 0 : subroutine solve_eos_given_DEgas( &
149 0 : s, k, xa, &
150 : logRho, egas, logT_guess, logT_tol, egas_tol, &
151 0 : logT, res, dres_dlnRho, dres_dlnT, dres_dxa, &
152 : ierr)
153 :
154 1472 : use eos_def
155 : use eos_lib, only: eosDT_get_T
156 :
157 : type (star_info), pointer :: s
158 : integer, intent(in) :: k ! 0 indicates not for a particular cell.
159 : real(dp), intent(in) :: &
160 : xa(:), logRho, egas, &
161 : logT_guess, logT_tol, egas_tol
162 : real(dp), intent(out) :: logT
163 : real(dp), dimension(num_eos_basic_results), intent(out) :: &
164 : res, dres_dlnRho, dres_dlnT
165 : real(dp), intent(out) :: dres_dxa(num_eos_d_dxa_results,s% species)
166 : integer, intent(out) :: ierr
167 :
168 : integer :: eos_calls
169 :
170 : include 'formats'
171 :
172 : ierr = 0
173 :
174 0 : if (s% doing_timing) s% timing_num_solve_eos_calls = s% timing_num_solve_eos_calls + 1
175 :
176 : call eosDT_get_T( &
177 : s% eos_handle, &
178 : s% species, s% chem_id, s% net_iso, xa, &
179 : logRho, i_egas, egas, logT_tol, egas_tol, MAX_ITER_FOR_SOLVE, logT_guess, &
180 : arg_not_provided, arg_not_provided, arg_not_provided, arg_not_provided, &
181 : logT, res, dres_dlnRho, dres_dlnT, &
182 0 : dres_dxa, eos_calls, ierr)
183 :
184 0 : end subroutine solve_eos_given_DEgas
185 :
186 :
187 : ! Solve for temperature & eos results data given density & pressure
188 :
189 0 : subroutine solve_eos_given_DP( &
190 0 : s, k, xa, &
191 : logRho, logP, logT_guess, logT_tol, logP_tol, &
192 0 : logT, res, dres_dlnRho, dres_dlnT, dres_dxa, &
193 : ierr)
194 :
195 0 : use eos_def
196 : use eos_lib, only: eosDT_get_T
197 :
198 : type (star_info), pointer :: s
199 : integer, intent(in) :: k ! 0 indicates not for a particular cell.
200 : real(dp), intent(in) :: &
201 : xa(:), logRho, logP, &
202 : logT_guess, logT_tol, logP_tol
203 : real(dp), intent(out) :: logT
204 : real(dp), dimension(num_eos_basic_results), intent(out) :: &
205 : res, dres_dlnRho, dres_dlnT
206 : real(dp), intent(out) :: dres_dxa(num_eos_d_dxa_results,s% species)
207 : integer, intent(out) :: ierr
208 :
209 : integer :: eos_calls
210 :
211 : include 'formats'
212 :
213 : ierr = 0
214 :
215 0 : if (s% doing_timing) s% timing_num_solve_eos_calls = s% timing_num_solve_eos_calls + 1
216 :
217 : call eosDT_get_T( &
218 : s% eos_handle, &
219 : s% species, s% chem_id, s% net_iso, xa, &
220 : logRho, i_logPtot, logP, logT_tol, logP_tol, MAX_ITER_FOR_SOLVE, logT_guess, &
221 : arg_not_provided, arg_not_provided, arg_not_provided, arg_not_provided, &
222 : logT, res, dres_dlnRho, dres_dlnT, &
223 0 : dres_dxa, eos_calls, ierr)
224 :
225 0 : end subroutine solve_eos_given_DP
226 :
227 :
228 : ! Solve for temperature & eos results data for a given density &
229 : ! entropy
230 :
231 0 : subroutine solve_eos_given_DS( &
232 0 : s, k, xa, &
233 : logRho, logS, logT_guess, logT_tol, logS_tol, &
234 0 : logT, res, dres_dlnRho, dres_dlnT, dres_dxa, &
235 : ierr)
236 :
237 0 : use eos_def
238 : use eos_lib, only: eosDT_get_T
239 :
240 : type (star_info), pointer :: s
241 : integer, intent(in) :: k ! 0 indicates not for a particular cell.
242 : real(dp), intent(in) :: &
243 : xa(:), logRho, logS, &
244 : logT_guess, logT_tol, logS_tol
245 : real(dp), intent(out) :: logT
246 : real(dp), dimension(num_eos_basic_results), intent(out) :: &
247 : res, dres_dlnRho, dres_dlnT
248 : real(dp), intent(out) :: dres_dxa(num_eos_d_dxa_results,s% species)
249 : integer, intent(out) :: ierr
250 :
251 : integer :: eos_calls
252 :
253 : include 'formats'
254 :
255 : ierr = 0
256 :
257 : call eosDT_get_T( &
258 : s% eos_handle, &
259 : s% species, s% chem_id, s% net_iso, xa, &
260 : logRho, i_lnS, logS*ln10, &
261 : logT_tol, logS_tol*ln10, MAX_ITER_FOR_SOLVE, logT_guess, &
262 : arg_not_provided, arg_not_provided, arg_not_provided, arg_not_provided, &
263 : logT, res, dres_dlnRho, dres_dlnT, dres_dxa, &
264 0 : eos_calls, ierr)
265 :
266 0 : if (s% doing_timing) s% timing_num_solve_eos_calls = s% timing_num_solve_eos_calls + eos_calls
267 :
268 0 : end subroutine solve_eos_given_DS
269 :
270 :
271 : ! Solve for density & eos results data given pressure & temperature
272 :
273 0 : subroutine solve_eos_given_PT( &
274 0 : s, k, xa, &
275 : logT, logP, logRho_guess, logRho_tol, logP_tol, &
276 0 : logRho, res, dres_dlnRho, dres_dlnT, dres_dxa, &
277 : ierr)
278 :
279 0 : use eos_def
280 : use eos_lib, only: eosDT_get_Rho
281 :
282 : type (star_info), pointer :: s
283 : integer, intent(in) :: k ! 0 indicates not for a particular cell.
284 : real(dp), intent(in) :: &
285 : xa(:), logT, logP, &
286 : logRho_guess, logRho_tol, logP_tol
287 : real(dp), intent(out) :: logRho
288 : real(dp), dimension(num_eos_basic_results), intent(out) :: &
289 : res, dres_dlnRho, dres_dlnT
290 : real(dp), intent(out) :: dres_dxa(num_eos_d_dxa_results,s% species)
291 : integer, intent(out) :: ierr
292 :
293 : integer :: eos_calls
294 :
295 : include 'formats'
296 :
297 : ierr = 0
298 :
299 0 : if (s% doing_timing) s% timing_num_solve_eos_calls = s% timing_num_solve_eos_calls + 1
300 :
301 : call eosDT_get_Rho( &
302 : s% eos_handle, &
303 : s% species, s% chem_id, s% net_iso, xa, &
304 : logT, i_logPtot, logP, logRho_tol, logP_tol, MAX_ITER_FOR_SOLVE, logRho_guess, &
305 : arg_not_provided, arg_not_provided, arg_not_provided, arg_not_provided, &
306 : logRho, res, dres_dlnRho, dres_dlnT, &
307 0 : dres_dxa, eos_calls, ierr)
308 :
309 0 : end subroutine solve_eos_given_PT
310 :
311 :
312 : ! Solve for density & eos results data given gas pressure &
313 : ! temperature
314 :
315 0 : subroutine solve_eos_given_PgasT( &
316 0 : s, k, xa, &
317 : logT, logPgas, logRho_guess, logRho_tol, logPgas_tol, &
318 0 : logRho, res, dres_dlnRho, dres_dlnT, dres_dxa, &
319 : ierr)
320 :
321 0 : use eos_def
322 : use eos_lib, only: eosDT_get_Rho
323 :
324 : type (star_info), pointer :: s
325 : integer, intent(in) :: k ! 0 indicates not for a particular cell.
326 : real(dp), intent(in) :: &
327 : xa(:), logT, logPgas, &
328 : logRho_guess, logRho_tol, logPgas_tol
329 : real(dp), intent(out) :: logRho
330 : real(dp), dimension(num_eos_basic_results), intent(out) :: &
331 : res, dres_dlnRho, dres_dlnT
332 : real(dp), intent(out) :: dres_dxa(num_eos_d_dxa_results,s% species)
333 : integer, intent(out) :: ierr
334 :
335 : integer :: eos_calls
336 :
337 : include 'formats'
338 :
339 : ierr = 0
340 :
341 : call eosDT_get_Rho( &
342 : s% eos_handle, &
343 : s% species, s% chem_id, s% net_iso, xa, &
344 : logT, i_lnPgas, logPgas*ln10, &
345 : logRho_tol, logPgas_tol*ln10, MAX_ITER_FOR_SOLVE, logRho_guess, &
346 : arg_not_provided, arg_not_provided, arg_not_provided, arg_not_provided, &
347 : logRho, res, dres_dlnRho, dres_dlnT, &
348 0 : dres_dxa, eos_calls, ierr)
349 0 : if (ierr /= 0 .and. s% report_ierr) then
350 0 : write(*,*) 'Call to eosDT_get_Rho failed in solve_eos_given_PgasT'
351 0 : write(*,2) 'logPgas', k, logPgas
352 0 : write(*,2) 'logT', k, logT
353 0 : write(*,2) 'logRho_guess', k, logRho_guess
354 : end if
355 :
356 0 : if (s% doing_timing) s% timing_num_solve_eos_calls = s% timing_num_solve_eos_calls + eos_calls
357 :
358 0 : end subroutine solve_eos_given_PgasT
359 :
360 :
361 : ! Solve for density & eos results data given gas pressure &
362 : ! temperature, with logRho_guess calculated automatically via an
363 : ! initial call to eos_gamma_PT_get
364 :
365 0 : subroutine solve_eos_given_PgasT_auto( &
366 0 : s, k, xa, &
367 : logT, logPgas, logRho_tol, logPgas_tol, &
368 0 : logRho, res, dres_dlnRho, dres_dlnT, dres_dxa, &
369 : ierr)
370 :
371 0 : use chem_lib, only: basic_composition_info
372 : use eos_def
373 : use eos_lib, only: eos_gamma_PT_get
374 :
375 : type (star_info), pointer :: s
376 : integer, intent(in) :: k ! 0 indicates not for a particular cell.
377 : real(dp), intent(in) :: &
378 : xa(:), logT, logPgas, &
379 : logRho_tol, logPgas_tol
380 : real(dp), intent(out) :: logRho
381 : real(dp), dimension(num_eos_basic_results), intent(out) :: &
382 : res, dres_dlnRho, dres_dlnT
383 : real(dp), intent(out) :: dres_dxa(num_eos_d_dxa_results,s% species)
384 : integer, intent(out) :: ierr
385 :
386 : real(dp) :: rho_guess, logRho_guess, gamma
387 :
388 : ! compute composition info
389 : real(dp) :: Y, Z, X, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
390 :
391 : call basic_composition_info( &
392 : s% species, s% chem_id, xa, X, Y, Z, &
393 0 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
394 :
395 0 : gamma = 5d0/3d0
396 : call eos_gamma_PT_get( &
397 : s% eos_handle, abar, exp10(logPgas), logPgas, exp10(logT), logT, gamma, &
398 : rho_guess, logRho_guess, res, dres_dlnRho, dres_dlnT, &
399 0 : ierr)
400 0 : if (ierr /= 0) then
401 : ierr = 0
402 0 : logRho_guess = arg_not_provided
403 : end if
404 :
405 : call solve_eos_given_PgasT( &
406 : s, k, xa, &
407 : logT, logPgas, logRho_guess, logRho_tol, logPgas_tol, &
408 : logRho, res, dres_dlnRho, dres_dlnT, dres_dxa, &
409 0 : ierr)
410 :
411 0 : end subroutine solve_eos_given_PgasT_auto
412 :
413 : end module eos_support
|
