Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2011-2020 Bill Paxton & 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 test_rates_support
21 :
22 : use rates_def
23 : use rates_lib
24 : use chem_lib
25 : use const_def, only: missing_value
26 : use const_lib, only: const_init
27 : use math_lib
28 : use utils_lib, only: mesa_error
29 :
30 : implicit none
31 :
32 : contains
33 :
34 3 : subroutine setup
35 : use chem_def
36 :
37 : integer :: ierr
38 : character(len=32) :: my_mesa_dir
39 :
40 : include 'formats'
41 :
42 : ierr = 0
43 :
44 1 : my_mesa_dir = '../..'
45 1 : call const_init(my_mesa_dir, ierr)
46 1 : if (ierr /= 0) then
47 0 : write (*, *) 'const_init failed'
48 0 : call mesa_error(__FILE__, __LINE__)
49 : end if
50 :
51 1 : call math_init()
52 :
53 1 : call chem_init('isotopes.data', ierr)
54 1 : if (ierr /= 0) then
55 0 : write (*, *) 'chem_init failed'
56 0 : call mesa_error(__FILE__, __LINE__)
57 : end if
58 :
59 : ! use special weak reaction data in test directory
60 :
61 : call rates_init('reactions.list', '', 'rate_tables', &
62 : .true., &
63 : .true., 'test_special.states', 'test_special.transitions', &
64 1 : '', ierr)
65 1 : if (ierr /= 0) then
66 0 : write (*, *) 'rates_init failed'
67 0 : call mesa_error(__FILE__, __LINE__)
68 : end if
69 :
70 1 : call rates_warning_init(.true., 10d0)
71 :
72 1 : call read_raw_rates_records(ierr)
73 1 : if (ierr /= 0) then
74 0 : write (*, *) 'read_raw_rates_records failed'
75 0 : call mesa_error(__FILE__, __LINE__)
76 : end if
77 :
78 1 : end subroutine setup
79 :
80 1 : subroutine do_test_rates()
81 :
82 : integer :: ierr
83 : type(T_Factors), target :: tf_rec
84 : type(T_Factors), pointer :: tf
85 : real(dp) :: logT, temp
86 : integer :: i, t
87 :
88 : integer :: nrates_to_eval
89 : integer, allocatable :: irs(:)
90 : real(dp), allocatable :: raw_rates(:)
91 : real(dp), dimension(9) :: temps
92 :
93 : logical, parameter :: dbg = .false.
94 :
95 : include 'formats'
96 :
97 1 : write (*, '(A)')
98 :
99 1 : temps = [6.0d0, 6.5d0, 7.0d0, 7.5d0, 8.0d0, 8.5d0, 9.0d0, 9.5d0, 10.0d0]
100 :
101 1 : tf => tf_rec
102 :
103 : if (dbg) then
104 :
105 : nrates_to_eval = 1
106 : allocate (irs(nrates_to_eval), raw_rates(nrates_to_eval))
107 :
108 : irs(1:nrates_to_eval) = [ &
109 : ir_s32_ga_si28 &
110 : ]
111 :
112 : else
113 :
114 1 : nrates_to_eval = num_predefined_reactions
115 1 : allocate (irs(nrates_to_eval), raw_rates(nrates_to_eval))
116 308 : do i = 1, nrates_to_eval
117 308 : irs(i) = i
118 : end do
119 :
120 : end if
121 :
122 10 : do t = 1, size(temps)
123 9 : logT = temps(t)
124 9 : temp = exp10(logT)
125 9 : call eval_tfactors(tf, logT, temp)
126 :
127 9 : write (*, 1) 'logT', logT
128 9 : write (*, 1) 'temp', temp
129 9 : write (*, '(A)')
130 :
131 2772 : raw_rates = missing_value
132 :
133 9 : call get_raw_rates(nrates_to_eval, irs, temp, tf, raw_rates, ierr)
134 9 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
135 :
136 2772 : do i = 1, nrates_to_eval
137 2763 : if (raw_rates(i) == missing_value) then
138 0 : write (*, *) 'missing value for '//trim(reaction_Name(irs(i)))
139 0 : call mesa_error(__FILE__, __LINE__)
140 : end if
141 2772 : write (*, 1) trim(reaction_Name(irs(i))), raw_rates(i)
142 : end do
143 19 : write (*, '(A)')
144 : end do
145 :
146 1 : write (*, *) 'done'
147 1 : write (*, '(A)')
148 :
149 1 : end subroutine do_test_rates
150 :
151 0 : subroutine test1
152 : integer :: ierr
153 : type(T_Factors), target :: tf_rec
154 : type(T_Factors), pointer :: tf
155 : real(dp) :: logT, temp, raw_rate, raw_rate1, raw_rate2
156 : integer :: ir
157 : logical, parameter :: dbg = .false.
158 :
159 : include 'formats'
160 :
161 0 : write (*, '(A)')
162 0 : write (*, *) 'test1'
163 :
164 0 : tf => tf_rec
165 :
166 0 : temp = 3d9
167 0 : logT = log10(temp)
168 0 : call eval_tfactors(tf, logT, temp)
169 :
170 0 : write (*, 1) 'logT', logT
171 0 : write (*, 1) 'temp', temp
172 0 : write (*, '(A)')
173 :
174 0 : ir = rates_reaction_id('r_ni56_wk_co56')
175 0 : if (ir == 0) then
176 0 : write (*, *) 'failed to find rate id'
177 0 : call mesa_error(__FILE__, __LINE__)
178 : end if
179 :
180 0 : call run1
181 0 : raw_rate1 = raw_rate
182 :
183 0 : temp = 3.000000001d9
184 0 : logT = log10(temp)
185 0 : call eval_tfactors(tf, logT, temp)
186 :
187 0 : write (*, 1) 'logT', logT
188 0 : write (*, 1) 'temp', temp
189 0 : write (*, 1) 'raw_rate1', raw_rate1
190 0 : write (*, '(A)')
191 :
192 0 : stop
193 :
194 : ir = rates_reaction_id('r_s32_ga_si28')
195 : call run1
196 : raw_rate2 = raw_rate
197 :
198 : write (*, 1) 'raw_rate2', raw_rate2
199 : write (*, 1) 'raw_rate2-raw_rate1', raw_rate2 - raw_rate1
200 :
201 : write (*, *) 'done'
202 : write (*, '(A)')
203 :
204 : contains
205 :
206 0 : subroutine run1
207 : include 'formats'
208 : call get_raw_rate(ir, temp, tf, raw_rate, ierr)
209 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
210 0 : write (*, 1) trim(reaction_Name(ir)), raw_rate
211 0 : write (*, '(A)')
212 0 : end subroutine run1
213 :
214 : end subroutine test1
215 :
216 1 : subroutine do_test_FL_epsnuc_3alf
217 : real(dp) :: T ! temperature
218 : real(dp) :: Rho ! density
219 : real(dp) :: Y ! helium mass fraction
220 : real(dp) :: UE ! electron molecular weight
221 : real(dp) :: eps_nuc ! eps_nuc in ergs/g/sec
222 : real(dp) :: deps_nuc_dT ! partial wrt temperature
223 : real(dp) :: deps_nuc_dRho ! partial wrt density
224 : include 'formats'
225 1 : T = 1d7
226 1 : Rho = 1d10
227 1 : Y = 1
228 1 : UE = 2
229 1 : call eval_FL_epsnuc_3alf(T, Rho, Y, UE, eps_nuc, deps_nuc_dT, deps_nuc_dRho)
230 1 : write (*, 1) 'FL_epsnuc_3alf', eps_nuc
231 1 : write (*, '(A)')
232 1 : end subroutine do_test_FL_epsnuc_3alf
233 :
234 1 : subroutine do_test_rate_table
235 : integer :: ierr
236 : type(T_Factors), target :: tf_rec
237 : type(T_Factors), pointer :: tf
238 : real(dp) :: logT, temp, raw_rate
239 : integer :: ir
240 : logical, parameter :: dbg = .false.
241 :
242 : include 'formats'
243 :
244 1 : write (*, '(A)')
245 1 : write (*, *) 'do_test_rate_table'
246 :
247 1 : tf => tf_rec
248 :
249 1 : temp = 9.0d8
250 1 : logT = log10(temp)
251 1 : call eval_tfactors(tf, logT, temp)
252 :
253 1 : write (*, 1) 'logT', logT
254 1 : write (*, 1) 'temp', temp
255 1 : write (*, '(A)')
256 :
257 1 : ir = rates_reaction_id('r3')
258 1 : call run1
259 :
260 1 : write (*, *) 'done'
261 1 : write (*, '(A)')
262 :
263 : contains
264 :
265 1 : subroutine run1
266 : include 'formats'
267 : call get_raw_rate(ir, temp, tf, raw_rate, ierr)
268 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
269 1 : write (*, 1) trim(reaction_Name(ir)), raw_rate
270 1 : write (*, '(A)')
271 1 : end subroutine run1
272 :
273 : end subroutine do_test_rate_table
274 :
275 1 : subroutine do_test_reverse_coefficients
276 : use chem_def, only: nuclide_data
277 : use const_def, only: pi, kB=>boltzm, NA=>avo, hbar
278 :
279 : integer, parameter :: nchecks = 6
280 :
281 : integer :: i
282 : integer, dimension(nchecks) :: chapters, inverse_exps
283 : type(nuclide_data), pointer :: nuclides
284 : real(dp) :: fac
285 :
286 : include 'formats'
287 :
288 : ! Check representative reverse rates with Ni-No > 0, = 0, and < 0.
289 : ! Include both 1 -> 1 and 2 -> 2 cases for Ni-No = 0.
290 : ! For each case, use the first non-weak reaclib reaction with the requested
291 : ! chapter and inverse_exp. The printed handle shows which reaction was used.
292 : ! If this test network does not include a representative for one case,
293 : ! print a skip line and continue.
294 : ! The mass ratio always carries a fixed 3/2 power. Ni-No enters through fac and inverse_exp.
295 1 : chapters = [r_one_one, r_three_one, r_two_one, r_two_two, r_one_two, r_one_three]
296 1 : inverse_exps = [0, 2, 1, 0, -1, -2]
297 :
298 1 : nuclides => reaclib_rates% nuclides
299 1 : if (.not.associated(nuclides)) call mesa_error(__FILE__, __LINE__)
300 :
301 1 : fac = pow(1d9*kB/(2d0*pi*hbar*hbar*NA),1.5d0)/NA
302 :
303 1 : write (*, '(A)')
304 1 : write (*, *) 'do_test_reverse_coefficients'
305 1 : write (*, *) 'inverse_coefficients(1) = log[(m_in/m_out)^(3/2)*(g_in/g_out)*fac^(Ni-No)]'
306 : write (*, '(a32, 2x, a3, 2x, a3, 2x, a5, 2x, a26, 2x, a)') &
307 1 : 'rate', 'Ni', 'No', 'Ni-No', 'inverse_coefficients(1)', 'note'
308 :
309 7 : do i = 1, nchecks
310 7 : call check1(chapters(i), inverse_exps(i))
311 : end do
312 :
313 1 : write (*, *) 'done'
314 1 : write (*, '(A)')
315 :
316 : contains
317 :
318 10 : subroutine check1(chapter, expected_inverse_exp)
319 : integer, intent(in) :: chapter, expected_inverse_exp
320 :
321 : integer :: ierr, i, j, lo, hi
322 : integer :: Ni, No, Nt, inverse_exp
323 : integer, dimension(max_species_per_reaction) :: ps
324 : character(len=max_id_length) :: handle
325 : real(dp), dimension(max_species_per_reaction) :: g, mass
326 : real(dp) :: tmp, expected_log_coeff, stored_log_coeff, tol
327 :
328 : include 'formats'
329 :
330 : ! Select the first non-weak reaclib reaction in this chapter with the
331 : ! requested inverse_exp. This keeps the test independent of branch-specific
332 : ! reaction handles while still checking one representative for each Ni, No case.
333 6 : lo = 0
334 197524 : do i = 1, reaclib_rates% nreactions
335 197522 : if (reaclib_rates% chapter(i) /= chapter) cycle
336 10004 : if (reaclib_rates% reaction_flag(i) == 'w' .or. reaclib_rates% reaction_flag(i) == 'e') cycle
337 4 : if (reaclib_rates% inverse_exp(i) /= expected_inverse_exp) cycle
338 4 : lo = i
339 197524 : exit
340 : end do
341 :
342 6 : if (lo <= 0) then
343 2 : Ni = Nin(chapter)
344 2 : No = Nout(chapter)
345 : write (*, '(a32, 2x, i3, 2x, i3, 2x, i5, 2x, a26, 2x, a)') &
346 2 : '(none)', Ni, No, expected_inverse_exp, '', 'skipped'
347 2 : return
348 : end if
349 :
350 4 : handle = trim(reaclib_rates% reaction_handle(lo))
351 4 : call reaclib_indices_for_reaction(handle, reaclib_rates, lo, hi, ierr)
352 4 : if (ierr /= 0 .or. lo <= 0 .or. hi < lo) then
353 0 : write (*, *) 'failed to find reaction '//trim(handle)
354 0 : call mesa_error(__FILE__, __LINE__)
355 : end if
356 :
357 4 : Ni = Nin(reaclib_rates% chapter(lo))
358 4 : No = Nout(reaclib_rates% chapter(lo))
359 4 : Nt = Ni + No
360 4 : inverse_exp = Ni - No
361 :
362 18 : ps(1:Nt) = reaclib_rates% pspecies(1:Nt, lo)
363 18 : g(1:Nt) = 2d0*nuclides% spin(ps(1:Nt)) + 1d0
364 18 : mass(1:Nt) = nuclides% W(ps(1:Nt))
365 :
366 18 : tmp = product(mass(1:Ni))/product(mass(Ni+1:Nt))
367 18 : expected_log_coeff = log(pow(tmp, 1.5d0)*(product(g(1:Ni))/product(g(Ni+1:Nt))))
368 4 : if (inverse_exp /= 0) expected_log_coeff = expected_log_coeff + dble(inverse_exp)*log(fac)
369 :
370 4 : tol = 1d-12*max(1d0, abs(expected_log_coeff))
371 :
372 10 : do j = lo, hi
373 6 : if (reaclib_rates% inverse_exp(j) /= inverse_exp) then
374 0 : write (*, *) 'bad inverse_exp for '//trim(handle)
375 0 : call mesa_error(__FILE__, __LINE__)
376 : end if
377 :
378 6 : stored_log_coeff = reaclib_rates% inverse_coefficients(1, j)
379 10 : if (abs(stored_log_coeff - expected_log_coeff) > tol) then
380 0 : write (*, *) 'bad inverse coefficient for '//trim(handle)
381 0 : write (*, 1) 'stored', stored_log_coeff
382 0 : write (*, 1) 'expected', expected_log_coeff
383 0 : call mesa_error(__FILE__, __LINE__)
384 : end if
385 : end do
386 :
387 : write (*, '(a32, 2x, i3, 2x, i3, 2x, i5, 2x, 1pe26.16, 2x, a)') &
388 4 : trim(handle), Ni, No, inverse_exp, expected_log_coeff, ''
389 :
390 : end subroutine check1
391 :
392 : end subroutine do_test_reverse_coefficients
393 :
394 0 : subroutine do_test2_FL_epsnuc_3alf
395 : real(dp) :: T ! temperature
396 : real(dp) :: Rho ! density
397 : real(dp) :: Y ! helium mass fraction
398 : real(dp) :: UE ! electron molecular weight
399 : real(dp) :: eps_nuc1, eps_nuc2 ! eps_nuc in ergs/g/sec
400 : real(dp) :: deps_nuc_dT ! partial wrt temperature
401 : real(dp) :: deps_nuc_dRho ! partial wrt density
402 : real(dp) :: dT, dRho
403 : include 'formats'
404 0 : T = 7.9432823472428218d+07
405 0 : dT = T*1d-8
406 0 : Rho = 3.1622776601683793d+09
407 0 : dRho = Rho*1d-8
408 0 : Y = 1
409 0 : UE = 2
410 0 : call eval_FL_epsnuc_3alf(T, Rho + dRho, Y, UE, eps_nuc1, deps_nuc_dT, deps_nuc_dRho)
411 0 : write (*, 1) 'FL_epsnuc_3alf 1', eps_nuc1
412 0 : write (*, '(A)')
413 0 : call eval_FL_epsnuc_3alf(T, Rho, Y, UE, eps_nuc2, deps_nuc_dT, deps_nuc_dRho)
414 0 : write (*, 1) 'FL_epsnuc_3alf 2', eps_nuc2
415 0 : write (*, '(A)')
416 0 : write (*, 1) 'analytic deps_nuc_dRho', deps_nuc_dRho
417 0 : write (*, 1) 'numerical deps_nuc_dRho', (eps_nuc1 - eps_nuc2)/dRho
418 0 : write (*, '(A)')
419 0 : write (*, 1) 'analytic dlneps_nuc_dlnRho', deps_nuc_dRho*Rho/eps_nuc2
420 0 : write (*, 1) 'numerical dlneps_nuc_dlnRho', (eps_nuc1 - eps_nuc2)/dRho*Rho/eps_nuc2
421 0 : write (*, '(A)')
422 0 : end subroutine do_test2_FL_epsnuc_3alf
423 :
424 1 : subroutine teardown
425 1 : call rates_shutdown
426 1 : end subroutine teardown
427 :
428 : end module test_rates_support
429 :
430 1 : program test_rates
431 :
432 1 : use test_screen
433 : use test_weak
434 : use test_ecapture
435 : use test_rates_support
436 :
437 : implicit none
438 :
439 1 : call setup
440 :
441 : !call do_test_rates(rates_JR_if_available); stop
442 : !call test1; stop
443 :
444 1 : call do_test_screen
445 :
446 1 : call do_test_weak
447 :
448 1 : call do_test_ecapture
449 :
450 1 : call do_test_rates()
451 1 : call do_test_FL_epsnuc_3alf()
452 1 : call do_test_rate_table
453 1 : call do_test_reverse_coefficients()
454 :
455 1 : call teardown
456 :
457 1 : end program test_rates
|
