Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2011-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 condint
21 :
22 : use const_def, only: dp, iln10
23 : use math_lib
24 : use utils_lib, only: mesa_error
25 :
26 : implicit none
27 :
28 : private
29 : public :: init_potekhin
30 : public :: do_electron_conduction
31 :
32 : integer, parameter :: num_logTs=29, num_logRhos=71, num_logzs=15
33 : !!! NB: These parameters must be consistent with the table "condtabl.d"!
34 : logical :: initialized = .false.
35 :
36 : real(dp) :: logTs(num_logTs), logRhos(num_logRhos), logzs(num_logzs)
37 : real(dp), target :: f_ary(4*num_logRhos*num_logTs*num_logzs) ! for bicubic splines
38 : real(dp), pointer :: f(:,:,:,:)
39 : integer :: ilinx(num_logzs), iliny(num_logzs)
40 :
41 : contains
42 :
43 11 : subroutine init_potekhin(ierr)
44 : use kap_def, only: kap_dir
45 : use interp_2d_lib_db, only: interp_mkbicub_db
46 : integer, intent(out) :: ierr
47 :
48 : character (len=256) :: filename
49 : integer :: read_err, iz, it, ir, shift
50 : integer :: ibcxmin ! bc flag for x=xmin
51 330 : real(dp) :: bcxmin(num_logTs) ! bc data vs. y at x=xmin
52 : integer :: ibcxmax ! bc flag for x=xmax
53 330 : real(dp) :: bcxmax(num_logTs) ! bc data vs. y at x=xmax
54 : integer :: ibcymin ! bc flag for y=ymin
55 792 : real(dp) :: bcymin(num_logRhos) ! bc data vs. x at y=ymin
56 : integer :: ibcymax ! bc flag for y=ymax
57 792 : real(dp) :: bcymax(num_logRhos) ! bc data vs. x at y=ymax
58 11 : real(dp) :: Z
59 11 : real(dp), pointer :: f1(:)
60 :
61 : include 'formats'
62 :
63 11 : ierr = 0
64 11 : if (initialized) return
65 :
66 5 : shift = 4*num_logRhos*num_logTs
67 5 : f(1:4,1:num_logRhos,1:num_logTs,1:num_logzs) => f_ary(1:shift*num_logzs)
68 :
69 5 : filename = trim(kap_dir) // '/condtabl.data'
70 5 : open(1,file=trim(filename),status='OLD',iostat=ierr)
71 5 : if (ierr /= 0) then
72 0 : write(*,'(A)')
73 0 : write(*,'(A)')
74 0 : write(*,'(A)')
75 0 : write(*,'(A)')
76 0 : write(*,'(A)')
77 0 : write(*,*) 'NOTICE: missing ' // trim(filename)
78 0 : write(*,*) 'Please remove the directory mesa/data/kap_data,'
79 0 : write(*,*) 'and rerun the mesa ./install script.'
80 0 : write(*,'(A)')
81 0 : call mesa_error(__FILE__,__LINE__)
82 : end if
83 : !print*,'Reading thermal conductivity data...'
84 5 : read_err = 0
85 5 : read(1,'(A)',iostat=read_err) ! skip the first line
86 5 : if (read_err /= 0) ierr = read_err
87 80 : do iz = 1, num_logzs
88 75 : read (1,*,iostat=read_err) z, (logTs(it),it=1,num_logTs)
89 75 : if (read_err /= 0) ierr = read_err
90 75 : if (z == 1d0) then
91 5 : logzs(iz) = 0d0
92 : else
93 70 : logzs(iz) = log10(z)
94 : end if
95 5405 : do ir = 1, num_logRhos
96 5325 : read(1,*,iostat=read_err) logRhos(ir), (f(1,ir,it,iz),it=1,num_logTs)
97 5400 : if (read_err /= 0) ierr = read_err
98 : end do
99 : end do
100 5 : close(1)
101 5 : if (ierr /= 0) then
102 0 : write(*,'(A)')
103 0 : write(*,'(A)')
104 0 : write(*,'(A)')
105 0 : write(*,'(A)')
106 0 : write(*,'(A)')
107 0 : write(*,*) 'NOTICE: error trying to read ' // trim(filename)
108 0 : write(*,*) 'Please remove the directory mesa/data/kap_data,'
109 0 : write(*,*) 'and rerun the mesa ./install script.'
110 0 : write(*,'(A)')
111 0 : call mesa_error(__FILE__,__LINE__)
112 : end if
113 : ! boundary condition is slope=0 at edges of tables
114 : ! this ensures continuous derivatives when we clip
115 5 : ibcxmin = 3; bcxmin(1:num_logTs) = 0d0
116 5 : ibcxmax = 3; bcxmax(1:num_logTs) = 0d0
117 5 : ibcymin = 3; bcymin(1:num_logRhos) = 0d0
118 5 : ibcymax = 3; bcymax(1:num_logRhos) = 0d0
119 80 : do iz = 1, num_logzs
120 75 : f1(1:shift) => f_ary(1+(iz-1)*shift:iz*shift)
121 : call interp_mkbicub_db( &
122 : logRhos, num_logRhos, logTs, num_logTs, f1, num_logRhos, &
123 : ibcxmin, bcxmin, ibcxmax, bcxmax, &
124 : ibcymin, bcymin, ibcymax, bcymax, &
125 75 : ilinx(iz), iliny(iz), ierr)
126 80 : if (ierr /= 0) then
127 0 : write(*,'(A)')
128 0 : write(*,'(A)')
129 0 : write(*,'(A)')
130 0 : write(*,'(A)')
131 0 : write(*,'(A)')
132 0 : write(*,*) 'NOTICE: error in ' // trim(filename)
133 0 : write(*,*) 'Please report the problem.'
134 0 : write(*,'(A)')
135 0 : call mesa_error(__FILE__,__LINE__)
136 : end if
137 : end do
138 5 : initialized = .true.
139 11 : end subroutine init_potekhin
140 :
141 :
142 86216 : subroutine do_electron_conduction_potekhin( &
143 : zbar, logRho_in, logT_in, kap, dlogkap_dlogRho, dlogkap_dlogT, ierr)
144 :
145 : use const_def, only: boltz_sigma
146 : real(dp), intent(in) :: zbar, logRho_in, logT_in
147 : real(dp), intent(out) :: kap, dlogkap_dlogRho, dlogkap_dlogT
148 : integer, intent(out) :: ierr
149 :
150 : integer :: iz, iz1, iz2, shift
151 86216 : real(dp) :: zlog, logRho, logT
152 86216 : real(dp) :: alfa, beta, &
153 86216 : logK1, kap1, dlogK1_dlogRho, dlogK1_dlogT, &
154 86216 : logK2, kap2, dlogK2_dlogRho, dlogK2_dlogT, &
155 86216 : logK, logkap, dlogK_dlogRho, dlogK_dlogT
156 86216 : real(dp), pointer :: f1(:)
157 :
158 : logical :: clipped_logRho, clipped_logT
159 :
160 : include 'formats'
161 :
162 86216 : ierr = 0
163 86216 : shift = 4*num_logRhos*num_logTs
164 :
165 86216 : if (logRho_in < logRhos(1)) then
166 22832 : logRho = logRhos(1)
167 22832 : clipped_logRho = .true.
168 63384 : else if (logRho_in > logRhos(num_logRhos)) then
169 0 : logRho = logRhos(num_logRhos)
170 0 : clipped_logRho = .true.
171 : else
172 63384 : logRho = logRho_in
173 63384 : clipped_logRho = .false.
174 : end if
175 :
176 86216 : if (logT_in < logTs(1)) then
177 0 : logT = logTs(1)
178 0 : clipped_logT = .true.
179 86216 : else if (logT_in > logTs(num_logTs)) then
180 0 : logT = logTs(num_logTs)
181 0 : clipped_logT = .true.
182 : else
183 86216 : logT = logT_in
184 86216 : clipped_logT = .false.
185 : end if
186 :
187 86216 : zlog = max(logzs(1),min(logzs(num_logzs),log10(max(1d-30,zbar))))
188 :
189 86216 : if (zlog <= logzs(1)) then ! use 1st
190 0 : call get1(1, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
191 86216 : else if (zlog >= logzs(num_logzs)) then ! use last
192 0 : call get1(num_logzs, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
193 : else ! interpolate
194 86216 : iz1 = -1
195 86218 : do iz = 2, num_logzs
196 86218 : if (zlog >= logzs(iz-1) .and. zlog <= logzs(iz)) then
197 86216 : iz1 = iz-1; iz2 = iz; exit
198 : end if
199 : end do
200 86216 : if (iz1 < 0) then
201 0 : write(*,2) 'num_logzs', num_logzs
202 0 : do iz = 1, num_logzs
203 0 : write(*,2) 'logzs(iz)', iz, logzs(iz)
204 : end do
205 0 : write(*,1) 'zlog', zlog
206 0 : write(*,*) 'confusion in do_electron_conduction'
207 0 : call mesa_error(__FILE__,__LINE__)
208 : end if
209 :
210 86216 : call get1(iz1, logK1, dlogK1_dlogRho, dlogK1_dlogT, ierr)
211 86216 : if (ierr /= 0) then
212 0 : write(*,*) 'interp failed for iz1 in do_electron_conduction', iz1, logRho, logT
213 0 : call mesa_error(__FILE__,__LINE__)
214 : end if
215 :
216 86216 : call get1(iz2, logK2, dlogK2_dlogRho, dlogK2_dlogT, ierr)
217 86216 : if (ierr /= 0) then
218 0 : write(*,*) 'interp failed for iz2 in do_electron_conduction', iz2, logRho, logT
219 0 : call mesa_error(__FILE__,__LINE__)
220 : end if
221 :
222 : ! linear interpolation in zlog
223 86216 : alfa = (zlog - logzs(iz1)) / (logzs(iz2) - logzs(iz1))
224 86216 : beta = 1d0-alfa
225 86216 : logK = alfa*logK2 + beta*logK1
226 86216 : if (clipped_logRho) then
227 22832 : dlogK_dlogRho = 0
228 : else
229 63384 : dlogK_dlogRho = alfa*dlogK2_dlogRho + beta*dlogK1_dlogRho
230 : end if
231 86216 : if (clipped_logT) then
232 0 : dlogK_dlogT = 0
233 : else
234 86216 : dlogK_dlogT = alfa*dlogK2_dlogT + beta*dlogK1_dlogT
235 : end if
236 : end if
237 :
238 : ! chi = thermal conductivity, = 10**logK (cgs units)
239 : ! conduction opacity kappa = 16*boltz_sigma*T^3 / (3*rho*chi)
240 : ! logkap = 3*logT - logRho - logK + log10(16*boltz_sigma/3)
241 :
242 86216 : logkap = 3d0*logT_in - logRho_in - logK + log10(16d0 * boltz_sigma / 3d0)
243 :
244 86216 : kap = exp10(logkap)
245 86216 : dlogkap_dlogRho = -1d0 - dlogK_dlogRho
246 86216 : dlogkap_dlogT = 3d0 - dlogK_dlogT
247 :
248 : contains
249 :
250 :
251 172432 : subroutine get1(iz, logK, dlogK_dlogRho, dlogK_dlogT, ierr)
252 : use kap_eval_support, only: Do_Kap_Interpolations
253 : integer, intent(in) :: iz
254 : real(dp), intent(out) :: logK, dlogK_dlogRho, dlogK_dlogT
255 : integer, intent(out) :: ierr
256 : logical, parameter :: dbg = .false.
257 : real(dp) :: fval, df_dx, df_dy, &
258 172432 : logRho0, logRho1, logT0, logT1
259 : integer :: i_logRho, j_logT, k
260 : include 'formats'
261 172432 : ierr = 0
262 172432 : f1(1:shift) => f_ary(1+(iz-1)*shift:iz*shift)
263 :
264 172432 : if (logRho < logRhos(2)) then
265 48472 : i_logRho = 1
266 : else
267 123960 : i_logRho = num_logRhos-1
268 2464396 : do k=2,num_logRhos-1
269 2464396 : if (logRho >= logRhos(k) .and. logRho < logRhos(k+1)) then
270 123960 : i_logRho = k; exit
271 : end if
272 : end do
273 : end if
274 172432 : logRho0 = logRhos(i_logRho)
275 172432 : logRho1 = logRhos(i_logRho+1)
276 :
277 172432 : if (logT < logTs(2)) then
278 32 : j_logT = 1
279 : else
280 172400 : j_logT = num_logTs-1
281 1995158 : do k=2,num_logTs-1
282 1995158 : if (logT >= logTs(k) .and. logT < logTs(k+1)) then
283 172400 : j_logT = k; exit
284 : end if
285 : end do
286 : end if
287 172432 : logT0 = logTs(j_logT)
288 172432 : logT1 = logTs(j_logT+1)
289 :
290 : call Do_Kap_Interpolations( &
291 : f1, num_logRhos, num_logTs, i_logRho, j_logT, logRho0, &
292 172432 : logRho, logRho1, logT0, logT, logT1, logK, dlogK_dlogRho, dlogK_dlogT)
293 172432 : if (ierr /= 0) return
294 :
295 172432 : end subroutine get1
296 :
297 :
298 : end subroutine do_electron_conduction_potekhin
299 :
300 :
301 86216 : subroutine do_electron_conduction_blouin( &
302 : zbar, logRho, logT, &
303 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
304 : use const_def, only: dp
305 : use auto_diff
306 : real(dp), intent(in) :: zbar ! average ionic charge (for electron conduction)
307 : real(dp), intent(in) :: logRho ! the density
308 : real(dp), intent(in) :: logT ! the temperature
309 : real(dp), intent(out) :: kap ! electron conduction opacity
310 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
311 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
312 : integer, intent(out) :: ierr ! 0 means AOK.
313 :
314 : ! this implements the correction formulae from Blouin et al. (2020)
315 : ! https://ui.adsabs.harvard.edu/abs/2020ApJ...899...46B/abstract
316 :
317 : real(dp), parameter :: alpha_H = -0.52d0, alpha_He = -0.46d0
318 : real(dp), parameter :: a_H = 2.0d0, a_He = 1.25d0
319 : real(dp), parameter :: b_H = 10.0d0, b_He = 2.5d0
320 : real(dp), parameter :: logRho0_H = 5.45d0, logRho0_He = 6.50d0
321 : real(dp), parameter :: logT0_H = 8.40d0, logT0_He = 8.57d0
322 : real(dp), parameter :: sigRho_H = 5.14d0, sigRho_He = 6.20d0
323 : real(dp), parameter :: sigT_H = 0.45d0, sigT_He = 0.55d0
324 :
325 : type(auto_diff_real_2var_order1) :: logRho_auto, logT_auto
326 : type(auto_diff_real_2var_order1) :: Rhostar, Tstar
327 : type(auto_diff_real_2var_order1) :: g_H, g_He, H_H, H_He
328 : type(auto_diff_real_2var_order1) :: log_correction, log_correction_H, log_correction_He
329 :
330 86216 : real(dp) :: alfa, frac_H, frac_He
331 :
332 : ! auto_diff
333 : ! var1: lnRho
334 : ! var2: lnT
335 :
336 86216 : logRho_auto = logRho
337 86216 : logRho_auto% d1val1 = iln10
338 86216 : logT_auto = logT
339 86216 : logT_auto% d1val2 = iln10
340 :
341 : ! call previous standard routines (Cassisi/Potekhin)
342 : call do_electron_conduction_potekhin( &
343 : zbar, logRho, logT, &
344 86216 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
345 :
346 : ! combined correction
347 : !
348 : ! The thermal conductivity is tabulated at Zbar = {1,2,3,4,6,...}
349 : ! and linear interpolation in logK vs logZbar is applied.
350 : !
351 : ! Therefore, we apply the Blouin+ 2020 corrections in a manner
352 : ! equivalent to individually correcting the Zbar = {1,2} tables.
353 :
354 86216 : if (Zbar <= 1d0) then ! all H
355 0 : frac_H = 1d0
356 0 : frac_He = 0d0
357 86216 : else if (Zbar <= 2d0) then ! mix H and He
358 86214 : alfa = (log10(Zbar) - log10(1d0)) / (log10(2d0) - log10(1d0))
359 86214 : frac_H = 1d0 - alfa
360 86214 : frac_He = alfa
361 2 : else if (Zbar <= 3d0) then ! mix He and no correction
362 2 : alfa = (log10(Zbar) - log10(2d0)) / (log10(3d0) - log10(2d0))
363 2 : frac_H = 0d0
364 2 : frac_He = 1d0 - alfa
365 : else ! no correction
366 : frac_H = 0d0
367 : frac_He = 0d0
368 : return
369 : end if
370 :
371 :
372 86216 : if (frac_H > 0) then
373 :
374 : ! correction for H
375 86214 : Rhostar = logRho_auto - logRho0_H
376 86214 : Tstar = logT_auto - logT0_H
377 :
378 : g_H = a_H * exp(&
379 : -pow2(Tstar*cos(alpha_H) + Rhostar*sin(alpha_H))/pow2(sigT_H) &
380 86214 : -pow2(Tstar*sin(alpha_H) - Rhostar*cos(alpha_H))/pow2(sigRho_H))
381 :
382 86214 : H_H = 0.5d0 * tanh(b_H*(g_H-0.5d0)) + 0.5d0
383 :
384 86214 : log_correction_H = -log(1d0 + g_H * H_H)
385 :
386 : else
387 :
388 2 : log_correction_H = 0d0
389 :
390 : end if
391 :
392 :
393 86216 : if (frac_He > 0) then
394 :
395 : ! correction for He
396 86216 : Rhostar = logRho_auto - logRho0_He
397 86216 : Tstar = logT_auto - logT0_He
398 :
399 : g_He = a_He * exp(&
400 : -pow2(Tstar*cos(alpha_He) + Rhostar*sin(alpha_He))/pow2(sigT_He) &
401 86216 : -pow2(Tstar*sin(alpha_He) - Rhostar*cos(alpha_He))/pow2(sigRho_He))
402 :
403 86216 : H_He = 0.5d0 * tanh(b_He*(g_He-0.5d0)) + 0.5d0
404 :
405 86216 : log_correction_He = -log(1d0 + g_He * H_He)
406 :
407 : else
408 :
409 0 : log_correction_He = 0d0
410 :
411 : end if
412 :
413 :
414 : ! blend H correction, He correction, and other correction (none)
415 86216 : log_correction = frac_H * log_correction_H + frac_He * log_correction_He
416 :
417 : ! apply correction factor
418 86216 : kap = exp(log(kap) + log_correction% val)
419 86216 : dlnkap_dlnRho = dlnkap_dlnRho + log_correction% d1val1
420 86216 : dlnkap_dlnT = dlnkap_dlnT + log_correction% d1val2
421 :
422 86216 : end subroutine do_electron_conduction_blouin
423 :
424 86216 : subroutine do_electron_conduction( &
425 : rq, zbar, logRho, logT, &
426 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
427 86216 : use kap_def, only: Kap_General_Info
428 : type (Kap_General_Info), pointer, intent(in) :: rq
429 : real(dp), intent(in) :: zbar ! average ionic charge (for electron conduction)
430 : real(dp), intent(in) :: logRho ! the density
431 : real(dp), intent(in) :: logT ! the temperature
432 : real(dp), intent(out) :: kap ! electron conduction opacity
433 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
434 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
435 : integer, intent(out) :: ierr ! 0 means AOK.
436 :
437 86216 : if (rq% use_blouin_conductive_opacities) then
438 : call do_electron_conduction_blouin( &
439 : zbar, logRho, logT, &
440 86216 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
441 : else
442 : call do_electron_conduction_potekhin( &
443 : zbar, logRho, logT, &
444 0 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
445 : end if
446 :
447 86216 : end subroutine do_electron_conduction
448 :
449 : end module condint
|
