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 1 : 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 3 : 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 0 : 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 2 : 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 0 : subroutine do_test2_FL_epsnuc_3alf
276 : real(dp) :: T ! temperature
277 : real(dp) :: Rho ! density
278 : real(dp) :: Y ! helium mass fraction
279 : real(dp) :: UE ! electron molecular weight
280 : real(dp) :: eps_nuc1, eps_nuc2 ! eps_nuc in ergs/g/sec
281 : real(dp) :: deps_nuc_dT ! partial wrt temperature
282 : real(dp) :: deps_nuc_dRho ! partial wrt density
283 : real(dp) :: dT, dRho
284 : include 'formats'
285 0 : T = 7.9432823472428218d+07
286 0 : dT = T*1d-8
287 0 : Rho = 3.1622776601683793d+09
288 0 : dRho = Rho*1d-8
289 0 : Y = 1
290 0 : UE = 2
291 0 : call eval_FL_epsnuc_3alf(T, Rho + dRho, Y, UE, eps_nuc1, deps_nuc_dT, deps_nuc_dRho)
292 0 : write (*, 1) 'FL_epsnuc_3alf 1', eps_nuc1
293 0 : write (*, '(A)')
294 0 : call eval_FL_epsnuc_3alf(T, Rho, Y, UE, eps_nuc2, deps_nuc_dT, deps_nuc_dRho)
295 0 : write (*, 1) 'FL_epsnuc_3alf 2', eps_nuc2
296 0 : write (*, '(A)')
297 0 : write (*, 1) 'analytic deps_nuc_dRho', deps_nuc_dRho
298 0 : write (*, 1) 'numerical deps_nuc_dRho', (eps_nuc1 - eps_nuc2)/dRho
299 0 : write (*, '(A)')
300 0 : write (*, 1) 'analytic dlneps_nuc_dlnRho', deps_nuc_dRho*Rho/eps_nuc2
301 0 : write (*, 1) 'numerical dlneps_nuc_dlnRho', (eps_nuc1 - eps_nuc2)/dRho*Rho/eps_nuc2
302 0 : write (*, '(A)')
303 0 : end subroutine do_test2_FL_epsnuc_3alf
304 :
305 1 : subroutine teardown
306 1 : call rates_shutdown
307 0 : end subroutine teardown
308 :
309 : end module test_rates_support
310 :
311 1 : program test_rates
312 :
313 1 : use test_screen
314 : use test_weak
315 : use test_ecapture
316 : use test_rates_support
317 :
318 : implicit none
319 :
320 1 : call setup
321 :
322 : !call do_test_rates(rates_JR_if_available); stop
323 : !call test1; stop
324 :
325 1 : call do_test_screen
326 :
327 1 : call do_test_weak
328 :
329 1 : call do_test_ecapture
330 :
331 1 : call do_test_rates()
332 1 : call do_test_FL_epsnuc_3alf()
333 1 : call do_test_rate_table
334 :
335 : call teardown
336 :
337 1 : end program test_rates
|
