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 eosPT_eval
21 : use eos_def
22 : use const_def, only: dp, crad, kerg, mp, ln10, arg_not_provided
23 : use math_lib
24 : use utils_lib, only: mesa_error
25 :
26 : implicit none
27 :
28 : integer, parameter :: doing_get_T = 1
29 : integer, parameter :: doing_get_Pgas = 2
30 :
31 : contains
32 :
33 24300 : subroutine Get_eosPT_Results(rq, &
34 : Z_in, X_in, abar, zbar, &
35 12150 : species, chem_id, net_iso, xa, &
36 : aPgas, alogPgas, atemp, alogtemp, &
37 : Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas, &
38 12150 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
39 : ierr)
40 : use utils_lib, only: is_bad
41 : type (EoS_General_Info), pointer :: rq
42 : real(dp), intent(in) :: Z_in, X_in, abar, zbar
43 : integer, intent(in) :: species
44 : integer, pointer :: chem_id(:), net_iso(:)
45 : real(dp), intent(in) :: xa(:)
46 : real(dp), intent(in) :: aPgas, alogPgas, atemp, alogtemp
47 : real(dp), intent(out) :: Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas
48 : real(dp), intent(inout) :: res(:), d_dlnRho_c_T(:), d_dlnT_c_Rho(:) ! (nv)
49 : real(dp), intent(inout) :: d_dxa_c_TRho(:,:) ! (nv, species)
50 : integer, intent(out) :: ierr
51 :
52 12150 : real(dp) :: X, Z, T, logT
53 12150 : real(dp) :: Pgas, logPgas, Prad, tiny
54 : logical, parameter :: dbg = .false.
55 :
56 : logical :: skip
57 :
58 : include 'formats'
59 :
60 12150 : ierr = 0
61 12150 : tiny = rq% tiny_fuzz
62 :
63 12150 : if (is_bad(X_in) .or. is_bad(Z_in)) then
64 0 : ierr = -1
65 0 : return
66 : end if
67 :
68 12150 : X = X_in; Z = Z_in
69 12150 : if (X < tiny) X = 0d0
70 12150 : if (Z < tiny) Z = 0d0
71 :
72 12150 : if (X > 1d0) then
73 0 : if (X > 1.0001D0) then
74 0 : write(*,1) 'Get_eosPT_Results: X bad', X
75 0 : ierr = -1
76 0 : return
77 : call mesa_error(__FILE__,__LINE__,'eosPT')
78 : end if
79 0 : X = 1d0
80 : end if
81 :
82 : call get_PT_args( &
83 12150 : aPgas, alogPgas, atemp, alogtemp, Pgas, logPgas, T, logT, ierr)
84 12150 : if (ierr /= 0) then
85 : if (dbg) write(*,*) 'error from get_PT_args'
86 : return
87 : end if
88 :
89 12150 : if (Pgas <= 0) then
90 0 : ierr = -1
91 0 : return
92 : end if
93 :
94 12150 : if (is_bad(Pgas) .or. is_bad(T)) then
95 0 : ierr = -1
96 0 : return
97 : end if
98 :
99 : call Get_PT_Results_using_DT( &
100 : rq, Z, X, abar, zbar, &
101 : species, chem_id, net_iso, xa, &
102 : Pgas, logPgas, T, logT, &
103 : Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas, &
104 12150 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, ierr)
105 :
106 12150 : if (ierr /= 0) then
107 : if (dbg) write(*,*) 'error from Get_PT_Results_using_DT'
108 : return
109 : end if
110 :
111 : end subroutine Get_eosPT_Results
112 :
113 :
114 12150 : subroutine get_PT_args( &
115 : aPg, alogPg, atemp, alogtemp, Pgas, logPgas, T, logT, ierr)
116 : real(dp), intent(in) :: aPg, alogPg
117 : real(dp), intent(in) :: atemp, alogtemp
118 : real(dp), intent(out) :: Pgas, logPgas, T, logT
119 : integer, intent(out) :: ierr
120 : include 'formats'
121 12150 : ierr = 0
122 12150 : T = atemp; logT = alogtemp
123 12150 : if (atemp == arg_not_provided .and. alogtemp == arg_not_provided) then
124 0 : ierr = -2; return
125 : end if
126 12150 : if (alogtemp == arg_not_provided) logT = log10(T)
127 12150 : if (atemp == arg_not_provided) T = exp10(logT)
128 12150 : if (T <= 0) then
129 0 : ierr = -1
130 0 : return
131 : end if
132 12150 : Pgas = aPg; logPgas = alogPg
133 12150 : if (Pgas == arg_not_provided .and. logPgas == arg_not_provided) then
134 0 : ierr = -3; return
135 : end if
136 12150 : if (logPgas == arg_not_provided) logPgas = log10(Pgas)
137 12150 : if (Pgas == arg_not_provided) Pgas = exp10(logPgas)
138 12150 : if (Pgas <= 0) then
139 0 : ierr = -1
140 0 : return
141 : end if
142 : end subroutine get_PT_args
143 :
144 :
145 12150 : subroutine Get_PT_Results_using_DT( &
146 : rq, Z, X, abar, zbar, &
147 12150 : species, chem_id, net_iso, xa, &
148 : Pgas, logPgas, T, logT, &
149 : Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas, &
150 12150 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, ierr)
151 : use eosDT_eval, only: get_Rho
152 : use utils_lib, only: is_bad
153 :
154 : type (EoS_General_Info), pointer :: rq ! general information about the request
155 : real(dp), intent(in) :: Z, X, abar, zbar
156 : integer, intent(in) :: species
157 : integer, pointer :: chem_id(:)
158 : integer, pointer :: net_iso(:)
159 : real(dp), intent(in) :: xa(:)
160 : real(dp), intent(inout) :: Pgas, logPgas, T, logT
161 : real(dp), intent(out) :: Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas
162 : real(dp), intent(inout) :: res(:) ! (nv)
163 : real(dp), intent(inout) :: d_dlnRho_c_T(:) ! (nv)
164 : real(dp), intent(inout) :: d_dlnT_c_Rho(:) ! (nv)
165 : real(dp), intent(inout) :: d_dxa_c_TRho(:,:) ! (nv, species)
166 : integer, intent(out) :: ierr
167 :
168 : integer:: i, eos_calls, max_iter, which_other
169 : real(dp) :: &
170 : logRho_guess, rho_guess, other, other_tol, logRho_tol, Prad, f, dfdx, &
171 12150 : logRho_bnd1, logRho_bnd2, other_at_bnd1, other_at_bnd2, logRho_result
172 :
173 : logical, parameter :: dbg = .false.
174 :
175 : include 'formats'
176 :
177 12150 : ierr = 0
178 :
179 12150 : which_other = i_lnPgas
180 12150 : other = logPgas*ln10
181 12150 : other_tol = 1d-8
182 12150 : logRho_tol = 1d-8
183 :
184 : ! guess based on fully ionized, ideal gas of ions and electrons
185 12150 : rho_guess = Pgas*abar*mp/(kerg*T*(1+zbar))
186 12150 : logRho_guess = log10(rho_guess)
187 :
188 12150 : logRho_bnd1 = arg_not_provided
189 12150 : logRho_bnd2 = arg_not_provided
190 12150 : other_at_bnd1 = arg_not_provided
191 12150 : other_at_bnd2 = arg_not_provided
192 :
193 12150 : max_iter = 20
194 12150 : eos_calls = 0
195 :
196 : if (dbg) write(*,1) 'rho_guess', rho_guess
197 : if (dbg) write(*,1) 'logRho_guess', logRho_guess
198 :
199 : call get_Rho( &
200 : rq% handle, Z, X, abar, zbar, &
201 : species, chem_id, net_iso, xa, &
202 : logT, which_other, other, &
203 : logRho_tol, other_tol, max_iter, logRho_guess, &
204 : logRho_bnd1, logRho_bnd2, other_at_bnd1, other_at_bnd2, &
205 : logRho_result, res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
206 12150 : eos_calls, ierr)
207 12150 : if (ierr /= 0) then
208 : if (dbg) then
209 : write(*,*) 'failed in get_Rho for Get_PT_Results_using_DT'
210 : write(*,1) 'Z = ', Z
211 : write(*,1) 'X = ', X
212 : write(*,1) 'abar = ', abar
213 : write(*,1) 'zbar = ', zbar
214 : write(*,1) 'logT = ', logT
215 : write(*,1) 'Pgas = ', Pgas
216 : write(*,1) 'logPgas = ', logPgas
217 : write(*,1) 'logRho_tol = ', logRho_tol
218 : write(*,1) 'other_tol = ', other_tol
219 : write(*,1) 'logRho_guess = ', logRho_guess
220 : write(*,'(A)')
221 : end if
222 : return
223 : end if
224 :
225 12150 : logRho = logRho_result
226 12150 : Rho = exp10(logRho)
227 :
228 : if (dbg) write(*,1) 'Rho', Rho
229 : if (dbg) write(*,1) 'logRho', logRho
230 : if (dbg) write(*,*)
231 : if (dbg) write(*,1) 'Pgas input', Pgas
232 : if (dbg) write(*,1) 'logPgas input', logPgas
233 : if (dbg) write(*,1) 'Pgas match', exp(res(i_lnPgas))
234 : if (dbg) write(*,1) 'logPgas match', res(i_lnPgas)/ln10
235 : if (dbg) write(*,*)
236 : if (dbg) write(*,1) 'get_Rho: grad_ad', res(i_grad_ad)
237 : if (dbg) write(*,*)
238 :
239 12150 : call do_partials
240 :
241 : contains
242 :
243 12150 : subroutine do_partials ! dlnRho_dlnPgas_c_T and dlnRho_dlnT_c_Pgas
244 12150 : real(dp) :: Prad, P, dP_dRho, dPgas_dRho, &
245 12150 : dP_dT, dPrad_dT, dPgas_dT, dRho_dPgas, dRho_dT
246 : include 'formats'
247 :
248 12150 : Prad = crad*T*T*T*T/3
249 12150 : P = Pgas + Prad
250 12150 : dP_dRho = res(i_chiRho)*P/Rho
251 12150 : dPgas_dRho = dP_dRho ! const T, so dP_dRho = dPgas_dRho
252 12150 : dRho_dPgas = 1/dPgas_dRho ! const T
253 12150 : dlnRho_dlnPgas_c_T = dRho_dPgas*Pgas/Rho ! const T
254 :
255 12150 : dPrad_dT = 4*crad*T*T*T/3
256 12150 : dP_dT = res(i_chiT)*P/T
257 12150 : dPgas_dT = dP_dT - dPrad_dT ! const Rho
258 12150 : dRho_dT = -dPgas_dT/dPgas_dRho ! const Pgas
259 12150 : dlnRho_dlnT_c_Pgas = dRho_dT*T/Rho
260 :
261 12150 : end subroutine do_partials
262 :
263 : end subroutine Get_PT_Results_using_DT
264 :
265 :
266 0 : subroutine get_T( &
267 : handle, Z, X, abar, zbar, &
268 0 : species, chem_id, net_iso, xa, &
269 : logPgas, which_other, other_value, &
270 : logT_tol, other_tol, max_iter, logT_guess, &
271 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
272 : logT_result, Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas, &
273 0 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
274 : eos_calls, ierr)
275 :
276 : integer, intent(in) :: handle
277 :
278 : real(dp), intent(in) :: Z ! the metals mass fraction
279 : real(dp), intent(in) :: X ! the hydrogen mass fraction
280 :
281 : real(dp), intent(in) :: abar, zbar
282 :
283 : integer, intent(in) :: species
284 : integer, pointer :: chem_id(:)
285 : integer, pointer :: net_iso(:)
286 : real(dp), intent(in) :: xa(:)
287 :
288 : real(dp), intent(in) :: logPgas ! log10 of density
289 : integer, intent(in) :: which_other
290 : real(dp), intent(in) :: other_value ! desired value for the other variable
291 : real(dp), intent(in) :: other_tol
292 :
293 : real(dp), intent(in) :: logT_tol
294 : integer, intent(in) :: max_iter ! max number of iterations
295 :
296 : real(dp), intent(in) :: logT_guess
297 : real(dp), intent(in) :: logT_bnd1, logT_bnd2 ! bounds for logT
298 : ! set to arg_not_provided if do not know bounds
299 : real(dp), intent(in) :: other_at_bnd1, other_at_bnd2 ! values at bounds
300 : ! if don't know these values, just set to arg_not_provided (defined in c_def)
301 :
302 : real(dp), intent(out) :: logT_result
303 : real(dp), intent(out) :: Rho, logRho ! density
304 : real(dp), intent(out) :: dlnRho_dlnPgas_c_T
305 : real(dp), intent(out) :: dlnRho_dlnT_c_Pgas
306 :
307 : real(dp), intent(inout) :: res(:) ! (nv)
308 : real(dp), intent(inout) :: d_dlnRho_c_T(:) ! (nv)
309 : real(dp), intent(inout) :: d_dlnT_c_Rho(:) ! (nv)
310 : real(dp), intent(inout) :: d_dxa_c_TRho(:,:) ! (nv, species)
311 :
312 : integer, intent(out) :: eos_calls
313 : integer, intent(out) :: ierr ! 0 means AOK.
314 :
315 : call do_safe_get_Pgas_T( &
316 : handle, Z, X, abar, zbar, &
317 : species, chem_id, net_iso, xa, &
318 : logPgas, which_other, other_value, doing_get_T, &
319 : logT_guess, logT_result, logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
320 : logT_tol, other_tol, max_iter, &
321 : Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas, &
322 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
323 0 : eos_calls, ierr)
324 :
325 0 : end subroutine get_T
326 :
327 :
328 0 : subroutine get_Pgas( &
329 : handle, Z, X, abar, zbar, &
330 0 : species, chem_id, net_iso, xa, &
331 : logT, which_other, other_value, &
332 : logPgas_tol, other_tol, max_iter, logPgas_guess, &
333 : logPgas_bnd1, logPgas_bnd2, other_at_bnd1, other_at_bnd2, &
334 : logPgas_result, Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas, &
335 0 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
336 : eos_calls, ierr)
337 :
338 : integer, intent(in) :: handle
339 :
340 : real(dp), intent(in) :: Z ! the metals mass fraction
341 : real(dp), intent(in) :: X ! the hydrogen mass fraction
342 :
343 : real(dp), intent(in) :: abar, zbar
344 :
345 : integer, intent(in) :: species
346 : integer, pointer :: chem_id(:)
347 : integer, pointer :: net_iso(:)
348 : real(dp), intent(in) :: xa(:)
349 :
350 : real(dp), intent(in) :: logT ! log10 of temperature
351 :
352 : integer, intent(in) :: which_other
353 : real(dp), intent(in) :: other_value ! desired value for the other variable
354 : real(dp), intent(in) :: other_tol
355 :
356 : real(dp), intent(in) :: logPgas_tol
357 :
358 : integer, intent(in) :: max_iter ! max number of Newton iterations
359 :
360 : real(dp), intent(in) :: logPgas_guess
361 : real(dp), intent(in) :: logPgas_bnd1, logPgas_bnd2 ! bounds for logPgas
362 : ! set to arg_not_provided if do not know bounds
363 : real(dp), intent(in) :: other_at_bnd1, other_at_bnd2 ! values at bounds
364 : ! if don't know these values, just set to arg_not_provided (defined in c_def)
365 :
366 : real(dp), intent(out) :: logPgas_result
367 : real(dp), intent(out) :: Rho, logRho ! density
368 : real(dp), intent(out) :: dlnRho_dlnPgas_c_T
369 : real(dp), intent(out) :: dlnRho_dlnT_c_Pgas
370 :
371 : real(dp), intent(inout) :: res(:) ! (nv)
372 : real(dp), intent(inout) :: d_dlnRho_c_T(:) ! (nv)
373 : real(dp), intent(inout) :: d_dlnT_c_Rho(:) ! (nv)
374 : real(dp), intent(inout) :: d_dxa_c_TRho(:,:) ! (nv, species)
375 :
376 : integer, intent(out) :: eos_calls
377 : integer, intent(out) :: ierr ! 0 means AOK.
378 :
379 : call do_safe_get_Pgas_T( &
380 : handle, Z, X, abar, zbar, &
381 : species, chem_id, net_iso, xa, &
382 : logT, which_other, other_value, doing_get_Pgas, &
383 : logPgas_guess, logPgas_result, logPgas_bnd1, logPgas_bnd2, other_at_bnd1, other_at_bnd2, &
384 : logPgas_tol, other_tol, max_iter, &
385 : Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas, &
386 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
387 0 : eos_calls, ierr)
388 :
389 0 : end subroutine get_Pgas
390 :
391 :
392 0 : subroutine do_safe_get_Pgas_T( &
393 : handle, Z, XH1, abar, zbar, &
394 0 : species, chem_id, net_iso, xa, &
395 : the_other_log, which_other, other_value, doing_which, &
396 : initial_guess, x, xbnd1, xbnd2, other_at_bnd1, other_at_bnd2, &
397 : xacc, yacc, ntry, &
398 : Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas, &
399 0 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
400 : eos_calls, ierr)
401 : use utils_lib, only: is_bad
402 : use num_lib, only: brent_safe_zero, look_for_brackets
403 : use chem_def, only: num_chem_isos
404 : integer, intent(in) :: handle
405 : real(dp), intent(in) :: Z, XH1, abar, zbar
406 : integer, intent(in) :: species
407 : integer, pointer :: chem_id(:)
408 : integer, pointer :: net_iso(:)
409 : real(dp), intent(in) :: xa(:)
410 : integer, intent(in) :: which_other
411 : real(dp), intent(in) :: other_value
412 : integer, intent(in) :: doing_which
413 : real(dp), intent(in) :: initial_guess ! for x
414 : real(dp), intent(out) :: x ! if doing_Pgas, then logPgas, else logT
415 : real(dp), intent(in) :: the_other_log
416 : real(dp), intent(in) :: xbnd1, xbnd2, other_at_bnd1, other_at_bnd2
417 : real(dp), intent(in) :: xacc, yacc ! tolerances
418 : integer, intent(in) :: ntry ! max number of iterations
419 : real(dp), intent(out) :: Rho, logRho ! density
420 : real(dp), intent(out) :: dlnRho_dlnPgas_c_T
421 : real(dp), intent(out) :: dlnRho_dlnT_c_Pgas
422 : real(dp), intent(inout) :: res(:) ! (nv)
423 : real(dp), intent(inout) :: d_dlnRho_c_T(:) ! (nv)
424 : real(dp), intent(inout) :: d_dlnT_c_Rho(:) ! (nv)
425 : real(dp), intent(inout) :: d_dxa_c_TRho(:,:) ! (nv, species)
426 : integer, intent(out) :: eos_calls, ierr
427 :
428 : integer :: i, j, lrpar, lipar, max_iter, irho, ix, iz
429 : real(dp), parameter :: dx = 0.1d0
430 0 : integer, pointer :: ipar(:)
431 0 : real(dp), pointer :: rpar(:)
432 0 : real(dp) :: Pgas, T, xb1, xb3, y1, y3, dfdx, f, logPgas, logT
433 : type (EoS_General_Info), pointer :: rq
434 :
435 : logical, parameter :: dbg = .false.
436 :
437 : include 'formats'
438 :
439 : ierr = 0
440 :
441 0 : call get_eos_ptr(handle, rq, ierr)
442 0 : if (ierr /= 0) then
443 0 : write(*, *) 'get_eos_ptr returned ierr', ierr
444 0 : return
445 : end if
446 :
447 0 : eos_calls = 0
448 0 : x = initial_guess
449 :
450 0 : if (doing_which /= doing_get_T) then
451 0 : Pgas = arg_not_provided
452 0 : T = exp10(the_other_log)
453 : else
454 0 : T = arg_not_provided
455 0 : Pgas = exp10(the_other_log)
456 : end if
457 :
458 0 : lipar = 0; nullify(ipar)
459 0 : lrpar = 0; nullify(rpar)
460 :
461 0 : xb1 = xbnd1; xb3 = xbnd2
462 0 : if (xb1 == arg_not_provided .or. xb3 == arg_not_provided .or. xb1 == xb3) then
463 :
464 : if (dbg) then
465 : write(*,'(A)')
466 : write(*,*) 'call look_for_brackets'
467 : write(*,2) 'ntry', ntry
468 : write(*,1) 'x', x
469 : write(*,1) 'dx', dx
470 : write(*,1) 'Z', Z
471 : write(*,1) 'XH1', XH1
472 : write(*,1) 'abar', abar
473 : write(*,1) 'zbar', zbar
474 : write(*,1) 'Pgas', Pgas
475 : write(*,1) 'T', T
476 : write(*,1) 'the_other_log', the_other_log
477 : write(*,'(A)')
478 : end if
479 :
480 : call look_for_brackets(x, dx, xb1, xb3, get_f_df, y1, y3, &
481 0 : ntry, lrpar, rpar, lipar, ipar, ierr)
482 0 : if (ierr /= 0) then
483 : if (dbg) then
484 : write(*, *) 'look_for_brackets returned ierr', ierr
485 : write(*,1) 'x', x
486 : write(*,1) 'dx', dx
487 : write(*,1) 'xb1', xb1
488 : write(*,1) 'xb3', xb3
489 : write(*,*) 'ntry', ntry
490 : write(*,*) 'lrpar', lrpar
491 : write(*,*) 'lipar', lipar
492 : end if
493 : return
494 : end if
495 : !write(*,*) 'done look_for_brackets'
496 : else
497 0 : if (other_at_bnd1 == arg_not_provided) then
498 0 : y1 = get_f_df(xb1, dfdx, lrpar, rpar, lipar, ipar, ierr)
499 0 : if (ierr /= 0) then
500 : return
501 : end if
502 : else
503 0 : y1 = other_at_bnd1 - other_value
504 : end if
505 0 : if (other_at_bnd2 == arg_not_provided) then
506 0 : y3 = get_f_df(xb3, dfdx, lrpar, rpar, lipar, ipar, ierr)
507 0 : if (ierr /= 0) then
508 : return
509 : end if
510 : else
511 0 : y3 = other_at_bnd2 - other_value
512 : end if
513 : end if
514 :
515 : if (dbg) then
516 : write(*,'(A)')
517 : write(*,*) 'call brent_safe_zero'
518 : write(*,1) 'xb1', xb1
519 : write(*,1) 'xb3', xb3
520 : write(*,1) 'y1', y1
521 : write(*,1) 'y3', y3
522 : end if
523 :
524 : x = brent_safe_zero( &
525 : xb1, xb3, 1d-14, 0.5d0*xacc, 0.5d0*yacc, get_f_df, y1, y3, &
526 0 : lrpar, rpar, lipar, ipar, ierr)
527 0 : if (ierr /= 0) then
528 : return
529 : end if
530 :
531 : contains
532 :
533 0 : real(dp) function get_f_df(x, dfdx, lrpar, rpar, lipar, ipar, ierr)
534 : integer, intent(in) :: lrpar, lipar
535 : real(dp), intent(in) :: x
536 : real(dp), intent(out) :: dfdx
537 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
538 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
539 : integer, intent(out) :: ierr
540 :
541 0 : real(dp) :: new, logT, logPgas
542 :
543 : include 'formats'
544 :
545 0 : ierr = 0
546 0 : get_f_df = 0
547 :
548 0 : dfdx = 0
549 :
550 0 : if (doing_which /= doing_get_T) then
551 0 : logPgas = x
552 0 : Pgas = exp10(logPgas)
553 0 : logT = the_other_log
554 0 : T = arg_not_provided
555 : else
556 0 : logT = x
557 0 : T = exp10(logT)
558 0 : logPgas = the_other_log
559 0 : Pgas = arg_not_provided
560 : end if
561 :
562 : ierr = 0
563 : call Get_eosPT_Results(rq, &
564 : Z, XH1, abar, zbar, &
565 : species, chem_id, net_iso, xa, &
566 : Pgas, logPgas, T, logT, &
567 : Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas, &
568 : res, d_dlnRho_c_T, d_dlnT_c_Rho, d_dxa_c_TRho, &
569 0 : ierr)
570 0 : if (ierr /= 0) then
571 : 22 format(a30, e26.16)
572 : if (dbg) then
573 : write(*, *) 'Get_eosPT_Results returned ierr', ierr
574 : write(*, 22) 'Z', Z
575 : write(*, 22) 'XH1', XH1
576 : write(*, 22) 'abar', abar
577 : write(*, 22) 'zbar', zbar
578 : write(*, 22) 'Pgas', Pgas
579 : write(*, 22) 'logPgas', logPgas
580 : write(*, 22) 'T', T
581 : write(*, 22) 'logT', logT
582 : write(*,'(A)')
583 : end if
584 : return
585 : end if
586 :
587 0 : eos_calls = eos_calls+1 ! count eos calls
588 :
589 0 : new = res(which_other)
590 0 : get_f_df = new - other_value
591 :
592 : ! f = f(lnRho(lnPgas,lnT),lnT)
593 0 : if (doing_which == doing_get_T) then
594 : dfdx = (d_dlnT_c_Rho(which_other) &
595 0 : + dlnRho_dlnT_c_Pgas*d_dlnRho_c_T(which_other))*ln10
596 0 : else if (doing_which == doing_get_Pgas) then
597 0 : dfdx = dlnRho_dlnPgas_c_T*d_dlnRho_c_T(which_other)*ln10
598 : else
599 0 : call mesa_error(__FILE__,__LINE__,'bad value for doing_which in eosPT_eval')
600 : end if
601 :
602 0 : end function get_f_df
603 :
604 : end subroutine do_safe_get_Pgas_T
605 :
606 : end module eosPT_eval
|
