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 2 : subroutine setup
35 : use chem_def
36 :
37 : integer :: ierr
38 : character(len=32) :: my_mesa_dir
39 :
40 : include 'formats'
41 :
42 2 : ierr = 0
43 :
44 2 : my_mesa_dir = '../..'
45 2 : call const_init(my_mesa_dir, ierr)
46 2 : if (ierr /= 0) then
47 0 : write (*, *) 'const_init failed'
48 0 : call mesa_error(__FILE__, __LINE__)
49 : end if
50 :
51 2 : call math_init()
52 :
53 2 : call chem_init('isotopes.data', ierr)
54 2 : 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 2 : '', ierr)
65 2 : if (ierr /= 0) then
66 0 : write (*, *) 'rates_init failed'
67 0 : call mesa_error(__FILE__, __LINE__)
68 : end if
69 :
70 2 : call rates_warning_init(.true., 10d0)
71 :
72 2 : call read_raw_rates_records(ierr)
73 2 : if (ierr /= 0) then
74 0 : write (*, *) 'read_raw_rates_records failed'
75 0 : call mesa_error(__FILE__, __LINE__)
76 : end if
77 :
78 2 : end subroutine setup
79 :
80 2 : subroutine do_test_rates()
81 :
82 : integer :: ierr
83 : type(T_Factors), target :: tf_rec
84 : type(T_Factors), pointer :: tf
85 2 : real(dp) :: logT, temp
86 : integer :: i, t
87 :
88 : integer :: nrates_to_eval
89 2 : integer, allocatable :: irs(:)
90 2 : real(dp), allocatable :: raw_rates(:)
91 20 : real(dp), dimension(9) :: temps
92 :
93 : logical, parameter :: dbg = .false.
94 :
95 : include 'formats'
96 :
97 2 : write (*, '(A)')
98 :
99 2 : temps = [6.0d0, 6.5d0, 7.0d0, 7.5d0, 8.0d0, 8.5d0, 9.0d0, 9.5d0, 10.0d0]
100 :
101 2 : 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 2 : nrates_to_eval = num_predefined_reactions
115 2 : allocate (irs(nrates_to_eval), raw_rates(nrates_to_eval))
116 616 : do i = 1, nrates_to_eval
117 616 : irs(i) = i
118 : end do
119 :
120 : end if
121 :
122 20 : do t = 1, size(temps)
123 18 : logT = temps(t)
124 18 : temp = exp10(logT)
125 18 : call eval_tfactors(tf, logT, temp)
126 :
127 18 : write (*, 1) 'logT', logT
128 18 : write (*, 1) 'temp', temp
129 18 : write (*, '(A)')
130 :
131 5544 : raw_rates = missing_value
132 :
133 18 : call get_raw_rates(nrates_to_eval, irs, temp, tf, raw_rates, ierr)
134 18 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
135 :
136 5544 : do i = 1, nrates_to_eval
137 5526 : 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 5544 : write (*, 1) trim(reaction_Name(irs(i))), raw_rates(i)
142 : end do
143 38 : write (*, '(A)')
144 : end do
145 :
146 2 : write (*, *) 'done'
147 2 : write (*, '(A)')
148 :
149 4 : 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 0 : 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 0 : 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 2 : subroutine do_test_FL_epsnuc_3alf
217 2 : real(dp) :: T ! temperature
218 2 : real(dp) :: Rho ! density
219 2 : real(dp) :: Y ! helium mass fraction
220 2 : real(dp) :: UE ! electron molecular weight
221 2 : real(dp) :: eps_nuc ! eps_nuc in ergs/g/sec
222 2 : real(dp) :: deps_nuc_dT ! partial wrt temperature
223 2 : real(dp) :: deps_nuc_dRho ! partial wrt density
224 : include 'formats'
225 2 : T = 1d7
226 2 : Rho = 1d10
227 2 : Y = 1
228 2 : UE = 2
229 2 : call eval_FL_epsnuc_3alf(T, Rho, Y, UE, eps_nuc, deps_nuc_dT, deps_nuc_dRho)
230 2 : write (*, 1) 'FL_epsnuc_3alf', eps_nuc
231 2 : write (*, '(A)')
232 2 : end subroutine do_test_FL_epsnuc_3alf
233 :
234 2 : subroutine do_test_rate_table
235 : integer :: ierr
236 : type(T_Factors), target :: tf_rec
237 : type(T_Factors), pointer :: tf
238 2 : real(dp) :: logT, temp, raw_rate
239 : integer :: ir
240 : logical, parameter :: dbg = .false.
241 :
242 : include 'formats'
243 :
244 2 : write (*, '(A)')
245 2 : write (*, *) 'do_test_rate_table'
246 :
247 2 : tf => tf_rec
248 :
249 2 : temp = 9.0d8
250 2 : logT = log10(temp)
251 2 : call eval_tfactors(tf, logT, temp)
252 :
253 2 : write (*, 1) 'logT', logT
254 2 : write (*, 1) 'temp', temp
255 2 : write (*, '(A)')
256 :
257 2 : ir = rates_reaction_id('r3')
258 2 : call run1
259 :
260 2 : write (*, *) 'done'
261 2 : write (*, '(A)')
262 :
263 : contains
264 :
265 2 : subroutine run1
266 : include 'formats'
267 2 : call get_raw_rate(ir, temp, tf, raw_rate, ierr)
268 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
269 2 : write (*, 1) trim(reaction_Name(ir)), raw_rate
270 2 : write (*, '(A)')
271 2 : end subroutine run1
272 :
273 : end subroutine do_test_rate_table
274 :
275 0 : subroutine do_test2_FL_epsnuc_3alf
276 0 : real(dp) :: T ! temperature
277 0 : real(dp) :: Rho ! density
278 0 : real(dp) :: Y ! helium mass fraction
279 0 : real(dp) :: UE ! electron molecular weight
280 0 : real(dp) :: eps_nuc1, eps_nuc2 ! eps_nuc in ergs/g/sec
281 0 : real(dp) :: deps_nuc_dT ! partial wrt temperature
282 0 : real(dp) :: deps_nuc_dRho ! partial wrt density
283 0 : 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 2 : subroutine teardown
306 2 : call rates_shutdown
307 2 : end subroutine teardown
308 :
309 : end module test_rates_support
310 :
311 2 : program test_rates
312 :
313 2 : use test_screen
314 : use test_weak
315 : use test_ecapture
316 : use test_rates_support
317 :
318 : implicit none
319 :
320 2 : call setup
321 :
322 : !call do_test_rates(rates_JR_if_available); stop
323 : !call test1; stop
324 :
325 2 : call do_test_screen
326 :
327 2 : call do_test_weak
328 :
329 2 : call do_test_ecapture
330 :
331 2 : call do_test_rates()
332 2 : call do_test_FL_epsnuc_3alf()
333 2 : call do_test_rate_table
334 :
335 2 : call teardown
336 :
337 2 : end program test_rates
|
