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 suzuki_tables
21 :
22 : use const_def, only: dp, mesa_data_dir
23 : use math_lib
24 : use utils_lib, only: mesa_error, is_bad, set_nan
25 : use rates_def
26 :
27 : implicit none
28 :
29 : integer :: num_suzuki_reactions
30 :
31 : integer, pointer, dimension(:) :: & ! (num_suzuki_reactions)
32 : suzuki_lhs_nuclide_id, suzuki_rhs_nuclide_id, suzuki_reaclib_id
33 : character(len=iso_name_length), dimension(:), pointer :: &
34 : suzuki_lhs_nuclide_name, suzuki_rhs_nuclide_name ! (num_suzuki_reactions)
35 : type (integer_dict), pointer :: suzuki_reactions_dict
36 :
37 : type(table_c), dimension(:), allocatable :: suzuki_reactions_tables
38 :
39 : type, extends(weak_rate_table) :: suzuki_rate_table
40 :
41 : ! density: log_{10}(\rho Y_e):
42 : ! temperature: log_{10}T:
43 : ! chemical potential: mu
44 : ! Coulomb effect parameters:
45 : ! shift in Q-value: dQ
46 : ! shift in chemical potential: Vs
47 : ! e-capture rate or beta-decay rate: e-cap-rate: log_{10}(rate)
48 : ! neutrino energy-loss rate: nu-energy-loss: log_{10}(rate)
49 : ! gamma-ray heating rate: gamma-energy: log_{10}(rate)
50 :
51 : ! dQ, Vs, and gamma are not needed by MESA
52 :
53 : integer :: num_T
54 : real(dp), allocatable :: logTs(:)
55 :
56 : integer :: &
57 : i_capture_mu = 1, &
58 : i_capture_dQ = 2, &
59 : i_capture_Vs = 3, &
60 : i_capture_rate = 4, &
61 : i_capture_nu = 5, &
62 : i_capture_gamma = 6, &
63 : i_decay_mu = 7, &
64 : i_decay_dQ = 8, &
65 : i_decay_Vs = 9, &
66 : i_decay_rate = 10, &
67 : i_decay_nu = 11, &
68 : i_decay_gamma = 12
69 :
70 : logical :: has_decay_data
71 : logical :: has_capture_data
72 :
73 : contains
74 :
75 : procedure :: setup => setup_suzuki_table
76 : procedure :: interpolate => interpolate_suzuki_table
77 :
78 : end type suzuki_rate_table
79 :
80 : interface suzuki_rate_table
81 : module procedure new_suzuki_rate_table
82 : end interface suzuki_rate_table
83 :
84 : contains
85 :
86 :
87 69 : function new_suzuki_rate_table(logTs, lYeRhos)
88 : real(dp), intent(in), dimension(:) :: logTs, lYeRhos
89 : type(suzuki_rate_table) :: new_suzuki_rate_table
90 :
91 69 : new_suzuki_rate_table% num_T = size(logTs)
92 69 : allocate(new_suzuki_rate_table% logTs(new_suzuki_rate_table% num_T))
93 2829 : new_suzuki_rate_table% logTs = logTs
94 :
95 69 : new_suzuki_rate_table% num_lYeRho = size(lYeRhos)
96 69 : allocate(new_suzuki_rate_table% lYeRhos(new_suzuki_rate_table% num_lYeRho))
97 10626 : new_suzuki_rate_table% lYeRhos = lYeRhos
98 :
99 207 : allocate(new_suzuki_rate_table% data(1, new_suzuki_rate_table% num_T, new_suzuki_rate_table% num_lYeRho, 12))
100 69 : end function new_suzuki_rate_table
101 :
102 :
103 69 : subroutine setup_suzuki_table(table, ierr)
104 : class(suzuki_rate_table), intent(inout) :: table
105 : integer, intent(out) :: ierr
106 :
107 69 : ierr = 0
108 69 : end subroutine setup_suzuki_table
109 :
110 71 : subroutine interpolate_suzuki_table(table, T9, lYeRho, &
111 : lambda, dlambda_dlnT, dlambda_dlnRho, &
112 : Qneu, dQneu_dlnT, dQneu_dlnRho, ierr)
113 : class(suzuki_rate_table), intent(inout) :: table
114 : real(dp), intent(in) :: T9, lYeRho
115 : real(dp), intent(out) :: lambda, dlambda_dlnT, dlambda_dlnRho
116 : real(dp), intent(out) :: Qneu, dQneu_dlnT, dQneu_dlnRho
117 : integer, intent(out) :: ierr
118 :
119 : integer :: ix, jy ! target cell in the spline data
120 : real(dp) :: x0, x1 ! x0 <= xget <= x1; x0 = xs(ix), x1 = xs(ix+1)
121 : real(dp) :: y0, y1 ! y0 <= yget <= y1; y0 = ys(jy), y1 = ys(jy+1)
122 :
123 : real(dp) :: logT
124 :
125 : real(dp) :: delta_logT, dlogT, dlYeRho, delta_lYeRho, y_alfa, y_beta, x_alfa, x_beta
126 :
127 : real(dp) :: ldecay, d_ldecay_dlogT, d_ldecay_dlYeRho, &
128 : lcapture, d_lcapture_dlogT, d_lcapture_dlYeRho, &
129 : ldecay_nu, d_ldecay_nu_dlogT, d_ldecay_nu_dlYeRho, &
130 : lcapture_nu, d_lcapture_nu_dlogT, d_lcapture_nu_dlYeRho
131 :
132 : real(dp) :: decay, capture, nu, decay_nu, capture_nu
133 :
134 : logical, parameter :: dbg = .false.
135 :
136 71 : logT = log10(T9) + 9d0
137 :
138 : ! xget = logT
139 : ! yget = lYeRho
140 :
141 71 : ierr = 0
142 :
143 : ! clip small values to edge of table
144 71 : if (logT < table % logTs(1)) &
145 0 : ierr = -1 !return !logT = table % logTs(1)
146 71 : if (lYeRho < table % lYeRhos(1)) &
147 69 : ierr = -1 !return !lYeRho = table % lYeRhos(1)
148 :
149 : ! clip large values to edge of table
150 71 : if (logT > table % logTs(table % num_T)) &
151 0 : ierr = -1 !return !logT = table % logTs(table % num_T)
152 71 : if (lYeRho > table % lYeRhos(table % num_lYeRho)) &
153 0 : ierr = -1 !return !lYeRho = table % lYeRhos(table % num_lYeRho)
154 :
155 71 : if (ierr /=0) return
156 :
157 2 : call setup_for_linear_interp
158 :
159 2 : if (table % has_decay_data) then
160 : call do_linear_interp( &
161 2 : table % data(:,1:table%num_T,1:table%num_lYeRho,table%i_decay_rate), &
162 3 : ldecay, d_ldecay_dlogT, d_ldecay_dlYeRho, ierr)
163 1 : decay = exp10(ldecay)
164 : else
165 1 : decay = 0d0
166 1 : d_ldecay_dlogT = 0d0
167 1 : d_ldecay_dlYeRho = 0d0
168 : end if
169 :
170 :
171 2 : if (table % has_capture_data) then
172 : call do_linear_interp( &
173 4 : table % data(:,1:table%num_T,1:table%num_lYeRho,table%i_capture_rate), &
174 6 : lcapture, d_lcapture_dlogT, d_lcapture_dlYeRho, ierr)
175 2 : capture = exp10(lcapture)
176 : else
177 0 : capture = 0d0
178 0 : d_lcapture_dlogT = 0d0
179 0 : d_lcapture_dlYeRho = 0d0
180 : end if
181 :
182 :
183 : ! set lambda
184 2 : lambda = decay + capture
185 2 : dlambda_dlnT = decay*d_ldecay_dlogT + capture*d_lcapture_dlogT
186 2 : dlambda_dlnRho = decay*d_ldecay_dlYeRho + capture*d_lcapture_dlYeRho
187 :
188 :
189 2 : if (table % has_decay_data) then
190 : call do_linear_interp( &
191 2 : table % data(:,1:table%num_T,1:table%num_lYeRho,table%i_decay_nu), &
192 3 : ldecay_nu, d_ldecay_nu_dlogT, d_ldecay_nu_dlYeRho, ierr)
193 1 : decay_nu = exp10(ldecay_nu)
194 : else
195 1 : decay_nu = 0
196 1 : d_ldecay_nu_dlogT = 0
197 1 : d_ldecay_nu_dlYeRho = 0
198 : end if
199 :
200 2 : if (table % has_capture_data) then
201 : call do_linear_interp( &
202 4 : table % data(:,1:table%num_T,1:table%num_lYeRho,table%i_capture_nu), &
203 6 : lcapture_nu, d_lcapture_nu_dlogT, d_lcapture_nu_dlYeRho, ierr)
204 2 : capture_nu = exp10(lcapture_nu)
205 : else
206 0 : capture_nu = 0d0
207 0 : d_lcapture_nu_dlogT = 0d0
208 0 : d_lcapture_nu_dlYeRho = 0d0
209 : end if
210 :
211 :
212 : if (dbg) then
213 : write(*,*) 'logT', logT
214 : write(*,*) 'lYeRho', lYeRho
215 : write(*,*) 'ldecay', ldecay
216 : write(*,*) 'lcapture', lcapture
217 : write(*,*) 'lambda', lambda
218 : end if
219 :
220 :
221 : ! set Qneu
222 : ! be careful; you don't want to get in to the situation where 1d-99/1d-99 = 1...
223 2 : if (lambda > 1d-30) then
224 1 : nu = capture_nu + decay_nu
225 1 : Qneu = nu / lambda
226 :
227 1 : dQneu_dlnT = 0d0
228 1 : dQneu_dlnRho = 0d0
229 :
230 1 : if (nu > 1d-30) then
231 :
232 1 : if (table % has_decay_data) then
233 0 : dQneu_dlnT = dQneu_dlnT + Qneu * ((decay_nu/nu)*d_ldecay_nu_dlogT - (decay/lambda)*d_ldecay_dlogT)
234 0 : dQneu_dlnRho = dQneu_dlnRho + Qneu * ((decay_nu/nu)*d_ldecay_nu_dlYeRho - (decay/lambda)*d_ldecay_dlYeRho)
235 : end if
236 :
237 1 : if (table % has_capture_data) then
238 1 : dQneu_dlnT = dQneu_dlnT + Qneu * ((capture_nu/nu)*d_lcapture_nu_dlogT - (capture/lambda)*d_lcapture_dlogT)
239 1 : dQneu_dlnRho = dQneu_dlnRho + Qneu * ((capture_nu/nu)*d_lcapture_nu_dlYeRho - (capture/lambda)*d_lcapture_dlYeRho)
240 : end if
241 :
242 : end if
243 :
244 : else
245 1 : Qneu = 0d0
246 1 : dQneu_dlnT = 0d0
247 1 : dQneu_dlnRho = 0d0
248 : end if
249 :
250 : contains
251 :
252 8 : subroutine find_location ! set ix, jy; x is logT; y is lYeRho
253 : integer :: i, j
254 : include 'formats'
255 : ! x0 <= logT <= x1
256 2 : ix = table % num_T-1 ! since weak_num_logT is small, just do a linear search
257 22 : do i = 2, table % num_T-1
258 22 : if (logT > table% logTs(i)) cycle
259 2 : ix = i-1
260 22 : exit
261 : end do
262 :
263 : ! y0 <= lYeRho <= y1
264 2 : jy = table % num_lYeRho-1 ! since weak_num_lYeRho is small, just do a linear search
265 152 : do j = 2, table % num_lYeRho-1
266 152 : if (lYeRho > table % lYeRhos(j)) cycle
267 2 : jy = j-1
268 152 : exit
269 : end do
270 :
271 2 : x0 = table % logTs(ix)
272 2 : x1 = table % logTs(ix+1)
273 2 : y0 = table % lYeRhos(jy)
274 2 : y1 = table % lYeRhos(jy+1)
275 :
276 2 : end subroutine find_location
277 :
278 2 : subroutine setup_for_linear_interp
279 : include 'formats'
280 :
281 2 : call find_location
282 :
283 2 : dlogT = logT - x0
284 2 : delta_logT = x1 - x0
285 2 : x_beta = dlogT / delta_logT ! fraction of x1 result
286 2 : x_alfa = 1 - x_beta ! fraction of x0 result
287 2 : if (x_alfa < 0 .or. x_alfa > 1) then
288 0 : write(*,1) 'suzuki: x_alfa', x_alfa
289 0 : write(*,1) 'logT', logT
290 0 : write(*,1) 'x0', x0
291 0 : write(*,1) 'x1', x1
292 0 : call mesa_error(__FILE__,__LINE__)
293 : end if
294 :
295 2 : dlYeRho = lYeRho - y0
296 2 : delta_lYeRho = y1 - y0
297 2 : y_beta = dlYeRho / delta_lYeRho ! fraction of y1 result
298 2 : y_alfa = 1 - y_beta ! fraction of y0 result
299 2 : if (is_bad(y_alfa) .or. y_alfa < 0 .or. y_alfa > 1) then
300 0 : write(*,1) 'suzuki: y_alfa', y_alfa
301 0 : write(*,1) 'logT', logT
302 0 : write(*,1) 'x0', x0
303 0 : write(*,1) 'dlogT', dlogT
304 0 : write(*,1) 'delta_logT', delta_logT
305 0 : write(*,1) 'lYeRho', lYeRho
306 0 : write(*,1) 'y0', y0
307 0 : write(*,1) 'dlYeRho', dlYeRho
308 0 : write(*,1) 'y1', y1
309 0 : write(*,1) 'delta_lYeRho', delta_lYeRho
310 0 : write(*,1) 'y_beta', y_beta
311 : !call mesa_error(__FILE__,__LINE__,'weak setup_for_linear_interp')
312 : end if
313 :
314 : if (dbg) then
315 : write(*,2) 'logT', ix, x0, logT, x1
316 : write(*,2) 'lYeRho', jy, y0, lYeRho, y1
317 : write(*,1) 'x_alfa, x_beta', x_alfa, x_beta
318 : write(*,1) 'y_alfa, y_beta', y_alfa, y_beta
319 : write(*,'(A)')
320 : end if
321 :
322 2 : end subroutine setup_for_linear_interp
323 :
324 6 : subroutine do_linear_interp(f, fval, df_dx, df_dy, ierr)
325 : use interp_1d_lib
326 : use utils_lib, only: is_bad
327 : real(dp), dimension(:,:,:) :: f ! (4, nx, ny)
328 : real(dp), intent(out) :: fval, df_dx, df_dy
329 : integer, intent(out) :: ierr
330 :
331 : real(dp) :: fx0, fx1, fy0, fy1
332 :
333 : include 'formats'
334 :
335 6 : ierr = 0
336 :
337 6 : fx0 = y_alfa*f(1,ix,jy) + y_beta*f(1,ix,jy+1)
338 6 : fx1 = y_alfa*f(1,ix+1,jy) + y_beta*f(1,ix+1,jy+1)
339 :
340 6 : fy0 = x_alfa*f(1,ix,jy) + x_beta*f(1,ix+1,jy)
341 6 : fy1 = x_alfa*f(1,ix,jy+1) + x_beta*f(1,ix+1,jy+1)
342 :
343 6 : fval = x_alfa*fx0 + x_beta*fx1
344 6 : df_dx = (fx1 - fx0)/(x1 - x0)
345 6 : df_dy = (fy1 - fy0)/(y1 - y0)
346 :
347 6 : if (is_bad(fval)) then
348 0 : ierr = -1
349 : return
350 :
351 : write(*,1) 'x_alfa', x_alfa
352 : write(*,1) 'x_beta', x_beta
353 : write(*,1) 'fx0', fx0
354 : write(*,1) 'fx1', fx1
355 : write(*,1) 'y_alfa', y_alfa
356 : write(*,1) 'y_beta', y_beta
357 : write(*,1) 'f(1,ix,jy)', f(1,ix,jy)
358 : write(*,1) 'f(1,ix,jy+1)', f(1,ix,jy+1)
359 : !call mesa_error(__FILE__,__LINE__,'weak do_linear_interp')
360 : end if
361 :
362 6 : end subroutine do_linear_interp
363 :
364 : end subroutine interpolate_suzuki_table
365 :
366 :
367 1 : subroutine private_load_suzuki_tables(ierr)
368 :
369 : use utils_lib
370 : use forum_m, only: hdf5io_t, OPEN_FILE_RO
371 : use chem_lib, only: chem_get_iso_id
372 : use chem_def, only: iso_name_length
373 :
374 : integer, intent(out) :: ierr
375 :
376 : character (len=256) :: filename
377 : character (len=256) :: suzuki_data_dir
378 1 : type(hdf5io_t) :: hi
379 1 : type(hdf5io_t) :: hi_rxn
380 1 : character(:), allocatable :: group_names(:)
381 : integer :: num_suzuki_reactions
382 : integer :: i
383 :
384 : logical, parameter :: dbg = .false.
385 :
386 : if (dbg) write(*,*) 'private_load_suzuki_tables'
387 :
388 1 : suzuki_data_dir = trim(mesa_data_dir) // '/rates_data'
389 1 : filename = trim(suzuki_data_dir) // '/suzuki/Suzuki2016.h5'
390 : if (dbg) then
391 : write(*,'(A)')
392 : write(*,*) 'read filename <' // trim(filename) // '>'
393 : write(*,'(A)')
394 : end if
395 :
396 : ! open file (read-only)
397 :
398 1 : hi = hdf5io_t(filename, OPEN_FILE_RO)
399 :
400 : ! get a list of group names
401 :
402 70 : group_names = hi%group_names()
403 1 : num_suzuki_reactions = SIZE(group_names)
404 :
405 : ! allocate space for all this data
406 :
407 1 : call alloc
408 1 : if (failed('allocate')) return
409 :
410 : ! load the data
411 :
412 1 : nullify(suzuki_reactions_dict)
413 :
414 70 : do i = 1, num_suzuki_reactions
415 69 : hi_rxn = hdf5io_t(hi, group_names(i))
416 69 : call load_suzuki_rxn(hi_rxn, i, 'r_'//group_names(i))
417 70 : call hi_rxn% final()
418 : end do
419 :
420 : ! close file
421 :
422 1 : call hi% final()
423 :
424 : ! set up reaction dictionary
425 :
426 1 : call integer_dict_create_hash(suzuki_reactions_dict, ierr)
427 1 : if (failed('integer_dict_create_hash')) return
428 :
429 : ! pre-construct interpolants
430 :
431 70 : do i = 1, num_suzuki_reactions
432 0 : associate(t => suzuki_reactions_tables(i) % t)
433 69 : if (ierr == 0) call t% setup(ierr)
434 : end associate
435 70 : if (failed('setup')) return
436 : end do
437 :
438 1 : if (dbg) write(*,*) 'finished load_suzuki_tables'
439 :
440 : contains
441 :
442 69 : subroutine load_suzuki_rxn(hi, rxn_idx, rxn_name)
443 :
444 1 : use weak_support, only : parse_weak_rate_name
445 :
446 : type(hdf5io_t), intent(inout) :: hi
447 : integer, intent(in) :: rxn_idx
448 : character(*), intent(in) :: rxn_name
449 :
450 69 : type(hdf5io_t) :: hi_data
451 69 : real(dp), allocatable, dimension(:) :: logTs, lYeRhos
452 69 : type(suzuki_rate_table) :: table
453 : character(len=iso_name_length) :: lhs, rhs
454 : character(len=2*iso_name_length+1) :: key
455 :
456 : ! parse the name
457 :
458 69 : call parse_weak_rate_name(rxn_name, lhs, rhs, ierr)
459 : if (dbg) write(*,*) 'parse_weak_rate_name gives ', trim(lhs), ' ', trim(rhs), ierr
460 :
461 4830 : suzuki_lhs_nuclide_id = chem_get_iso_id(lhs)
462 69 : suzuki_lhs_nuclide_name(rxn_idx) = lhs
463 4830 : suzuki_rhs_nuclide_id = chem_get_iso_id(rhs)
464 69 : suzuki_rhs_nuclide_name(rxn_idx) = rhs
465 69 : call create_weak_dict_key(lhs, rhs, key)
466 69 : call integer_dict_define(suzuki_reactions_dict, key, rxn_idx, ierr)
467 69 : if (failed('integer_dict_define')) return
468 :
469 : ! read table axes
470 :
471 69 : call hi% alloc_read_dset('logTs', logTs)
472 : if (dbg) write(*,*) "num logTs", SIZE(logTs)
473 :
474 69 : call hi% alloc_read_dset('lYeRhos', lYeRhos)
475 : if (dbg) write(*,*) "num lYeRhos", SIZE(lYeRhos)
476 :
477 69 : table = suzuki_rate_table(logTs, lYeRhos)
478 69 : call set_nan(table% data)
479 :
480 : ! read capture data
481 :
482 69 : table% has_capture_data = hi% group_exists('capture')
483 :
484 69 : if (table% has_capture_data) then
485 :
486 : if (dbg) write(*,*) "found capture group; reading..."
487 :
488 63 : hi_data = hdf5io_t(hi, 'capture')
489 :
490 63 : call hi_data% read_dset('mu', table% data(1,:,:,table%i_capture_mu))
491 63 : call hi_data% read_dset('dQ', table% data(1,:,:,table%i_capture_dQ))
492 63 : call hi_data% read_dset('Vs', table% data(1,:,:,table%i_capture_Vs))
493 63 : call hi_data% read_dset('rate', table% data(1,:,:,table%i_capture_rate))
494 63 : call hi_data% read_dset('nu', table% data(1,:,:,table%i_capture_nu))
495 63 : call hi_data% read_dset('gamma', table% data(1,:,:,table%i_capture_gamma))
496 :
497 63 : call hi_data% final()
498 :
499 : end if
500 :
501 : ! read decay data
502 :
503 69 : table% has_decay_data = hi% group_exists('decay')
504 :
505 69 : if (table% has_decay_data) then
506 :
507 : if (dbg) write(*,*) "found decay group; reading..."
508 :
509 40 : hi_data = hdf5io_t(hi, 'decay')
510 :
511 40 : call hi_data% read_dset('mu', table% data(1,:,:,table%i_decay_mu))
512 40 : call hi_data% read_dset('dQ', table% data(1,:,:,table%i_decay_dQ))
513 40 : call hi_data% read_dset('Vs', table% data(1,:,:,table%i_decay_Vs))
514 40 : call hi_data% read_dset('rate', table% data(1,:,:,table%i_decay_rate))
515 40 : call hi_data% read_dset('nu', table% data(1,:,:,table%i_decay_nu))
516 40 : call hi_data% read_dset('gamma', table% data(1,:,:,table%i_decay_gamma))
517 :
518 40 : call hi_data% final()
519 :
520 : end if
521 :
522 : ! assign the table
523 :
524 69 : allocate(suzuki_reactions_tables(rxn_idx)% t, source=table)
525 :
526 : ! finish
527 :
528 69 : end subroutine load_suzuki_rxn
529 :
530 1 : subroutine alloc
531 :
532 : allocate( &
533 : suzuki_reaclib_id(num_suzuki_reactions), &
534 : suzuki_lhs_nuclide_name(num_suzuki_reactions), &
535 : suzuki_rhs_nuclide_name(num_suzuki_reactions), &
536 : suzuki_lhs_nuclide_id(num_suzuki_reactions), &
537 : suzuki_rhs_nuclide_id(num_suzuki_reactions), &
538 : suzuki_reactions_tables(num_suzuki_reactions), &
539 70 : stat=ierr)
540 :
541 1 : end subroutine alloc
542 :
543 140 : logical function failed(str)
544 : character (len=*) :: str
545 140 : failed = (ierr /= 0)
546 140 : if (failed) then
547 0 : write(*,*) 'failed: ' // trim(str)
548 : end if
549 140 : end function failed
550 :
551 : end subroutine private_load_suzuki_tables
552 :
553 69 : end module suzuki_tables
|
