Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2009-2019 Bill Paxton, Frank Timmes & 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_net_support
21 :
22 : use chem_def
23 : use chem_lib, only: chem_init
24 : use math_lib
25 : use net_def
26 : use net_lib
27 : use const_def
28 : use rates_def
29 : use test_net_do_one
30 : use utils_lib, only: mesa_error
31 :
32 : implicit none
33 :
34 : integer, parameter :: max_files = 20, max_cnt = 100000
35 :
36 : character(len=256) :: cache_suffix
37 : integer :: num_reactions
38 :
39 : integer, dimension(:), pointer :: net_iso, chem_id, isos_to_show
40 :
41 : integer, pointer :: reaction_table(:)
42 : integer, pointer :: rates_to_show(:)
43 :
44 : real(dp), dimension(:), pointer :: rho_vector, T_vector
45 :
46 : integer :: nrates_to_show, nisos_to_show, net_handle
47 :
48 : contains
49 :
50 14 : subroutine do_test_net(do_plots, symbolic)
51 : logical, intent(in) :: do_plots, symbolic
52 14 : call set_composition(species, xin)
53 14 : eta = 0
54 14 : if (do_plots) then
55 0 : call Create_Plot_Files(species, xin)
56 : else
57 14 : call Do_One_Net(symbolic)
58 : end if
59 14 : end subroutine do_test_net
60 :
61 2 : subroutine load_libs
62 : use const_lib, only: const_init
63 : use const_def, only: mesa_dir
64 : use chem_lib
65 : use rates_lib, only: rates_init, rates_warning_init
66 : integer :: ierr
67 : character(len=64) :: my_mesa_dir
68 :
69 2 : my_mesa_dir = '../..'
70 2 : call const_init(my_mesa_dir, ierr)
71 2 : if (ierr /= 0) then
72 0 : write (*, *) 'const_init failed'
73 0 : call mesa_error(__FILE__, __LINE__)
74 : end if
75 :
76 2 : call math_init()
77 :
78 2 : ierr = 0
79 2 : call chem_init('isotopes.data', ierr)
80 2 : if (ierr /= 0) then
81 0 : write (*, *) 'chem_init failed'
82 0 : call mesa_error(__FILE__, __LINE__)
83 : end if
84 :
85 : call rates_init('reactions.list', '', 'rate_tables', .false., .false., &
86 2 : '', '', '', ierr)
87 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
88 :
89 2 : call rates_warning_init(.true., 10d0)
90 :
91 2 : end subroutine load_libs
92 :
93 14 : subroutine test_net_setup(net_file_in)
94 : character(len=*), intent(in) :: net_file_in
95 : type(Net_General_Info), pointer :: g
96 :
97 : integer :: info, ierr
98 :
99 : include 'formats'
100 :
101 14 : net_file = net_file_in
102 :
103 14 : allocate (net_iso(num_chem_isos), isos_to_show(num_chem_isos), chem_id(num_chem_isos))
104 :
105 14 : call net_init(ierr)
106 14 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
107 :
108 14 : net_handle = alloc_net_handle(ierr)
109 14 : if (ierr /= 0) then
110 0 : write (*, *) 'alloc_net_handle failed'
111 0 : call mesa_error(__FILE__, __LINE__)
112 : end if
113 :
114 14 : call net_start_def(net_handle, ierr)
115 14 : if (ierr /= 0) then
116 0 : write (*, *) 'net_start_def failed'
117 0 : call mesa_error(__FILE__, __LINE__)
118 : end if
119 :
120 14 : call read_net_file(net_file, net_handle, ierr)
121 14 : if (ierr /= 0) then
122 0 : write (*, *) 'read_net_file failed ', trim(net_file)
123 0 : call mesa_error(__FILE__, __LINE__)
124 : end if
125 :
126 14 : call net_finish_def(net_handle, ierr)
127 14 : if (ierr /= 0) then
128 0 : write (*, *) 'net_finish_def failed'
129 0 : call mesa_error(__FILE__, __LINE__)
130 : end if
131 :
132 14 : allocate (reaction_id(rates_reaction_id_max), reaction_table(rates_reaction_id_max))
133 14 : allocate (rates_to_show(rates_reaction_id_max))
134 :
135 14 : cache_suffix = ''
136 14 : call net_setup_tables(net_handle, cache_suffix, info)
137 14 : if (ierr /= 0) then
138 0 : write (*, *) 'net_setup_tables failed'
139 0 : call mesa_error(__FILE__, __LINE__)
140 : end if
141 :
142 14 : call net_ptr(net_handle, g, ierr)
143 14 : if (ierr /= 0) then
144 0 : write (*, *) 'net_ptr failed'
145 0 : call mesa_error(__FILE__, __LINE__)
146 : end if
147 :
148 14 : species = g%num_isos
149 14 : num_reactions = g%num_reactions
150 :
151 14 : call get_chem_id_table(net_handle, species, chem_id, ierr)
152 14 : if (ierr /= 0) then
153 0 : write (*, *) 'get_chem_id_table failed'
154 0 : call mesa_error(__FILE__, __LINE__)
155 : end if
156 :
157 14 : call get_net_iso_table(net_handle, net_iso, ierr)
158 14 : if (ierr /= 0) then
159 0 : write (*, *) 'get_net_iso_table failed'
160 0 : call mesa_error(__FILE__, __LINE__)
161 : end if
162 :
163 14 : call get_reaction_id_table(net_handle, num_reactions, reaction_id, ierr)
164 14 : if (ierr /= 0) then
165 0 : write (*, *) 'get_reaction_id_table failed'
166 0 : call mesa_error(__FILE__, __LINE__)
167 : end if
168 :
169 14 : call get_net_reaction_table(net_handle, reaction_table, ierr)
170 14 : if (ierr /= 0) then
171 0 : write (*, *) 'get_net_reaction_table failed'
172 0 : call mesa_error(__FILE__, __LINE__)
173 : end if
174 :
175 14 : call do_test_net_alloc(species)
176 :
177 100 : end subroutine test_net_setup
178 :
179 14 : subroutine do_test_net_alloc(species)
180 : integer, intent(in) :: species
181 : allocate ( &
182 : xin(species), xin_copy(species), d_eps_nuc_dx(species), dxdt(species), &
183 14 : d_dxdt_dRho(species), d_dxdt_dT(species), d_dxdt_dx(species, species))
184 14 : end subroutine do_test_net_alloc
185 :
186 14 : subroutine Do_One_Net(symbolic)
187 : logical, intent(in) :: symbolic
188 14 : call do1_net(net_handle, symbolic)
189 14 : end subroutine Do_One_Net
190 :
191 4 : subroutine Setup_eos(handle)
192 : ! allocate and load the eos tables
193 : use eos_def
194 : use eos_lib
195 : integer, intent(out) :: handle
196 : integer :: ierr
197 : logical, parameter :: use_cache = .true.
198 4 : call eos_init('', use_cache, ierr)
199 4 : if (ierr /= 0) then
200 0 : write (*, *) 'eos_init failed in Setup_eos'
201 0 : call mesa_error(__FILE__, __LINE__)
202 : end if
203 4 : handle = alloc_eos_handle(ierr)
204 4 : if (ierr /= 0) then
205 0 : write (*, *) 'failed trying to allocate eos handle'
206 0 : call mesa_error(__FILE__, __LINE__)
207 : end if
208 4 : end subroutine Setup_eos
209 :
210 14 : subroutine test_net_cleanup
211 14 : call do_test_net_cleanup
212 14 : call free_net_handle(net_handle)
213 4 : end subroutine test_net_cleanup
214 :
215 14 : subroutine do_test_net_cleanup
216 14 : deallocate (xin)
217 14 : deallocate (d_eps_nuc_dx)
218 14 : deallocate (dxdt)
219 14 : deallocate (d_dxdt_dRho)
220 14 : deallocate (d_dxdt_dT)
221 14 : deallocate (d_dxdt_dx)
222 14 : end subroutine do_test_net_cleanup
223 :
224 12 : subroutine change_net(new_net_file)
225 : character(len=*), intent(in) :: new_net_file
226 12 : call test_net_cleanup
227 12 : call test_net_setup(new_net_file)
228 12 : end subroutine change_net
229 :
230 0 : subroutine Create_Plot_Files(species, xin)
231 : use utils_lib, only: mkdir
232 : integer, intent(in) :: species
233 : real(dp) :: xin(species)
234 0 : type(net_info) :: n
235 :
236 : integer:: i, j, k, ierr, num_reactions
237 0 : real(dp) :: T, logT, logT_min, logT_max
238 0 : real(dp) :: logRho_min, logRho_max, dlogT, dlogRho
239 : integer :: logT_points, logRho_points
240 : integer :: io, io_first, io_last, io_rho, io_tmp, io_params, num_out
241 : character(len=256) :: fname, dir
242 :
243 0 : real(dp), allocatable :: output_values(:, :, :)
244 : type(Net_General_Info), pointer :: g
245 :
246 : ! full range
247 0 : logT_max = 9.1d0
248 0 : logT_min = 6d0
249 0 : logRho_min = -3d0
250 0 : logRho_max = 10d0
251 :
252 : ! oxygen burning range
253 0 : logT_max = 9.5d0
254 0 : logT_min = 8d0
255 0 : logRho_min = -3d0
256 0 : logRho_max = 12d0
257 :
258 : ! test FL
259 0 : logT_max = 9d0
260 0 : logT_min = 7d0
261 0 : logRho_min = 5d0
262 0 : logRho_max = 10.2d0
263 :
264 : ! test FL
265 0 : logT_max = 8.4d0
266 0 : logT_min = 7.8d0
267 0 : logRho_min = 2d0
268 0 : logRho_max = 6d0
269 :
270 : ! test C+C
271 0 : logT_max = 10d0
272 0 : logT_min = 7.5d0
273 0 : logRho_min = 6d0
274 0 : logRho_max = 12d0
275 :
276 0 : logT_points = 251
277 0 : logRho_points = 251
278 :
279 0 : dir = 'plot_data'
280 0 : call mkdir(dir)
281 0 : write (*, *) trim(dir)
282 :
283 : 01 format(E30.22)
284 :
285 0 : dlogT = (logT_max - logT_min)/(logT_points - 1)
286 0 : dlogRho = (logRho_max - logRho_min)/(logRho_points - 1)
287 :
288 0 : io_params = 40
289 0 : io_rho = 41
290 0 : io_tmp = 42
291 0 : io_first = 43
292 :
293 0 : fname = trim(dir)//'/'//'params.data'
294 0 : open (unit=io_params, file=trim(fname))
295 0 : write (io_params, '(6f16.6)') &
296 0 : xin(net_iso(ih1)), xin(net_iso(ihe4)), xin(net_iso(ic12)), &
297 0 : xin(net_iso(in14)), xin(net_iso(io16))
298 0 : close (io_params)
299 :
300 0 : fname = trim(dir)//'/'//'rho.data'
301 0 : open (unit=io_rho, file=trim(fname))
302 :
303 0 : fname = trim(dir)//'/'//'tmp.data'
304 0 : open (unit=io_tmp, file=trim(fname))
305 :
306 0 : io = io_first - 1
307 0 : io_last = Open_Files(io, dir)
308 0 : num_out = io_last - io_first + 1
309 :
310 0 : call get_net_ptr(net_handle, g, ierr)
311 0 : if (ierr /= 0) return
312 :
313 0 : num_reactions = g%num_reactions
314 :
315 0 : allocate (output_values(logRho_points, logT_points, num_out))
316 :
317 : !xx$OMP PARALLEL DO PRIVATE(logT, T, j)
318 0 : do j = 1, logT_points
319 0 : logT = logT_min + dlogT*(j - 1)
320 0 : T = exp10(logT)
321 :
322 : call do_inner_loop(species, num_reactions, logT, T, j, output_values, n, xin, &
323 0 : logRho_points, logRho_min, dlogRho)
324 :
325 : end do
326 : !xx$OMP END PARALLEL DO
327 :
328 0 : write (*, *) 'write the files'
329 :
330 : ! write out the results
331 0 : do j = 1, logRho_points
332 0 : write (io_rho, 01) logRho_min + dlogRho*(j - 1)
333 : end do
334 0 : close (io_rho)
335 :
336 0 : do i = 1, logT_points
337 0 : write (io_tmp, 01) logT_min + dlogT*(i - 1)
338 : end do
339 0 : close (io_tmp)
340 :
341 0 : !$OMP PARALLEL DO PRIVATE(k)
342 : do k = 1, num_out
343 : write (*, *) k
344 : write (io_first + k - 1, '(e14.6)') output_values(:, :, k)
345 : end do
346 : !$OMP END PARALLEL DO
347 :
348 0 : do io = io_first, io_last
349 0 : close (io)
350 : end do
351 :
352 0 : end subroutine Create_Plot_Files
353 :
354 0 : subroutine do_inner_loop(species, num_reactions, logT, T, j, output_values, n, xin, &
355 : logRho_points, logRho_min, dlogRho)
356 : integer, intent(in) :: species, num_reactions
357 : real(dp), intent(in) :: logT, T
358 : integer, intent(in) :: j, logRho_points
359 : type(Net_info) :: n
360 : real(dp), intent(OUT) :: output_values(:, :, :)
361 : real(dp), intent(in) :: xin(species), logRho_min, dlogRho
362 :
363 : integer :: i
364 0 : real(dp) :: logRho, Rho
365 :
366 0 : do i = 1, logRho_points
367 0 : logRho = logRho_min + dlogRho*(i - 1)
368 0 : Rho = exp10(logRho)
369 : call do_one_net_eval(species, num_reactions, logT, T, logRho, Rho, &
370 0 : i, j, output_values, n, xin)
371 : end do
372 :
373 0 : end subroutine do_inner_loop
374 :
375 0 : subroutine do_one_net_eval(species, num_reactions, logT, T, logRho, Rho, &
376 0 : i, j, output_values, n, xin)
377 : use chem_lib, only: composition_info
378 : integer, intent(in) :: species, num_reactions
379 : real(dp), intent(in) :: logT, T, logRho, Rho
380 : integer, intent(in) :: i, j
381 : type(Net_General_Info), pointer :: g => null()
382 : type(net_info) :: n
383 : real(dp), intent(OUT) :: output_values(:, :, :)
384 : real(dp), intent(in) :: xin(species)
385 :
386 0 : real(dp) :: z, abar, zbar, z2bar, z53bar, ye, sum, mx, weak_rate_factor
387 0 : real(dp), target :: rate_factors_a(num_reactions)
388 0 : real(dp), pointer :: rate_factors(:)
389 :
390 0 : real(dp) :: eps_nuc
391 0 : real(dp) :: d_eps_nuc_dT
392 0 : real(dp) :: d_eps_nuc_dRho
393 0 : real(dp) :: d_eps_nuc_dx(species)
394 : ! partial derivatives wrt mass fractions
395 :
396 0 : real(dp) :: dxdt(species)
397 : ! rate of change of mass fractions caused by nuclear reactions
398 0 : real(dp) :: d_dxdt_dRho(species)
399 0 : real(dp) :: d_dxdt_dT(species)
400 0 : real(dp) :: d_dxdt_dx(species, species)
401 0 : real(dp) :: eps_nuc_categories(num_categories)
402 :
403 0 : integer :: info, k, h1, he4, chem_id(species)
404 0 : real(dp) :: xh, xhe, mass_correction
405 0 : real(dp), dimension(species) :: dabar_dx, dzbar_dx, dmc_dx
406 : logical :: skip_jacobian
407 :
408 0 : rate_factors => rate_factors_a
409 :
410 0 : h1 = net_iso(ih1)
411 0 : he4 = net_iso(ihe4)
412 :
413 0 : g => n%g
414 :
415 0 : call get_chem_id_table(net_handle, species, chem_id, info)
416 0 : if (info /= 0) call mesa_error(__FILE__, __LINE__)
417 :
418 : call composition_info( &
419 0 : species, chem_id, xin, xh, xhe, z, abar, zbar, z2bar, z53bar, ye, &
420 0 : mass_correction, sum, dabar_dx, dzbar_dx, dmc_dx)
421 :
422 0 : rate_factors(:) = 1
423 0 : weak_rate_factor = 1
424 0 : skip_jacobian = .false.
425 :
426 : call net_get(net_handle, skip_jacobian, n, species, num_reactions, &
427 0 : xin, T, logT, Rho, logRho, &
428 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
429 : rate_factors, weak_rate_factor, &
430 : std_reaction_Qs, std_reaction_neuQs, &
431 0 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
432 0 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
433 : screening_mode, &
434 : eps_nuc_categories, eps_neu_total, &
435 0 : info)
436 0 : if (info /= 0) then
437 0 : write (*, *) 'bad result from net_get'
438 0 : call mesa_error(__FILE__, __LINE__)
439 : end if
440 :
441 0 : sum = 0d0
442 0 : mx = 0d0
443 0 : do k = species, 1, -1
444 0 : if (abs(dxdt(k)) > mx) mx = abs(dxdt(k))
445 0 : sum = sum + dxdt(k)
446 : end do
447 0 : if (mx <= 0) mx = 1d0
448 0 : if (sum - sum /= 0) then
449 0 : write (*, *) logRho, logT
450 0 : call mesa_error(__FILE__, __LINE__)
451 : end if
452 :
453 0 : k = 1; output_values(i, j, k) = safe_log10(eps_nuc)
454 :
455 0 : k = k + 1; output_values(i, j, k) = sum/max(1d-20, mx)
456 :
457 0 : k = k + 1; output_values(i, j, k) = d_eps_nuc_dT*T/max(1d-20, eps_nuc)
458 0 : k = k + 1; output_values(i, j, k) = d_eps_nuc_dRho*Rho/max(1d-20, eps_nuc)
459 0 : k = k + 1; output_values(i, j, k) = safe_log10(abs(d_eps_nuc_dx(h1)))
460 0 : k = k + 1; output_values(i, j, k) = safe_log10(abs(d_eps_nuc_dx(he4)))
461 :
462 0 : if (k > max_files) then
463 0 : write (*, *) 'need to enlarge max_files'
464 0 : call mesa_error(__FILE__, __LINE__)
465 : end if
466 :
467 0 : end subroutine do_one_net_eval
468 :
469 0 : integer function Open_Files(io_start, dir)
470 : integer, intent(in) :: io_start
471 : character(len=256), intent(in) :: dir
472 : character(len=256) :: fname
473 : integer :: io
474 0 : io = io_start
475 :
476 0 : fname = trim(dir)//'/'//'log_net_eps.data'
477 0 : io = io + 1; open (unit=io, file=trim(fname))
478 :
479 0 : fname = trim(dir)//'/'//'sum_dxdt_div_max_dxdt.data'
480 0 : io = io + 1; open (unit=io, file=trim(fname))
481 :
482 0 : fname = trim(dir)//'/'//'d_lneps_dlnT.data'
483 0 : io = io + 1; open (unit=io, file=trim(fname))
484 :
485 0 : fname = trim(dir)//'/'//'d_lneps_dlnRho.data'
486 0 : io = io + 1; open (unit=io, file=trim(fname))
487 :
488 0 : fname = trim(dir)//'/'//'d_eps_nuc_dxh1.data'
489 0 : io = io + 1; open (unit=io, file=trim(fname))
490 :
491 0 : fname = trim(dir)//'/'//'d_eps_nuc_dxhe4.data'
492 0 : io = io + 1; open (unit=io, file=trim(fname))
493 :
494 0 : Open_Files = io
495 :
496 0 : end function Open_Files
497 :
498 14 : subroutine set_composition(species, xin)
499 : integer, intent(in) :: species
500 : real(dp), intent(OUT) :: xin(species)
501 :
502 14 : real(dp) :: sum
503 : integer :: i, adjustment_iso
504 :
505 14 : eta = 0
506 :
507 14 : adjustment_iso = net_iso(img24)
508 :
509 14 : if (net_file == 'basic.net') then
510 :
511 2 : adjustment_iso = net_iso(img24)
512 :
513 18 : xin = 0
514 2 : xin(net_iso(ih1)) = 0.655186E+00 ! h1
515 2 : xin(net_iso(ihe4)) = 0.31002164D+00 ! he4
516 2 : xin(net_iso(ic12)) = 0.002725D-01 ! c12
517 2 : xin(net_iso(in14)) = 0.203101D-01 + 0.612124D-06 + 0.109305D-02 + 0.356004D-04
518 2 : xin(net_iso(io16)) = 0.094000D-01 ! o16
519 2 : xin(net_iso(ine20)) = 0.162163D-02 ! ne20
520 2 : xin(net_iso(img24)) = 0.658226D-25 ! mg24
521 2 : xin(net_iso(ihe3)) = 0.201852D-02 ! he3
522 :
523 12 : else if (net_file == 'o18_and_ne22.net') then
524 :
525 2 : adjustment_iso = net_iso(img24)
526 :
527 22 : xin = 0
528 2 : xin(net_iso(ih1)) = 0.655186E+00
529 2 : xin(net_iso(ihe4)) = 0.31002164D+00
530 2 : xin(net_iso(ic12)) = 0.002725D-01
531 2 : xin(net_iso(in14)) = 0.203101D-01 + 0.612124D-06 + 0.109305D-02
532 2 : xin(net_iso(io16)) = 0.094000D-01
533 2 : xin(net_iso(io18)) = 1d-20
534 2 : xin(net_iso(ine20)) = 0.162163D-02
535 2 : xin(net_iso(img24)) = 0.658226D-25
536 2 : xin(net_iso(ine22)) = 0.201852D-02
537 :
538 10 : else if (net_file == 'pp_extras.net') then
539 :
540 2 : adjustment_iso = net_iso(img24)
541 :
542 26 : xin = 0
543 2 : xin(net_iso(ih1)) = 0.655186E+00 ! h1
544 2 : xin(net_iso(ihe4)) = 0.31002164D+00 ! he4
545 2 : xin(net_iso(ic12)) = 0.002725D-01 ! c12
546 2 : xin(net_iso(in14)) = 0.203101D-01 + 0.612124D-06 + 0.109305D-02 + 0.356004D-04
547 2 : xin(net_iso(io16)) = 0.094000D-01 ! o16
548 2 : xin(net_iso(ine20)) = 0.162163D-02 ! ne20
549 2 : xin(net_iso(img24)) = 0.658226D-25 ! mg24
550 :
551 2 : xin(net_iso(ih2)) = 0.632956D-17 ! h2
552 2 : xin(net_iso(ihe3)) = 0.201852D-02 ! he3
553 2 : xin(net_iso(ili7)) = 0.664160D-15 ! li7
554 2 : xin(net_iso(ibe7)) = 0.103866D-15 ! be7
555 :
556 8 : else if (net_file == 'cno_extras.net') then
557 :
558 2 : adjustment_iso = net_iso(img24)
559 44 : xin = 0
560 2 : xin(net_iso(ih1)) = 0.173891680788D-01 ! h1
561 2 : xin(net_iso(ihe4)) = 0.963245225401D+00 ! he4
562 2 : xin(net_iso(ic12)) = 0.238935745993D-03 ! c12
563 2 : xin(net_iso(in14)) = 0.134050688300D-01 ! n14
564 2 : xin(net_iso(io16)) = 0.268791618452D-03 ! o16
565 2 : xin(net_iso(ine20)) = 0.180001692845D-02 ! ne20
566 2 : xin(net_iso(img24)) = 0.353667702698D-02 ! mg24
567 :
568 2 : xin(net_iso(ic13)) = 0.717642727071D-04 ! c13
569 2 : xin(net_iso(in13)) = 0.370732258156D-09 ! n13
570 2 : xin(net_iso(in15)) = 0.450484708137D-06 ! n15
571 2 : xin(net_iso(io14)) = 0.100000000000D-49 ! o14
572 2 : xin(net_iso(io15)) = 0.874815374966D-10 ! o15
573 2 : xin(net_iso(if17)) = 0.100000000000D-49 ! f17
574 2 : xin(net_iso(if18)) = 0.100000000000D-49 ! f18
575 2 : xin(net_iso(ine18)) = 0.100000000000D-49 ! ne18
576 2 : xin(net_iso(ine19)) = 0.100000000000D-49 ! ne19
577 2 : xin(net_iso(img22)) = 0.439011547696D-04 ! mg22
578 :
579 6 : else if (net_file == 'pp_cno_extras_o18_ne22.net') then
580 :
581 0 : adjustment_iso = net_iso(img24)
582 0 : xin = 0
583 0 : xin(net_iso(ih1)) = 0.173891680788D-01 ! h1
584 0 : xin(net_iso(ihe4)) = 0.963245225401D+00 ! he4
585 0 : xin(net_iso(ic12)) = 0.238935745993D-03 ! c12
586 0 : xin(net_iso(in14)) = 0.134050688300D-01 ! n14
587 0 : xin(net_iso(io16)) = 0.268791618452D-03 ! o16
588 0 : xin(net_iso(ine20)) = 0.180001692845D-02 ! ne20
589 0 : xin(net_iso(img24)) = 0.353667702698D-02 ! mg24
590 :
591 0 : xin(net_iso(ic13)) = 0.717642727071D-04 ! c13
592 0 : xin(net_iso(in13)) = 0.370732258156D-09 ! n13
593 0 : xin(net_iso(in15)) = 0.450484708137D-06 ! n15
594 0 : xin(net_iso(io14)) = 0.100000000000D-49 ! o14
595 0 : xin(net_iso(io15)) = 0.874815374966D-10 ! o15
596 0 : xin(net_iso(if17)) = 0.100000000000D-49 ! f17
597 0 : xin(net_iso(if18)) = 0.100000000000D-49 ! f18
598 0 : xin(net_iso(ine18)) = 0.100000000000D-49 ! ne18
599 0 : xin(net_iso(ine19)) = 0.100000000000D-49 ! ne19
600 0 : xin(net_iso(img22)) = 0.439011547696D-04 ! mg22
601 :
602 6 : else if (net_file == 'approx21_cr60_plus_co56.net') then
603 :
604 2 : adjustment_iso = net_iso(img24)
605 46 : xin = 0
606 2 : xin(net_iso(ihe4)) = 3.4555392534813939D-01
607 2 : xin(net_iso(io16)) = 1.9367778420430937D-01
608 2 : xin(net_iso(ih1)) = 1.9052931367172501D-01
609 2 : xin(net_iso(isi28)) = 9.1023936032064601D-02
610 2 : xin(net_iso(ini56)) = 5.8492459653295005D-02
611 2 : xin(net_iso(is32)) = 4.1416642476999908D-02
612 2 : xin(net_iso(ine20)) = 2.0927782332477558D-02
613 2 : xin(net_iso(ic12)) = 2.0589617104448312D-02
614 2 : xin(net_iso(img24)) = 1.8975914108406104D-02
615 2 : xin(net_iso(iar36)) = 7.0600338785939956D-03
616 2 : xin(net_iso(ica40)) = 4.9182171981353726D-03
617 2 : xin(net_iso(ife52)) = 3.1618820806005280D-03
618 2 : xin(net_iso(in14)) = 1.9878460082760952D-03
619 2 : xin(net_iso(ife56)) = 7.1668776621130190D-04
620 2 : xin(net_iso(ico56)) = 4.4974971099921181D-04
621 2 : xin(net_iso(ife54)) = 3.3432423890988422D-04
622 2 : xin(net_iso(icr48)) = 1.7602626907769762D-04
623 2 : xin(net_iso(ihe3)) = 5.1650097191631545D-06
624 2 : xin(net_iso(iti44)) = 2.6452422203389906D-06
625 2 : xin(net_iso(icr60)) = 4.7653353095937982D-08
626 2 : xin(net_iso(iprot)) = 1.2739022407246250D-11
627 2 : xin(net_iso(ineut)) = 0.0000000000000000D+00
628 :
629 4 : else if (net_file == 'approx21_cr60_plus_fe53_fe55_co56.net') then
630 :
631 0 : adjustment_iso = net_iso(img24)
632 0 : xin = 0
633 0 : xin(net_iso(ihe4)) = 3.4555392534813939D-01
634 0 : xin(net_iso(io16)) = 1.9367778420430937D-01
635 0 : xin(net_iso(ih1)) = 1.9052931367172501D-01
636 0 : xin(net_iso(isi28)) = 9.1023936032064601D-02
637 0 : xin(net_iso(ini56)) = 5.8492459653295005D-02
638 0 : xin(net_iso(is32)) = 4.1416642476999908D-02
639 0 : xin(net_iso(ine20)) = 2.0927782332477558D-02
640 0 : xin(net_iso(ic12)) = 2.0589617104448312D-02
641 0 : xin(net_iso(img24)) = 1.8975914108406104D-02
642 0 : xin(net_iso(iar36)) = 7.0600338785939956D-03
643 0 : xin(net_iso(ica40)) = 4.9182171981353726D-03
644 0 : xin(net_iso(ife52)) = 3.1618820806005280D-03
645 0 : xin(net_iso(in14)) = 1.9878460082760952D-03
646 0 : xin(net_iso(ife56)) = 7.1668776621130190D-04
647 0 : xin(net_iso(ico56)) = 4.4974971099921181D-04
648 0 : xin(net_iso(ife54)) = 3.3432423890988422D-04
649 0 : xin(net_iso(icr48)) = 1.7602626907769762D-04
650 0 : xin(net_iso(ihe3)) = 5.1650097191631545D-06
651 0 : xin(net_iso(iti44)) = 2.6452422203389906D-06
652 0 : xin(net_iso(icr60)) = 4.7653353095937982D-08
653 0 : xin(net_iso(iprot)) = 1.2739022407246250D-11
654 0 : xin(net_iso(ife53)) = 0.0000000000000000D+00
655 0 : xin(net_iso(ife55)) = 0.0000000000000000D+00
656 0 : xin(net_iso(ineut)) = 0.0000000000000000D+00
657 :
658 : else if (net_file == 'approx21_new.net' .or. &
659 : net_file == 'approx21_old.net' .or. &
660 4 : net_file == 'approx21_plus_co56.net' .or. &
661 : net_file == 'approx21.net') then
662 :
663 4 : adjustment_iso = net_iso(img24)
664 90 : xin = 0
665 4 : xin(net_iso(ife56)) = 8.0387021484318166D-01
666 4 : xin(net_iso(ife54)) = 1.6096648736760832D-01
667 4 : xin(net_iso(icr56)) = 2.9480945535920525D-02
668 4 : xin(net_iso(ihe4)) = 4.8624637161320565D-03
669 4 : xin(net_iso(ini56)) = 2.8376270731890360D-04
670 4 : xin(net_iso(isi28)) = 1.5018628906135739D-04
671 4 : xin(net_iso(is32)) = 1.1613271635573457D-04
672 4 : xin(net_iso(iprot)) = 1.1139431633673653D-04
673 4 : xin(net_iso(ica40)) = 5.3688377473494185D-05
674 4 : xin(net_iso(iar36)) = 5.2702831822567062D-05
675 4 : xin(net_iso(ife52)) = 3.7866504131935185D-05
676 4 : xin(net_iso(icr48)) = 5.8401123037667974D-06
677 4 : xin(net_iso(ineut)) = 4.9141227703118397D-06
678 4 : xin(net_iso(iti44)) = 1.5085038561154746D-06
679 4 : xin(net_iso(io16)) = 9.5384019609049255D-07
680 4 : xin(net_iso(img24)) = 6.3808207717725580D-07
681 4 : xin(net_iso(ic12)) = 2.9048656673991868D-07
682 4 : xin(net_iso(ine20)) = 9.6468865023609427D-09
683 4 : xin(net_iso(ihe3)) = 4.6203603862263096D-80
684 4 : xin(net_iso(in14)) = 7.5867472225841235D-99
685 :
686 : else
687 :
688 0 : write (*, *) 'net_file '//trim(net_file)
689 0 : call mesa_error(__FILE__, __LINE__, 'set_composition: do not recognize net_file')
690 :
691 : end if
692 :
693 14 : sum = 0d0
694 246 : do i = 1, species
695 232 : if (xin(i) < 1d-99) xin(i) = 1d-99
696 246 : sum = sum + xin(i)
697 : end do
698 14 : if (abs(1d0 - sum) > 1d-4) write (*, *) 'change abundance sum by', 1d0 - sum
699 14 : xin(adjustment_iso) = xin(adjustment_iso) + (1d0 - sum)
700 14 : if (xin(adjustment_iso) < 0d0) call mesa_error(__FILE__, __LINE__, 'error in sum of abundances')
701 :
702 14 : end subroutine set_composition
703 :
704 0 : subroutine read_test_data(filename, n, rho_vec, T_vec, ierr)
705 : ! the data files have columns of mass, radius, density, temp
706 : use utils_lib
707 : character(len=*), intent(in) :: filename
708 : integer, intent(out) :: n
709 : real(dp), dimension(:), pointer :: rho_vec, T_vec ! to be allocated and filled
710 : integer, intent(out) :: ierr
711 :
712 : integer :: iounit, i
713 0 : real(dp) :: junk
714 :
715 0 : ierr = 0
716 0 : open (newunit=iounit, file=trim(filename), action='read', iostat=ierr)
717 0 : if (ierr /= 0) then
718 0 : write (*, *) 'failed to open ', trim(filename)
719 0 : return
720 : end if
721 :
722 0 : i = 0
723 0 : do
724 0 : read (unit=iounit, fmt=*, iostat=ierr) junk, junk, junk, junk
725 0 : if (ierr == 0) then
726 0 : i = i + 1; cycle
727 : end if
728 0 : ierr = 0
729 0 : n = i
730 0 : exit
731 : end do
732 0 : rewind (iounit)
733 :
734 0 : allocate (rho_vec(n), T_vec(n), stat=ierr); if (ierr /= 0) return
735 :
736 0 : do i = 1, n
737 0 : read (iounit, *) junk, junk, rho_vec(i), T_vec(i)
738 : end do
739 :
740 0 : close (iounit)
741 :
742 0 : end subroutine read_test_data
743 :
744 0 : subroutine Do_One_Test(net_file, do_timing)
745 : use chem_lib, only: composition_info
746 : use rates_lib
747 : character(len=*), intent(in) :: net_file
748 : logical, intent(in) :: do_timing
749 0 : call Do_One_Testcase(net_file, do_timing, .false.)
750 0 : end subroutine Do_One_Test
751 :
752 0 : subroutine Do_One_Test_and_show_Qs(net_file, do_timing)
753 0 : use chem_lib, only: composition_info
754 : use rates_lib
755 : character(len=*), intent(in) :: net_file
756 : logical, intent(in) :: do_timing
757 0 : call Do_One_Testcase(net_file, do_timing, .true.)
758 0 : end subroutine Do_One_Test_and_show_Qs
759 :
760 0 : subroutine Do_One_Testcase(net_file, do_timing, show_Qs)
761 0 : use chem_lib, only: composition_info
762 : use rates_lib
763 : character(len=*), intent(in) :: net_file
764 : logical, intent(in) :: do_timing, show_Qs
765 :
766 0 : real(dp) :: logRho, logT, Rho, T, xsum, Q1, &
767 0 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, weak_rate_factor, &
768 0 : dvardx, dvardx_0, dx_0, err, var_0, xdum, &
769 0 : eps_nuc_categories(num_categories), xh, xhe, mass_correction !approx_abar, approx_zbar
770 : integer :: i, j, info, ierr
771 : integer :: j_dx, j_dx_sink
772 0 : real(dp), dimension(:), pointer :: d_eps_nuc_dx, dabar_dx, dzbar_dx, dmc_dx
773 0 : real(dp), pointer :: rate_factors(:), &
774 0 : actual_Qs(:), actual_neuQs(:)
775 0 : logical, pointer :: from_weaklib(:)
776 : logical :: skip_jacobian, doing_d_dlnd, doing_dx
777 : real(dp), dimension(:), pointer :: &
778 0 : rate_raw, rate_raw_dT, &
779 0 : rate_screened
780 0 : type(net_info) :: n
781 :
782 : include 'formats'
783 :
784 0 : write (*, *) 'Do_One_Test '//trim(net_file)
785 :
786 0 : if (do_timing) call mesa_error(__FILE__, __LINE__, 'no support for do_timing')
787 :
788 0 : call test_net_setup(net_file)
789 :
790 0 : ierr = 0
791 :
792 0 : write (*, *) 'species', species
793 :
794 0 : allocate (d_eps_nuc_dx(species), dabar_dx(species), dzbar_dx(species), dmc_dx(species))
795 :
796 0 : info = 0
797 :
798 : allocate ( &
799 : rate_factors(num_reactions), &
800 : actual_Qs(num_reactions), &
801 : actual_neuQs(num_reactions), &
802 : from_weaklib(num_reactions), &
803 0 : stat=info)
804 0 : if (info /= 0) call mesa_error(__FILE__, __LINE__)
805 :
806 0 : rate_factors(:) = 1
807 0 : weak_rate_factor = 1
808 :
809 : if (.false.) then ! get neutrino Q
810 :
811 : Q1 = eval_neutrino_Q(img22, is30)
812 : write (*, 1) 'Qneu mg22->s30', Q1
813 :
814 : Q1 = eval_neutrino_Q(is30, ini56)
815 : write (*, 1) 'Qneu s30->ni56', Q1
816 :
817 : Q1 = eval_neutrino_Q(ica38, ini56)
818 : write (*, 1) 'Qneu ca38->ni56', Q1
819 :
820 : Q1 = eval_neutrino_Q(ini56, ige64)
821 : write (*, 1) 'Qneu ni56->ge64', Q1
822 :
823 : Q1 = eval_neutrino_Q(ige64, ise68)
824 : write (*, 1) 'Qneu ge64->se68', Q1
825 :
826 : Q1 = eval_neutrino_Q(ise68, ikr72)
827 : write (*, 1) 'Qneu se68->kr72', Q1
828 :
829 : Q1 = eval_neutrino_Q(ikr72, isr76)
830 : write (*, 1) 'Qneu kr72->sr76', Q1
831 :
832 : Q1 = eval_neutrino_Q(isr76, isn104)
833 : write (*, 1) 'Qneu sr76->sn104', Q1
834 :
835 : stop
836 : end if
837 :
838 : if (.false.) then ! get reaction Q
839 : ! co55 -> fe55
840 : Q1 = isoB(ife55) - isoB(ico55)
841 : write (*, 1) 'Q co55->fe55', Q1
842 : stop
843 : end if
844 :
845 0 : xin = 0
846 0 : eta = 0
847 :
848 0 : if (net_file == 'mesa_201.net') then
849 :
850 0 : nrates_to_show = 2
851 :
852 0 : rates_to_show(1:nrates_to_show) = [ &
853 : ir_ar36_ag_ca40, &
854 0 : ir_ca40_ga_ar36]
855 :
856 0 : xin(net_iso(ife56)) = 6.3551174304779179D-01
857 0 : xin(net_iso(icr52)) = 1.0518849309100423D-01
858 0 : xin(net_iso(ini60)) = 6.8579809304547934D-02
859 0 : xin(net_iso(ife55)) = 3.0800958229563902D-02
860 0 : xin(net_iso(imn55)) = 2.1373990699858084D-02
861 0 : xin(net_iso(ico57)) = 2.0785551124804885D-02
862 0 : xin(net_iso(ife57)) = 1.6083543077536507D-02
863 0 : xin(net_iso(imn53)) = 1.5959192021165327D-02
864 0 : xin(net_iso(ico59)) = 1.3346191916961030D-02
865 0 : xin(net_iso(ife58)) = 1.2810477227656202D-02
866 0 : xin(net_iso(ife54)) = 1.2428871023625330D-02
867 0 : xin(net_iso(ini62)) = 1.1831182533246063D-02
868 0 : xin(net_iso(icr53)) = 8.1041573656346795D-03
869 0 : xin(net_iso(ini61)) = 5.1462544919734727D-03
870 0 : xin(net_iso(imn54)) = 4.8562102089927551D-03
871 0 : xin(net_iso(icr54)) = 3.4271335428296490D-03
872 0 : xin(net_iso(ini59)) = 2.9811021846265560D-03
873 0 : xin(net_iso(ini58)) = 2.8713349650366189D-03
874 0 : xin(net_iso(iv51)) = 2.0988150935152424D-03
875 0 : xin(net_iso(ico58)) = 2.0282210582857861D-03
876 0 : xin(net_iso(icr51)) = 1.2727926750247761D-03
877 0 : xin(net_iso(icr50)) = 3.5421790633561727D-04
878 0 : xin(net_iso(icu63)) = 3.0335040211002022D-04
879 0 : xin(net_iso(iti48)) = 2.8512639104984498D-04
880 0 : xin(net_iso(iti50)) = 2.8394752519368671D-04
881 0 : xin(net_iso(ico56)) = 2.4310701594699283D-04
882 0 : xin(net_iso(ico60)) = 1.5319574812323961D-04
883 0 : xin(net_iso(ihe4)) = 1.3780017854631601D-04
884 0 : xin(net_iso(imn56)) = 1.1066944504563417D-04
885 0 : xin(net_iso(iv50)) = 8.7723703257792468D-05
886 0 : xin(net_iso(iv49)) = 8.0539507322740820D-05
887 0 : xin(net_iso(icu61)) = 6.0898214714637941D-05
888 0 : xin(net_iso(iti49)) = 5.8026224632063015D-05
889 0 : xin(net_iso(imn52)) = 4.2516037186521133D-05
890 0 : xin(net_iso(ico61)) = 4.1332712278653746D-05
891 0 : xin(net_iso(ini63)) = 3.9285182071536010D-05
892 0 : xin(net_iso(ico55)) = 3.3984664667118780D-05
893 0 : xin(net_iso(ife59)) = 3.1874001320244153D-05
894 0 : xin(net_iso(izn64)) = 2.7904240321969555D-05
895 0 : xin(net_iso(izn66)) = 2.3416686728782276D-05
896 0 : xin(net_iso(icu62)) = 2.1144678609352833D-05
897 0 : xin(net_iso(ini57)) = 1.1679099363693715D-05
898 0 : xin(net_iso(iv52)) = 1.0016381776869645D-05
899 0 : xin(net_iso(ini64)) = 8.9564547480773243D-06
900 0 : xin(net_iso(iti47)) = 8.7562205406651916D-06
901 0 : xin(net_iso(icu65)) = 8.3830901771893844D-06
902 0 : xin(net_iso(iti46)) = 6.5019528109988749D-06
903 0 : xin(net_iso(icu64)) = 6.1556773376380492D-06
904 0 : xin(net_iso(iar38)) = 5.0886468271422554D-06
905 0 : xin(net_iso(ife53)) = 4.4893983930970075D-06
906 0 : xin(net_iso(is34)) = 4.2369862519232918D-06
907 0 : xin(net_iso(izn65)) = 4.1809380230467465D-06
908 0 : xin(net_iso(icr55)) = 3.5647078269917385D-06
909 0 : xin(net_iso(imn51)) = 1.3849206094166064D-06
910 0 : xin(net_iso(ife60)) = 1.0919469947619407D-06
911 0 : xin(net_iso(ih1)) = 9.8211241034612710D-07
912 0 : xin(net_iso(ica44)) = 6.4721481111312556D-07
913 0 : xin(net_iso(isi30)) = 6.0482126630410841D-07
914 0 : xin(net_iso(ica42)) = 5.8926684360031471D-07
915 0 : xin(net_iso(isi28)) = 5.5342250413192408D-07
916 0 : xin(net_iso(iv48)) = 5.4911502105605942D-07
917 0 : xin(net_iso(iv53)) = 5.4571869339657211D-07
918 0 : xin(net_iso(icu60)) = 4.5761289432308086D-07
919 0 : xin(net_iso(izn63)) = 4.4652045082610858D-07
920 0 : xin(net_iso(isc47)) = 4.4099241256913694D-07
921 0 : xin(net_iso(ini56)) = 3.9547421607331126D-07
922 0 : xin(net_iso(iti51)) = 3.6358450091693021D-07
923 0 : xin(net_iso(izn62)) = 3.0268870369662295D-07
924 0 : xin(net_iso(is32)) = 3.0184105591459131D-07
925 0 : xin(net_iso(icr49)) = 2.7706803242086906D-07
926 0 : xin(net_iso(isc45)) = 2.5466905171154329D-07
927 0 : xin(net_iso(ik39)) = 1.9906885458096793D-07
928 0 : xin(net_iso(icl35)) = 1.5826438708200070D-07
929 0 : xin(net_iso(is33)) = 1.2794750349676913D-07
930 0 : xin(net_iso(ip31)) = 1.1805533268957108D-07
931 0 : xin(net_iso(icl37)) = 1.0108890802543261D-07
932 0 : xin(net_iso(ica43)) = 9.5766190955127930D-08
933 0 : xin(net_iso(iar36)) = 8.3287350202069049D-08
934 0 : xin(net_iso(isi29)) = 7.6157023642649312D-08
935 0 : xin(net_iso(isc49)) = 7.2645544618960897D-08
936 0 : xin(net_iso(icu59)) = 6.9036259006119691D-08
937 0 : xin(net_iso(ica40)) = 6.4387988887989985D-08
938 0 : xin(net_iso(isc46)) = 6.1294200120692868D-08
939 0 : xin(net_iso(iar37)) = 5.0049708949178339D-08
940 0 : xin(net_iso(ineut)) = 4.7236564548991118D-08
941 0 : xin(net_iso(isc48)) = 3.7392395931746378D-08
942 0 : xin(net_iso(ife52)) = 3.0982331086000098D-08
943 0 : xin(net_iso(icr56)) = 2.9477366114308358D-08
944 0 : xin(net_iso(ica41)) = 2.7100443580454983D-08
945 0 : xin(net_iso(is35)) = 2.6527306613430685D-08
946 0 : xin(net_iso(ica45)) = 2.5532508173462626D-08
947 0 : xin(net_iso(ico62)) = 2.3608043613616947D-08
948 0 : xin(net_iso(iar39)) = 2.3503192609361190D-08
949 0 : xin(net_iso(ica46)) = 2.2455595751043146D-08
950 0 : xin(net_iso(iv47)) = 1.9789662501752072D-08
951 0 : xin(net_iso(icu66)) = 1.9285467280104737D-08
952 0 : xin(net_iso(icl36)) = 1.8710358753573163D-08
953 0 : xin(net_iso(is36)) = 1.6433794658525574D-08
954 0 : xin(net_iso(ik41)) = 1.1021812817088955D-08
955 0 : xin(net_iso(ini65)) = 1.0432196346423500D-08
956 0 : xin(net_iso(ip33)) = 8.9625871046594568D-09
957 0 : xin(net_iso(ik40)) = 7.2106087354006743D-09
958 0 : xin(net_iso(iar40)) = 7.1174073521523365D-09
959 0 : xin(net_iso(iti45)) = 5.6870435962784455D-09
960 0 : xin(net_iso(ip32)) = 4.9537776441010461D-09
961 0 : xin(net_iso(icr48)) = 3.0709436939798925D-09
962 0 : xin(net_iso(isc44)) = 2.7862949914482315D-09
963 0 : xin(net_iso(isc43)) = 1.4283233379580965D-09
964 0 : xin(net_iso(isi31)) = 1.3747008133379580D-09
965 0 : xin(net_iso(iti52)) = 1.2866447687101849D-09
966 0 : xin(net_iso(ico63)) = 1.0456195723572430D-09
967 0 : xin(net_iso(iti44)) = 7.1579364938842059D-10
968 0 : xin(net_iso(io16)) = 5.6818053126812287D-10
969 0 : xin(net_iso(ial27)) = 5.6117840295305725D-10
970 0 : xin(net_iso(ica47)) = 5.4798226854137171D-10
971 0 : xin(net_iso(img24)) = 3.8679049700264050D-10
972 0 : xin(net_iso(ini66)) = 2.8640604416758339D-10
973 0 : xin(net_iso(izn61)) = 2.4995195933116682D-10
974 0 : xin(net_iso(ica48)) = 1.9626419275146166D-10
975 0 : xin(net_iso(ife61)) = 1.8711288748163173D-10
976 0 : xin(net_iso(ik42)) = 1.4923884785773920D-10
977 0 : xin(net_iso(isi32)) = 1.4443919428884894D-10
978 0 : xin(net_iso(ip30)) = 1.3908837795296053D-10
979 0 : xin(net_iso(ic12)) = 1.3185610207511085D-10
980 0 : xin(net_iso(ik43)) = 1.1518951971920992D-10
981 0 : xin(net_iso(iv54)) = 8.8440974843643898D-11
982 0 : xin(net_iso(img26)) = 8.5260453209638504D-11
983 0 : xin(net_iso(icl38)) = 2.5853881572793256D-11
984 0 : xin(net_iso(isc50)) = 1.7292118708137885D-11
985 0 : xin(net_iso(iar41)) = 1.0665700916421081D-11
986 0 : xin(net_iso(img25)) = 9.2051027558630546D-12
987 0 : xin(net_iso(ial28)) = 8.9602126799277990D-12
988 0 : xin(net_iso(izn60)) = 7.0658622923470296D-12
989 0 : xin(net_iso(ip34)) = 4.2003981867966896D-12
990 0 : xin(net_iso(icr57)) = 2.8195149736768269D-12
991 0 : xin(net_iso(ine20)) = 1.1829006560529443D-12
992 0 : xin(net_iso(ife62)) = 1.0149415562457171D-12
993 0 : xin(net_iso(ik44)) = 5.5654687896999145D-13
994 0 : xin(net_iso(is31)) = 4.0845373817884839D-13
995 0 : xin(net_iso(iv55)) = 2.8610531204397973D-13
996 0 : xin(net_iso(ina23)) = 2.5122330767236196D-13
997 0 : xin(net_iso(is37)) = 2.1529809224524208D-13
998 0 : xin(net_iso(iti53)) = 1.4129020276964862D-13
999 0 : xin(net_iso(iar35)) = 1.2934962278023034D-13
1000 0 : xin(net_iso(ial26)) = 1.2407881461251650D-13
1001 0 : xin(net_iso(in15)) = 1.1945765698183132D-13
1002 0 : xin(net_iso(ico64)) = 8.1590671587313667D-14
1003 0 : xin(net_iso(img27)) = 6.7969591881380018D-14
1004 0 : xin(net_iso(ih2)) = 5.6613884110319004D-14
1005 0 : xin(net_iso(ini67)) = 5.1506647175988674D-14
1006 0 : xin(net_iso(ica39)) = 3.8945192023978035D-14
1007 0 : xin(net_iso(ini55)) = 3.0883773851523300D-14
1008 0 : xin(net_iso(ine22)) = 1.3808913342478939D-14
1009 0 : xin(net_iso(ica49)) = 1.0542673783683999D-14
1010 0 : xin(net_iso(isc51)) = 9.2083686560605166D-15
1011 0 : xin(net_iso(isi27)) = 8.3464191225256989D-15
1012 0 : xin(net_iso(ine21)) = 6.6697557229646203D-15
1013 0 : xin(net_iso(ife51)) = 6.5787879505834196D-15
1014 0 : xin(net_iso(io17)) = 3.8661065919908615D-15
1015 0 : xin(net_iso(ic13)) = 2.2604411883637647D-15
1016 0 : xin(net_iso(icr58)) = 2.1539938862813166D-15
1017 0 : xin(net_iso(isi33)) = 1.4888165334138879D-15
1018 0 : xin(net_iso(ina24)) = 7.2312125147882022D-16
1019 0 : xin(net_iso(ial25)) = 5.3576177397850227D-16
1020 0 : xin(net_iso(ico65)) = 4.9566812556376405D-16
1021 0 : xin(net_iso(icr47)) = 3.0194388465619186D-16
1022 0 : xin(net_iso(in14)) = 2.6891785216435022D-16
1023 0 : xin(net_iso(ina22)) = 2.5931749891034273D-16
1024 0 : xin(net_iso(io15)) = 2.5261003710407275D-16
1025 0 : xin(net_iso(ini68)) = 2.2419710841076782D-16
1026 0 : xin(net_iso(io18)) = 1.8206042351246347D-16
1027 0 : xin(net_iso(iti43)) = 9.9895985562055471D-17
1028 0 : xin(net_iso(if19)) = 7.0445166521693986D-17
1029 0 : xin(net_iso(ihe3)) = 6.8591249543794866D-17
1030 0 : xin(net_iso(iti54)) = 3.3955810406819116D-17
1031 0 : xin(net_iso(ife63)) = 3.0222399040043984D-17
1032 0 : xin(net_iso(in13)) = 2.5097133179904447D-17
1033 0 : xin(net_iso(izn59)) = 2.3782224732332168D-17
1034 0 : xin(net_iso(img23)) = 2.2784900370197838D-17
1035 0 : xin(net_iso(if17)) = 1.0052207505944401D-17
1036 0 : xin(net_iso(if18)) = 6.3013547340360942D-18
1037 0 : xin(net_iso(iv56)) = 2.1677062516769839D-18
1038 0 : xin(net_iso(ina21)) = 1.9109919103480182D-18
1039 0 : xin(net_iso(ine23)) = 1.2143041459039416D-18
1040 0 : xin(net_iso(ili6)) = 1.4203908924439747D-19
1041 0 : xin(net_iso(ib11)) = 1.1320699694157424D-19
1042 0 : xin(net_iso(if20)) = 4.4936577179455616D-20
1043 0 : xin(net_iso(ine19)) = 2.6387620496069204D-20
1044 0 : xin(net_iso(ife64)) = 1.1841283303587735D-20
1045 0 : xin(net_iso(ib10)) = 1.0494357182634838D-20
1046 0 : xin(net_iso(ibe9)) = 7.8292004126082962D-21
1047 0 : xin(net_iso(ico66)) = 6.4514234785580282D-21
1048 0 : xin(net_iso(ili7)) = 6.4231341717191173D-21
1049 0 : xin(net_iso(in16)) = 4.7443318701592232D-21
1050 0 : xin(net_iso(ibe7)) = 3.3532048357084064D-22
1051 0 : xin(net_iso(io19)) = 3.2062683754970905D-22
1052 0 : xin(net_iso(ibe10)) = 6.5331631260779713D-23
1053 0 : xin(net_iso(ico67)) = 4.9629777598135805D-24
1054 0 : xin(net_iso(ife65)) = 3.7335326045384740D-26
1055 0 : xin(net_iso(ife66)) = 8.6772104841406824D-30
1056 0 : xin(net_iso(ib8)) = 4.1439050548710376D-31
1057 :
1058 0 : write (*, *) 'sum xin', sum(xin(:))
1059 :
1060 0 : logT = 9.6532818288064650D+00
1061 0 : logRho = 7.9479966082179185D+00
1062 0 : eta = 2.7403163311838425D+00
1063 :
1064 0 : screening_mode = extended_screening
1065 :
1066 0 : call net_set_logTcut(net_handle, 0d0, 0d0, info)
1067 0 : if (info /= 0) then
1068 0 : write (*, *) 'failed in net_set_logTcut'
1069 0 : call mesa_error(__FILE__, __LINE__)
1070 : end if
1071 :
1072 0 : else if (net_file == 'approx21_cr60_plus_co56.net') then
1073 :
1074 0 : nrates_to_show = 2
1075 :
1076 0 : rates_to_show(1:nrates_to_show) = [ &
1077 : irco56ec_to_fe56, &
1078 0 : irni56ec_to_co56]
1079 :
1080 0 : xin(net_iso(ineut)) = 1.9615092621881698D-06
1081 0 : xin(net_iso(ih1)) = 0.0000000000000000D+00
1082 0 : xin(net_iso(iprot)) = 3.3544919533130298D-04
1083 0 : xin(net_iso(ihe3)) = 0.0000000000000000D+00
1084 0 : xin(net_iso(ihe4)) = 9.6491906054115388D-03
1085 0 : xin(net_iso(ic12)) = 4.4455299753403973D-07
1086 0 : xin(net_iso(in14)) = 0.0000000000000000D+00
1087 0 : xin(net_iso(io16)) = 7.0439210492179401D-07
1088 0 : xin(net_iso(ine20)) = 4.5369892002182717D-09
1089 0 : xin(net_iso(img24)) = 8.0018732551576689D-07
1090 0 : xin(net_iso(isi28)) = 3.3697411043955434D-04
1091 0 : xin(net_iso(is32)) = 2.7239844165292807D-04
1092 0 : xin(net_iso(iar36)) = 1.4240894544492388D-04
1093 0 : xin(net_iso(ica40)) = 1.5156321095874170D-04
1094 0 : xin(net_iso(iti44)) = 4.3979509377388036D-06
1095 0 : xin(net_iso(icr48)) = 2.8113145151721771D-05
1096 0 : xin(net_iso(icr60)) = 4.1810609055173498D-03
1097 0 : xin(net_iso(ife52)) = 1.8928784597602586D-04
1098 0 : xin(net_iso(ife54)) = 2.6996971418369603D-01
1099 0 : xin(net_iso(ife56)) = 7.0702963220409232D-01
1100 0 : xin(net_iso(ico56)) = 6.8333792315189530D-03
1101 0 : xin(net_iso(ini56)) = 8.7251484519146338D-04
1102 :
1103 0 : write (*, *) 'sum xin', sum(xin(:))
1104 :
1105 0 : logT = 9.8200000000000003D+00
1106 0 : logRho = 8.2586740078478176D+00
1107 0 : eta = 0d0
1108 :
1109 0 : screening_mode = extended_screening
1110 :
1111 0 : call net_set_logTcut(net_handle, 0d0, 0d0, info)
1112 0 : if (info /= 0) then
1113 0 : write (*, *) 'failed in net_set_logTcut'
1114 0 : call mesa_error(__FILE__, __LINE__)
1115 : end if
1116 :
1117 0 : else if (net_file == 'approx21_cr60_plus_fe53_fe55_co56.net') then
1118 :
1119 0 : nrates_to_show = 2
1120 :
1121 0 : rates_to_show(1:nrates_to_show) = [ &
1122 : irco56ec_to_fe56, &
1123 0 : irni56ec_to_co56]
1124 :
1125 0 : xin(net_iso(ihe4)) = 3.4555392534813939D-01
1126 0 : xin(net_iso(io16)) = 1.9367778420430937D-01
1127 0 : xin(net_iso(ih1)) = 1.9052931367172501D-01
1128 0 : xin(net_iso(isi28)) = 9.1023936032064601D-02
1129 0 : xin(net_iso(ini56)) = 5.8492459653295005D-02
1130 0 : xin(net_iso(is32)) = 4.1416642476999908D-02
1131 0 : xin(net_iso(ine20)) = 2.0927782332477558D-02
1132 0 : xin(net_iso(ic12)) = 2.0589617104448312D-02
1133 0 : xin(net_iso(img24)) = 1.8975914108406104D-02
1134 0 : xin(net_iso(iar36)) = 7.0600338785939956D-03
1135 0 : xin(net_iso(ica40)) = 4.9182171981353726D-03
1136 0 : xin(net_iso(ife52)) = 3.1618820806005280D-03
1137 0 : xin(net_iso(in14)) = 1.9878460082760952D-03
1138 0 : xin(net_iso(ife56)) = 7.1668776621130190D-04
1139 0 : xin(net_iso(ico56)) = 4.4974971099921181D-04
1140 0 : xin(net_iso(ife54)) = 3.3432423890988422D-04
1141 0 : xin(net_iso(icr48)) = 1.7602626907769762D-04
1142 0 : xin(net_iso(ihe3)) = 5.1650097191631545D-06
1143 0 : xin(net_iso(iti44)) = 2.6452422203389906D-06
1144 0 : xin(net_iso(icr60)) = 4.7653353095937982D-08
1145 0 : xin(net_iso(iprot)) = 1.2739022407246250D-11
1146 0 : xin(net_iso(ife53)) = 0.0000000000000000D+00
1147 0 : xin(net_iso(ife55)) = 0.0000000000000000D+00
1148 0 : xin(net_iso(ineut)) = 0.0000000000000000D+00
1149 :
1150 0 : write (*, *) 'sum xin', sum(xin(:))
1151 :
1152 0 : logT = 4.6233007922659333D+00
1153 0 : logRho = -1.0746410107891649D+01
1154 0 : eta = -2.2590260158215202D+01
1155 :
1156 0 : screening_mode = extended_screening
1157 :
1158 0 : call net_set_logTcut(net_handle, 0d0, 0d0, info)
1159 0 : if (info /= 0) then
1160 0 : write (*, *) 'failed in net_set_logTcut'
1161 0 : call mesa_error(__FILE__, __LINE__)
1162 : end if
1163 :
1164 0 : else if (net_file == 'basic.net') then
1165 :
1166 0 : nrates_to_show = 1
1167 :
1168 0 : if (rates_reaction_id('rc12_to_n14') <= 0) call mesa_error(__FILE__, __LINE__, 'bad reaction')
1169 0 : write (*, *) 'rc12_to_n14', rates_reaction_id('rc12_to_n14')
1170 :
1171 0 : rates_to_show(1:nrates_to_show) = [ &
1172 0 : rates_reaction_id('rc12_to_n14')]
1173 :
1174 0 : xin(net_iso(ihe4)) = 9.8119124177708650D-01
1175 0 : xin(net_iso(in14)) = 9.8369547495994036D-03
1176 0 : xin(net_iso(io16)) = 2.9223115895360822D-03
1177 0 : xin(net_iso(ine20)) = 2.0337034688681288D-03
1178 0 : xin(net_iso(ihe3)) = 0.0000000000000000D+00
1179 0 : xin(net_iso(ih1)) = 0.0000000000000000D+00
1180 :
1181 0 : write (*, *) 'when 1st'
1182 0 : xin(net_iso(ic12)) = 2.3551202768735954D-04 ! 2.3551179217556737D-04
1183 0 : xin(net_iso(img24)) = 3.7802763872225075D-03 ! 3.7802766227342998D-03
1184 :
1185 : !write(*,*) 'when 2nd'
1186 : !xin(net_iso(ic12))= 2.3551179217556737D-04
1187 : !xin(net_iso(img24))= 3.7802766227342998D-03
1188 :
1189 0 : logT = 7.8820961200821831D+00
1190 0 : logRho = 5.6124502722887746D+00
1191 0 : eta = 1.2629003275927920D+01
1192 :
1193 0 : write (*, *) 'sum xin', sum(xin(:))
1194 :
1195 0 : screening_mode = extended_screening
1196 :
1197 0 : call net_set_logTcut(net_handle, 0d0, 0d0, info)
1198 0 : if (info /= 0) then
1199 0 : write (*, *) 'failed in net_set_logTcut'
1200 0 : call mesa_error(__FILE__, __LINE__)
1201 : end if
1202 :
1203 0 : else if (net_file == 'agb.net') then
1204 :
1205 0 : nrates_to_show = 4
1206 :
1207 0 : rates_to_show(1:nrates_to_show) = [ &
1208 : ir_h1_h1_wk_h2, &
1209 : ir_c13_an_o16, &
1210 : ir_f19_ap_ne22, &
1211 0 : ir_he3_ag_be7]
1212 :
1213 0 : xin(net_iso(ih1)) = 1
1214 :
1215 0 : write (*, *) 'sum xin', sum(xin(:))
1216 :
1217 0 : logT = 8.6864273893515023D+00
1218 0 : logRho = 2.0591020210828619D+00
1219 0 : eta = -1.4317150417353590D+01
1220 :
1221 0 : screening_mode = extended_screening
1222 :
1223 0 : call net_set_logTcut(net_handle, 0d0, 0d0, info)
1224 0 : if (info /= 0) then
1225 0 : write (*, *) 'failed in net_set_logTcut'
1226 0 : call mesa_error(__FILE__, __LINE__)
1227 : end if
1228 :
1229 : if (.false.) then
1230 : if (.false.) then
1231 : rate_factors(:) = 0
1232 : i = reaction_table(ir_h2_pg_he3)
1233 : if (i == 0) call mesa_error(__FILE__, __LINE__)
1234 : rate_factors(i) = 1
1235 : else
1236 : i = reaction_table(ir_ne20_ag_mg24)
1237 : if (i == 0) call mesa_error(__FILE__, __LINE__)
1238 : rate_factors(i) = 0
1239 : i = reaction_table(ir_o16_ag_ne20)
1240 : if (i == 0) call mesa_error(__FILE__, __LINE__)
1241 : rate_factors(i) = 0
1242 : end if
1243 : end if
1244 :
1245 0 : else if (net_file == 'pp_and_cno_extras.net') then
1246 :
1247 0 : nrates_to_show = 8
1248 :
1249 0 : rates_to_show(1:nrates_to_show) = [ &
1250 : rates_reaction_id('r_n13_wk_c13'), &
1251 : rates_reaction_id('r_o15_wk_n15'), &
1252 : rates_reaction_id('r_f17_wk_o17'), &
1253 : rates_reaction_id('r_f18_wk_o18'), &
1254 : rates_reaction_id('r_o14_wk_n14'), &
1255 : rates_reaction_id('r_ne18_wk_f18'), &
1256 : rates_reaction_id('r_ne19_wk_f19'), &
1257 0 : ir_he4_he4_he4_to_c12]
1258 :
1259 0 : xin = 0
1260 0 : xin(net_iso(ih1)) = 7.2265805432969643D-01
1261 0 : xin(net_iso(ihe3)) = 6.7801726921522655D-04
1262 0 : xin(net_iso(ihe4)) = 2.6667042876319019D-01
1263 0 : xin(net_iso(ic12)) = 1.9056849943017622D-03
1264 0 : xin(net_iso(in13)) = 1.4791107757148081D-04
1265 0 : xin(net_iso(in14)) = 6.0253770632534619D-04
1266 0 : xin(net_iso(in15)) = 1.9263919190065488D-07
1267 0 : xin(net_iso(io14)) = 9.5977897394247582D-09
1268 0 : xin(net_iso(io15)) = 7.6666060610240002D-08
1269 0 : xin(net_iso(io16)) = 5.5886173952684358D-03
1270 0 : xin(net_iso(io17)) = 2.0449560316023342D-06
1271 0 : xin(net_iso(io18)) = 1.1313355448541673D-05
1272 0 : xin(net_iso(if19)) = 2.1343308067891499D-07
1273 0 : xin(net_iso(ine20)) = 1.2563102679570999D-03
1274 0 : xin(net_iso(img24)) = 4.7858754879924638D-04
1275 :
1276 0 : logT = 9.6d0
1277 0 : logRho = 6.0d0
1278 0 : rho = 7.8571498592117219D+00
1279 0 : T = 8.5648111120065376D+06
1280 0 : abar = 1.2655060647252907D+00
1281 0 : zbar = 1.0901664301076275D+00
1282 0 : z2bar = 1.3036906023574921D+00
1283 0 : ye = 8.6144702146826535D-01
1284 0 : eta = -3.4387570967781595D+00
1285 :
1286 0 : screening_mode = extended_screening
1287 :
1288 0 : else if (net_file == 'approx21_plus_co56.net') then
1289 :
1290 0 : nrates_to_show = 5
1291 :
1292 0 : rates_to_show(1:nrates_to_show) = [ &
1293 : ir_v47_pa_ti44, &
1294 : ir_mn51_pa_cr48, &
1295 : ir_ar36_ap_k39, &
1296 : ir_co55_pa_fe52, &
1297 : ir_s32_ap_cl35 &
1298 0 : ]
1299 0 : xin = 0
1300 :
1301 0 : xin(net_iso(ife54)) = 7.8234742556602999D-01
1302 0 : xin(net_iso(isi28)) = 7.8210084821085060D-02
1303 0 : xin(net_iso(ini56)) = 5.2306555890846963D-02
1304 0 : xin(net_iso(is32)) = 5.1718599300602117D-02
1305 0 : xin(net_iso(iar36)) = 1.3861579457866108D-02
1306 0 : xin(net_iso(ica40)) = 1.2347894507489101D-02
1307 0 : xin(net_iso(ife56)) = 4.8313904464514874D-03
1308 0 : xin(net_iso(ife52)) = 3.9677414859767722D-03
1309 0 : xin(net_iso(icr48)) = 2.9870241199304463D-04
1310 0 : xin(net_iso(ihe4)) = 4.0616997047809087D-05
1311 0 : xin(net_iso(iti44)) = 3.3724322385639036D-05
1312 0 : xin(net_iso(iprot)) = 1.7517084938626562D-05
1313 0 : xin(net_iso(img24)) = 1.1248868356795997D-05
1314 0 : xin(net_iso(io16)) = 6.1628768742062784D-06
1315 0 : xin(net_iso(ic12)) = 7.4857838132808694D-07
1316 0 : xin(net_iso(ine20)) = 5.6045205217740520D-09
1317 0 : xin(net_iso(icr56)) = 1.7593106525983871D-09
1318 0 : xin(net_iso(ineut)) = 1.9978632764577629D-11
1319 0 : xin(net_iso(in14)) = 0.0000000000000000D+00
1320 0 : xin(net_iso(ihe3)) = 0.0000000000000000D+00
1321 0 : xin(net_iso(ih1)) = 0.0000000000000000D+00
1322 0 : write (*, *) 'test case sum xin', sum(xin(1:species))
1323 :
1324 0 : logT = 9.5806070583042597D+00
1325 0 : logRho = 7.1251356937727763D+00
1326 0 : eta = 6.3504500633747440D-01
1327 :
1328 0 : T = 3.8072119794259882D+09
1329 0 : rho = 1.3339381513098311D+07
1330 0 : abar = 4.8254908015575154D+01
1331 0 : zbar = 2.3420437256420026D+01
1332 0 : z2bar = 5.7180065624152712D+02
1333 0 : ye = 4.8534829345982072D-01
1334 0 : screening_mode = extended_screening
1335 :
1336 : else
1337 :
1338 0 : write (*, *) 'need to define setup for net_file ', trim(net_file)
1339 0 : call mesa_error(__FILE__, __LINE__, 'Do_One_Test')
1340 :
1341 : end if
1342 :
1343 0 : Rho = exp10(logRho)
1344 0 : T = exp10(logT)
1345 :
1346 0 : write (*, *)
1347 0 : write (*, *)
1348 :
1349 0 : info = 0
1350 :
1351 0 : ierr = 0
1352 : call composition_info( &
1353 : species, chem_id, xin, xh, xhe, z, abar, zbar, z2bar, z53bar, ye, &
1354 0 : mass_correction, xsum, dabar_dx, dzbar_dx, dmc_dx)
1355 :
1356 0 : write (*, '(a40,d26.16)') 'xh', xh
1357 0 : write (*, '(a40,d26.16)') 'xhe', xhe
1358 0 : write (*, '(a40,d26.16)') 'abar', abar
1359 0 : write (*, '(a40,d26.16)') 'zbar', zbar
1360 0 : write (*, '(a40,d26.16)') 'z2bar', z2bar
1361 0 : write (*, '(a40,d26.16)') 'ye', ye
1362 0 : do i = 1, species
1363 0 : write (*, '(a40,i6,d26.16)') 'init x '//trim(chem_isos%name(chem_id(i))), i, xin(i)
1364 : end do
1365 0 : write (*, '(A)')
1366 0 : write (*, '(a40,d26.16)') 'logT', logT
1367 0 : write (*, '(a40,d26.16)') 'logRho', logRho
1368 0 : write (*, '(a40,d26.16)') 'eta', eta
1369 :
1370 0 : skip_jacobian = .false.
1371 :
1372 0 : n%screening_mode = screening_mode
1373 :
1374 : if (.false.) then
1375 : write (*, *) 'call net_get_rates_only'
1376 : call net_get_rates_only( &
1377 : net_handle, n, species, num_reactions, &
1378 : xin, T, logT, Rho, logRho, &
1379 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
1380 : rate_factors, weak_rate_factor, &
1381 : std_reaction_Qs, std_reaction_neuQs, &
1382 : screening_mode, &
1383 : ierr)
1384 : call mesa_error(__FILE__, __LINE__, 'net_get_rates_only')
1385 : end if
1386 :
1387 : call net_get_with_Qs( &
1388 : net_handle, skip_jacobian, n, species, num_reactions, &
1389 : xin, T, logT, Rho, logRho, &
1390 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
1391 : rate_factors, weak_rate_factor, &
1392 : std_reaction_Qs, std_reaction_neuQs, &
1393 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
1394 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
1395 : screening_mode, &
1396 : eps_nuc_categories, eps_neu_total, &
1397 0 : actual_Qs, actual_neuQs, from_weaklib, info)
1398 0 : if (info /= 0) then
1399 0 : write (*, 1) 'logT', logT
1400 0 : write (*, 1) 'logRho', logRho
1401 0 : write (*, *) 'bad return from net_get'
1402 0 : call mesa_error(__FILE__, __LINE__)
1403 : end if
1404 :
1405 : if (.true.) then ! dfridr tests for partials
1406 :
1407 : !doing_d_dlnd = .true.
1408 0 : doing_d_dlnd = .false.
1409 0 : doing_dx = .false.
1410 :
1411 0 : j_dx = 22
1412 0 : var_0 = dxdt(j_dx)
1413 :
1414 0 : if (doing_dx) then
1415 0 : j_dx = 20 ! fe56
1416 0 : j_dx_sink = 17 ! cr56
1417 0 : dx_0 = xin(j_dx)*1d-6
1418 0 : dvardx_0 = d_eps_nuc_dx(j_dx)
1419 0 : write (*, 1) 'xin(fe56)', xin(j_dx)
1420 0 : write (*, 1) 'xin(cr56)', xin(j_dx_sink)
1421 0 : write (*, 1) 'eps_nuc', eps_nuc
1422 0 : write (*, '(A)')
1423 : !call mesa_error(__FILE__,__LINE__,'testing')
1424 0 : else if (doing_d_dlnd) then
1425 0 : dx_0 = max(1d-14, abs(logRho*ln10*1d-6))
1426 0 : dvardx_0 = d_eps_nuc_dRho*Rho ! d_dxdt_dRho(1)*Rho ! analytic value of partial
1427 : else
1428 0 : dx_0 = max(1d-14, abs(logT*ln10*1d-6))
1429 : !dvardx_0 = d_eps_nuc_dT*T
1430 0 : dvardx_0 = d_dxdt_dT(j_dx)*T
1431 : end if
1432 0 : err = 0d0
1433 0 : dvardx = dfridr(dx_0, err)
1434 0 : xdum = (dvardx - dvardx_0)/max(abs(dvardx_0), 1d-50)
1435 0 : write (*, 1) 'analytic, numeric, est err in numeric, rel diff', &
1436 0 : dvardx_0, dvardx, err, xdum
1437 0 : if (doing_dx) then
1438 0 : write (*, *) 'doing d_dx '//trim(chem_isos%name(chem_id(j_dx)))
1439 0 : else if (doing_d_dlnd) then
1440 0 : write (*, *) 'doing dlnd'
1441 : else ! doing d_dlnT
1442 0 : write (*, *) 'doing dlnT'
1443 : end if
1444 0 : write (*, *) 'test net'
1445 0 : write (*, '(A)')
1446 :
1447 0 : call mesa_error(__FILE__, __LINE__, 'test net')
1448 :
1449 : end if
1450 :
1451 0 : if (show_Qs) then
1452 0 : write (*, '(A)')
1453 0 : write (*, 1) 'logT', logT
1454 0 : write (*, 1) 'logRho', logRho
1455 0 : write (*, '(A)')
1456 0 : write (*, '(30x,4a20)') 'Q total', 'Q neutrino', 'Q total-neutrino'
1457 0 : do i = 1, num_reactions
1458 0 : if (from_weaklib(i)) then
1459 0 : write (*, '(a30,99f20.10)') 'weaklib '//trim(reaction_Name(reaction_id(i))), &
1460 0 : actual_Qs(i), actual_neuQs(i), actual_Qs(i) - actual_neuQs(i)
1461 : else
1462 0 : write (*, '(a30,99f20.10)') trim(reaction_Name(reaction_id(i))), &
1463 0 : actual_Qs(i), actual_neuQs(i), actual_Qs(i) - actual_neuQs(i)
1464 : end if
1465 : end do
1466 0 : write (*, '(A)')
1467 0 : stop
1468 : end if
1469 :
1470 0 : write (*, 2) 'screening_mode', screening_mode
1471 :
1472 : if (.true.) then
1473 0 : write (*, 1) 'logT', logT
1474 0 : write (*, 1) 'logRho', logRho
1475 0 : write (*, '(A)')
1476 :
1477 0 : write (*, 1) 'eps_nuc', eps_nuc
1478 0 : write (*, 1) 'd_epsnuc_dlnd', d_eps_nuc_dRho*Rho
1479 0 : write (*, 1) 'd_epsnuc_dlnT', d_eps_nuc_dT*T
1480 0 : write (*, '(A)')
1481 :
1482 0 : if (eps_nuc > 0) then
1483 0 : write (*, 1) 'log eps_nuc', log10(eps_nuc)
1484 0 : write (*, 1) 'd_lnepsnuc_dlnd', d_eps_nuc_dRho*Rho/eps_nuc
1485 0 : write (*, 1) 'd_lnepsnuc_dlnT', d_eps_nuc_dT*T/eps_nuc
1486 0 : write (*, '(A)')
1487 : end if
1488 :
1489 : !stop
1490 :
1491 : if (.false.) then
1492 :
1493 : do i = 1, species
1494 : write (*, 1) 'd_eps_nuc_dx '//trim(chem_isos%name(chem_id(i))), d_eps_nuc_dx(i)
1495 : end do
1496 : write (*, '(A)')
1497 :
1498 : do i = 1, species
1499 : write (*, 1) 'd_dxdt_dlnRho '//trim(chem_isos%name(chem_id(i))), d_dxdt_dRho(i)*rho
1500 : end do
1501 : write (*, '(A)')
1502 :
1503 : do i = 1, species
1504 : write (*, 1) 'd_dxdt_dlnT '//trim(chem_isos%name(chem_id(i))), d_dxdt_dT(i)*T
1505 : end do
1506 : write (*, '(A)')
1507 :
1508 : if (.false.) then
1509 : do i = 1, species
1510 : write (*, 1) 'd_dxdt_dx(:,neut) '// &
1511 : trim(chem_isos%name(chem_id(i))), d_dxdt_dx(i, net_iso(ineut))
1512 : end do
1513 : write (*, '(A)')
1514 : end if
1515 :
1516 : do i = 1, species
1517 : write (*, 1) 'dxdt '//trim(chem_isos%name(chem_id(i))), dxdt(i)
1518 : end do
1519 : write (*, 1) 'sum(dxdt)', sum(dxdt(1:species))
1520 :
1521 : do i = 1, nrates_to_show
1522 : j = rates_to_show(i)
1523 : if (j == 0) cycle
1524 : write (*, 1) 'd_rate_raw_dT '//trim(reaction_Name(j)), &
1525 : rate_raw_dT(reaction_table(j))
1526 : end do
1527 : write (*, '(A)')
1528 :
1529 : end if
1530 :
1531 0 : do i = 1, nrates_to_show
1532 0 : j = rates_to_show(i)
1533 0 : if (j == 0) cycle
1534 0 : write (*, 1) 'rate_raw '//trim(reaction_Name(j)), &
1535 0 : rate_raw(reaction_table(j))
1536 : end do
1537 0 : write (*, '(A)')
1538 :
1539 0 : do i = 1, nrates_to_show
1540 0 : j = rates_to_show(i)
1541 0 : if (j == 0) cycle
1542 0 : write (*, 1) 'rate_screened '//trim(reaction_Name(j)), &
1543 0 : rate_screened(reaction_table(j))
1544 : end do
1545 0 : write (*, '(A)')
1546 :
1547 0 : do i = 1, nrates_to_show
1548 0 : j = rates_to_show(i)
1549 0 : if (j == 0) cycle
1550 0 : write (*, 1) 'Q total '//trim(reaction_Name(j)), &
1551 0 : actual_Qs(reaction_table(j))
1552 : end do
1553 0 : write (*, '(A)')
1554 :
1555 0 : do i = 1, nrates_to_show
1556 0 : j = rates_to_show(i)
1557 0 : if (j == 0) cycle
1558 0 : write (*, 1) 'Q neutrino '//trim(reaction_Name(j)), &
1559 0 : actual_neuQs(reaction_table(j))
1560 : end do
1561 0 : write (*, '(A)')
1562 :
1563 0 : do i = 1, nrates_to_show
1564 0 : j = rates_to_show(i)
1565 0 : if (j == 0) cycle
1566 0 : write (*, 1) 'Q total-neutrino '//trim(reaction_Name(j)), &
1567 0 : actual_Qs(reaction_table(j)) - actual_neuQs(reaction_table(j))
1568 : end do
1569 0 : write (*, '(A)')
1570 :
1571 : if (.false.) then
1572 :
1573 : do i = 1, species
1574 : write (*, 1) 'x '//trim(chem_isos%name(chem_id(i))), xin(i)
1575 : end do
1576 : write (*, '(A)')
1577 :
1578 : do i = 1, species
1579 : if (-dxdt(i) > 1d-90) &
1580 : write (*, 1) 'x/dxdt '//trim(chem_isos%name(chem_id(i))), xin(i)/dxdt(i)
1581 : end do
1582 : write (*, '(A)')
1583 :
1584 : do i = 1, num_categories
1585 : if (abs(eps_nuc_categories(i)) < 1d-20) cycle
1586 : write (*, 1) 'eps_nuc_cat '//trim(category_name(i)), eps_nuc_categories(i)
1587 : end do
1588 : write (*, '(A)')
1589 :
1590 : end if
1591 :
1592 0 : stop
1593 : end if
1594 :
1595 : write (*, '(A)')
1596 : write (*, '(A)')
1597 : write (*, *) 'net_name ', trim(net_file)
1598 : write (*, *) 'species', species
1599 : write (*, 1) 'abar =', abar
1600 : write (*, 1) 'zbar =', zbar
1601 : write (*, 1) 'z2bar =', z2bar
1602 : write (*, 1) 'ye =', ye
1603 : write (*, *)
1604 : do i = 1, nrates_to_show
1605 : j = rates_to_show(i)
1606 : if (j == 0) cycle
1607 : if (reaction_table(j) == 0) then
1608 : write (*, *) 'missing reaction_table(j) for '//trim(reaction_Name(j))
1609 : stop
1610 : end if
1611 : end do
1612 : write (*, '(A)')
1613 : write (*, 1) 'eps_nuc', eps_nuc
1614 : write (*, 1) 'd_eps_nuc_dRho', d_eps_nuc_dRho
1615 : write (*, 1) 'd_eps_nuc_dT', d_eps_nuc_dT
1616 : do j = 1, species
1617 : write (*, 1) 'd_eps_nuc_dx '//trim(chem_isos%name(chem_id(j))), d_eps_nuc_dx(j)
1618 : end do
1619 : write (*, '(A)')
1620 : do j = 1, species
1621 : write (*, 1) 'dxdt '//trim(chem_isos%name(chem_id(j))), dxdt(j)
1622 : end do
1623 : write (*, '(A)')
1624 : do j = 1, species
1625 : write (*, 1) 'd_dxdt_dRho '//trim(chem_isos%name(chem_id(j))), d_dxdt_dRho(j)
1626 : end do
1627 : write (*, '(A)')
1628 : do j = 1, species
1629 : write (*, 1) 'd_dxdt_dT '//trim(chem_isos%name(chem_id(j))), d_dxdt_dT(j)
1630 : end do
1631 : write (*, '(A)')
1632 : do j = 1, species
1633 : write (*, 1) 'd_dxdt_dx(1,:) '//trim(chem_isos%name(chem_id(j))), d_dxdt_dx(1, j)
1634 : end do
1635 : write (*, *)
1636 : do j = 1, num_categories
1637 : write (*, 1) trim(category_name(j)), eps_nuc_categories(j)
1638 : end do
1639 : write (*, '(A)')
1640 : write (*, 1) 'eta =', eta
1641 : write (*, 1) 'logT =', logT
1642 : write (*, 1) 'logRho =', logRho
1643 : write (*, *) 'screening_mode =', screening_mode
1644 : write (*, '(A)')
1645 : do j = 1, species
1646 : write (*, 1) 'xin(net_iso(i'//trim(chem_isos%name(chem_id(j)))//'))= ', xin(j)
1647 : end do
1648 : write (*, '(A)')
1649 : write (*, 1) 'sum(xin(1:species))', sum(xin(1:species))
1650 : write (*, 1) '1 - sum(xin(1:species))', 1 - sum(xin(1:species))
1651 : write (*, '(A)')
1652 : write (*, 1) 'eps_nuc', eps_nuc
1653 : write (*, 1) 'eps_nuc_neu_total', eps_neu_total
1654 :
1655 0 : deallocate (rate_factors, actual_Qs, actual_neuQs, from_weaklib)
1656 :
1657 : contains
1658 :
1659 0 : real(dp) function dfridr_func(delta_x) result(val)
1660 : real(dp), intent(in) :: delta_x
1661 : integer :: ierr
1662 0 : real(dp) :: var, log_var
1663 : include 'formats'
1664 0 : ierr = 0
1665 :
1666 0 : if (doing_dx) then
1667 :
1668 0 : xin_copy = xin
1669 0 : xin_copy(j_dx) = xin_copy(j_dx) + delta_x
1670 0 : xin_copy(j_dx_sink) = xin_copy(j_dx_sink) - delta_x
1671 :
1672 : call net_get_with_Qs( &
1673 : net_handle, skip_jacobian, n, species, num_reactions, &
1674 : xin_copy, T, logT, Rho, logRho, &
1675 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
1676 : rate_factors, weak_rate_factor, &
1677 : std_reaction_Qs, std_reaction_neuQs, &
1678 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
1679 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
1680 : screening_mode, &
1681 : eps_nuc_categories, eps_neu_total, &
1682 0 : actual_Qs, actual_neuQs, from_weaklib, ierr)
1683 :
1684 0 : write (*, 1) 'xin(j)', xin_copy(j_dx)
1685 0 : write (*, 1) 'xin(j_sink)', xin_copy(j_dx_sink)
1686 0 : write (*, 1) 'eps_nuc', eps_nuc
1687 0 : write (*, '(A)')
1688 :
1689 0 : else if (doing_d_dlnd) then
1690 :
1691 0 : log_var = logRho + delta_x/ln10
1692 0 : var = exp10(log_var)
1693 :
1694 : call net_get_with_Qs( &
1695 : net_handle, skip_jacobian, n, species, num_reactions, &
1696 : xin, T, logT, var, log_var, &
1697 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
1698 : rate_factors, weak_rate_factor, &
1699 : std_reaction_Qs, std_reaction_neuQs, &
1700 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
1701 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
1702 : screening_mode, &
1703 : eps_nuc_categories, eps_neu_total, &
1704 0 : actual_Qs, actual_neuQs, from_weaklib, ierr)
1705 :
1706 : else
1707 :
1708 0 : log_var = logT + delta_x/ln10
1709 0 : var = exp10(log_var)
1710 :
1711 : call net_get_with_Qs( &
1712 : net_handle, skip_jacobian, n, species, num_reactions, &
1713 : xin, var, log_var, Rho, logRho, &
1714 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
1715 : rate_factors, weak_rate_factor, &
1716 : std_reaction_Qs, std_reaction_neuQs, &
1717 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
1718 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
1719 : screening_mode, &
1720 : eps_nuc_categories, eps_neu_total, &
1721 0 : actual_Qs, actual_neuQs, from_weaklib, ierr)
1722 :
1723 : end if
1724 :
1725 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'failed in call on net_get_with_Qs')
1726 : !val = eps_nuc ! dxdt(1)
1727 0 : val = dxdt(j_dx)
1728 :
1729 0 : end function dfridr_func
1730 :
1731 0 : real(dp) function dfridr(hx, err) ! from Frank
1732 : real(dp), intent(in) :: hx
1733 : real(dp), intent(out) :: err
1734 : ! this routine returns the first derivative of a function func(x)
1735 : ! at the point x, by ridders method of polynomial extrapolation.
1736 : ! value hx is the initial step size;
1737 : ! it should be an increment for which func changes substantially.
1738 : ! an estimate of the error in the first derivative is returned in err.
1739 : integer, parameter :: ntab = 20
1740 : integer :: i, j
1741 0 : real(dp) :: errt, fac, hh, a(ntab, ntab)
1742 : real(dp), parameter :: con2 = 2d0, con = sqrt(con2), big = 1d50, safe = 2d0
1743 : include 'formats'
1744 0 : dfridr = 0d0
1745 0 : hh = hx
1746 : ! 2nd order central difference
1747 0 : a(1, 1) = (dfridr_func(hh) - dfridr_func(-hh))/(2d0*hh)
1748 0 : write (*, 2) 'dfdx hh', 1, a(1, 1), hh
1749 0 : err = big
1750 : ! successive columns in the neville tableu will go to smaller stepsizes
1751 : ! and higher orders of extrapolation
1752 0 : do i = 2, ntab
1753 0 : hh = hh/con
1754 0 : a(1, i) = (dfridr_func(hh) - dfridr_func(-hh))/(2d0*hh)
1755 0 : write (*, 2) 'dfdx hh', i, a(1, i), hh
1756 : ! compute extrapolations of various orders; the error strategy is to compare
1757 : ! each new extrapolation to one order lower but both at the same stepsize
1758 : ! and at the previous stepsize
1759 0 : fac = con2
1760 0 : do j = 2, i
1761 0 : a(j, i) = (a(j - 1, i)*fac - a(j - 1, i - 1))/(fac - 1d0)
1762 0 : fac = con2*fac
1763 0 : errt = max(abs(a(j, i) - a(j - 1, i)), abs(a(j, i) - a(j - 1, i - 1)))
1764 : ! write(*,3) 'a(j,i)', j, i, a(j,i), errt
1765 0 : if (errt <= err) then
1766 0 : err = errt
1767 0 : dfridr = a(j, i)
1768 : !write(*,3) 'dfridr err', i, j, dfridr, err
1769 : end if
1770 : end do
1771 : ! if higher order is worse by a significant factor safe, then bail
1772 0 : if (abs(a(i, i) - a(i - 1, i - 1)) >= safe*err) then
1773 0 : write (*, 1) 'higher order is worse', err, a(i, i), a(i - 1, i - 1)
1774 0 : return
1775 : end if
1776 : end do
1777 0 : end function dfridr
1778 :
1779 : end subroutine Do_One_Testcase
1780 :
1781 0 : subroutine test_neutrino_Q
1782 : real(dp), parameter :: Qnu_n13 = 0.714440d0 !..13n(e+nu)13c
1783 : real(dp), parameter :: Qnu_o15 = 1.005513d0 !..15o(e+nu)15n
1784 : real(dp), parameter :: Qnu_f17 = 1.009145d0 !..17f(e+nu)17o
1785 : real(dp), parameter :: Qnu_f18 = 0.393075d0 !..18f(e+nu)18o
1786 : real(dp), parameter :: Qnu_o14 = 2.22d0 !..14o(e+nu)14n
1787 : real(dp), parameter :: Qnu_ne18 = 1.87d0 !..18ne(e+nu)18f
1788 : real(dp), parameter :: Qnu_ne19 = 1.25d0 !..19ne(e+nu)19f
1789 : !real(dp), parameter :: Qnu_mg21 = 6.2d0 !..mg21(e+nu)na21
1790 : real(dp), parameter :: Qnu_mg22 = 2.1d0 !..mg22(e+nu)na22
1791 :
1792 : 1 format(a40, 1pd26.16)
1793 :
1794 0 : write (*, 1) 'expected Q for 13n(e+nu)13c', Qnu_n13
1795 0 : write (*, 1) 'calculated Q for 13n(e+nu)13c', eval_neutrino_Q(in13, ic13)
1796 0 : write (*, *)
1797 0 : write (*, 1) 'expected Q for 15o(e+nu)15n', Qnu_o15
1798 0 : write (*, 1) 'calculated Q for 15o(e+nu)15n', eval_neutrino_Q(io15, in15)
1799 0 : write (*, *)
1800 0 : write (*, 1) 'expected Q for 17f(e+nu)17o', Qnu_f17
1801 0 : write (*, 1) 'calculated Q for 17f(e+nu)17o', eval_neutrino_Q(if17, io17)
1802 0 : write (*, *)
1803 0 : write (*, 1) 'expected Q for 18f(e+nu)18o', Qnu_f18
1804 0 : write (*, 1) 'calculated Q for 18f(e+nu)18o', eval_neutrino_Q(if18, io18)
1805 0 : write (*, *)
1806 0 : write (*, 1) 'expected Q for 14o(e+nu)14n', Qnu_o14
1807 0 : write (*, 1) 'calculated Q for 14o(e+nu)14n', eval_neutrino_Q(io14, in14)
1808 0 : write (*, *)
1809 0 : write (*, 1) 'expected Q for 18ne(e+nu)18f', Qnu_ne18
1810 0 : write (*, 1) 'calculated Q for 18ne(e+nu)18f', eval_neutrino_Q(ine18, if18)
1811 0 : write (*, *)
1812 0 : write (*, 1) 'expected Q for 19ne(e+nu)19f', Qnu_ne19
1813 0 : write (*, 1) 'calculated Q for 19ne(e+nu)19f', eval_neutrino_Q(ine19, if19)
1814 0 : write (*, *)
1815 0 : write (*, 1) 'expected Q for mg22(e+nu)na22', Qnu_mg22
1816 0 : write (*, 1) 'calculated Q for mg22(e+nu)na22', eval_neutrino_Q(img22, ina22)
1817 0 : write (*, *)
1818 0 : stop
1819 :
1820 : end subroutine test_neutrino_Q
1821 :
1822 : end module test_net_support
|
