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_do_one
21 :
22 : use chem_def
23 : use chem_lib
24 : use net_def
25 : use net_lib
26 : use const_def
27 : use rates_def
28 : use utils_lib, only: mesa_error
29 :
30 : implicit none
31 :
32 : logical, parameter :: extended_set = .false.
33 : logical, parameter :: sorted = .true.
34 :
35 : logical :: qt
36 : character(len=64) :: net_file
37 : integer :: species
38 : real(dp) :: z, abar, zbar, z2bar, z53bar, ye, &
39 : eta, d_eta_dlnT, d_eta_dlnRho, eps_neu_total
40 : integer :: screening_mode
41 : real(dp) :: test_logT, test_logRho
42 : integer, pointer :: reaction_id(:)
43 : real(dp), dimension(:), pointer :: &
44 : xin, xin_copy, d_eps_nuc_dx, dxdt, d_dxdt_dRho, d_dxdt_dT
45 : real(dp), pointer :: d_dxdt_dx(:, :)
46 :
47 : contains
48 :
49 14 : subroutine do1_net(handle, symbolic)
50 : use chem_lib, only: composition_info
51 : integer, intent(in) :: handle
52 : logical, intent(in) :: symbolic
53 :
54 14 : real(dp) :: logRho, logT, Rho, T, sum, mass_correction, weak_rate_factor, &
55 28 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, xh, xhe, rate_limit
56 :
57 28 : integer :: info, i, j, k, chem_id(species), num_reactions
58 710 : real(dp), dimension(species) :: dabar_dx, dzbar_dx, dmc_dx
59 14 : real(dp), dimension(:), pointer :: eps_nuc_categories
60 14 : real(dp), pointer :: rate_factors(:)
61 : type(Net_General_Info), pointer :: g
62 14 : type(Net_Info) :: n
63 : logical :: skip_jacobian
64 :
65 14 : info = 0
66 14 : call get_net_ptr(handle, g, info)
67 14 : if (info /= 0) call mesa_error(__FILE__, __LINE__)
68 :
69 14 : num_reactions = g%num_reactions
70 :
71 14 : logRho = test_logRho
72 14 : logT = test_logT
73 :
74 14 : if (.not. qt) write (*, *)
75 :
76 : allocate ( &
77 : rate_factors(num_reactions), &
78 : eps_nuc_categories(num_categories), &
79 14 : stat=info)
80 14 : if (info /= 0) call mesa_error(__FILE__, __LINE__)
81 :
82 14 : call get_chem_id_table(handle, species, chem_id, info)
83 14 : if (info /= 0) call mesa_error(__FILE__, __LINE__)
84 :
85 : call composition_info( &
86 0 : species, chem_id, xin, xh, xhe, z, abar, zbar, z2bar, z53bar, ye, &
87 14 : mass_correction, sum, dabar_dx, dzbar_dx, dmc_dx)
88 :
89 14 : Rho = exp10(logRho)
90 14 : T = exp10(logT)
91 :
92 14 : if (net_file == '19_to_ni56.net') then
93 0 : logT = 9D+00
94 0 : logRho = 8D+00
95 0 : eta = 3D+00
96 : end if
97 :
98 14 : if (net_file == 'approx21_cr60_plus_co56.net') then
99 2 : logT = 4.6233007922659333D+00
100 2 : logRho = -1.0746410107891649D+01
101 2 : eta = -2.2590260158215202D+01
102 : end if
103 :
104 922 : rate_factors(:) = 1
105 14 : weak_rate_factor = 1
106 14 : rate_limit = 0d0
107 14 : skip_jacobian = .false.
108 14 : d_eta_dlnT = 0d0
109 14 : d_eta_dlnRho = 0d0
110 :
111 14 : if (symbolic) then
112 : call net_get_symbolic_d_dxdt_dx(handle, n, species, num_reactions, &
113 : xin, T, logT, Rho, logRho, &
114 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
115 : rate_factors, weak_rate_factor, &
116 : std_reaction_Qs, std_reaction_neuQs, &
117 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
118 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
119 : screening_mode, &
120 : eps_nuc_categories, eps_neu_total, &
121 0 : info)
122 : else
123 : call net_get(handle, skip_jacobian, n, species, num_reactions, &
124 : xin, T, logT, Rho, logRho, &
125 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
126 : rate_factors, weak_rate_factor, &
127 : std_reaction_Qs, std_reaction_neuQs, &
128 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
129 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
130 : screening_mode, &
131 : eps_nuc_categories, eps_neu_total, &
132 14 : info)
133 : end if
134 14 : if (info /= 0) then
135 0 : write (*, *) 'bad return from net_get'
136 0 : call mesa_error(__FILE__, __LINE__)
137 : end if
138 :
139 14 : if (symbolic .and. .not. qt) then
140 0 : write (*, *) 'nonzero d_dxdt_dx entries'
141 0 : k = 0
142 0 : do j = 1, species
143 0 : do i = 1, species
144 0 : if (d_dxdt_dx(i, j) /= 0) then
145 0 : k = k + 1
146 : write (*, '(a50,2i5)') &
147 0 : trim(chem_isos%name(chem_id(i)))// &
148 0 : ' '//trim(chem_isos%name(chem_id(j))), i, j
149 : end if
150 : end do
151 : end do
152 0 : write (*, '(A)')
153 0 : write (*, '(a50,i5)') 'num non zeros', k
154 0 : write (*, '(A)')
155 14 : else if (.not. qt) then
156 : call show_results( &
157 : g, n, logT, logRho, species, num_reactions, xin, &
158 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
159 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
160 7 : n%eps_nuc_categories, extended_set, sorted)
161 : end if
162 :
163 14 : deallocate (rate_factors, eps_nuc_categories)
164 :
165 14 : return
166 :
167 : write (*, *)
168 : 1 format(a40, 6x, e25.10)
169 : 2 format(a40, a6, e25.10)
170 : write (*, 1) 'abar', abar
171 : do i = 1, species
172 : write (*, 2) 'dabar_dx', trim(chem_isos%name(chem_id(i))), dabar_dx(i)
173 : end do
174 : write (*, *)
175 : write (*, 1) 'zbar', zbar
176 : do i = 1, species
177 : write (*, 2) 'dzbar_dx', trim(chem_isos%name(chem_id(i))), dzbar_dx(i)
178 : end do
179 : write (*, *)
180 :
181 56 : end subroutine do1_net
182 :
183 7 : subroutine show_results( &
184 7 : g, n, logT, logRho, species, num_reactions, xin, &
185 7 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
186 7 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
187 : eps_nuc_categories, &
188 : extended_set, sorted)
189 : type(Net_General_Info), pointer :: g
190 : type(net_info) :: n
191 : real(dp), intent(in) :: logT, logRho
192 : integer, intent(in) :: species, num_reactions
193 : real(dp), intent(in) :: xin(species)
194 : real(dp), intent(in) :: eps_nuc
195 : real(dp), intent(in) :: d_eps_nuc_dT
196 : real(dp), intent(in) :: d_eps_nuc_dRho
197 : real(dp), intent(in) :: d_eps_nuc_dx(species)
198 : real(dp), intent(in) :: dxdt(species)
199 : real(dp), intent(in) :: d_dxdt_dRho(species)
200 : real(dp), intent(in) :: d_dxdt_dT(species)
201 : real(dp), intent(in) :: d_dxdt_dx(species, species)
202 : real(dp), intent(in), dimension(num_categories) :: eps_nuc_categories
203 : logical, intent(in) :: extended_set
204 : logical, intent(in) :: sorted
205 :
206 : integer :: j
207 :
208 : include 'formats'
209 :
210 7 : if (net_file == 'approx21_cr60_plus_co56.net') then
211 1 : write (*, *)
212 1 : write (*, '(a40, f20.9)') 'log temp', logT
213 1 : write (*, '(a40, f20.9)') 'log rho', logRho
214 1 : write (*, *)
215 1 : j = irco56ec_to_fe56
216 1 : write (*, 1) 'rate_raw rco56ec_to_fe56', &
217 2 : n%rate_raw(g%net_reaction(j))
218 1 : j = irni56ec_to_co56
219 1 : write (*, 1) 'rate_raw rni56ec_to_co56', &
220 2 : n%rate_raw(g%net_reaction(j))
221 1 : return
222 : end if
223 :
224 6 : write (*, *)
225 : call show_summary_results(logT, logRho, &
226 6 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx)
227 :
228 6 : if (extended_set) then
229 0 : write (*, *)
230 : call show_all_rates( &
231 : g, num_reactions, n, &
232 0 : logT, logRho, sorted)
233 : end if
234 :
235 6 : write (*, *)
236 : call show_by_category( &
237 : g, num_reactions, eps_nuc_categories, &
238 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
239 6 : sorted)
240 :
241 6 : if (.not. extended_set) return
242 :
243 0 : write (*, *)
244 0 : call show_dx_dt(g, species, xin, dxdt, sorted)
245 :
246 0 : write (*, *)
247 0 : call show_d_eps_nuc_dx(g, species, xin, d_eps_nuc_dx, sorted)
248 :
249 0 : write (*, *)
250 :
251 14 : end subroutine show_results
252 :
253 6 : subroutine show_summary_results(logT, logRho, &
254 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx)
255 : real(dp), intent(in) :: logT, logRho
256 : real(dp), intent(in) :: eps_nuc
257 : real(dp), intent(in) :: d_eps_nuc_dT
258 : real(dp), intent(in) :: d_eps_nuc_dRho
259 : real(dp), intent(in) :: d_eps_nuc_dx(species)
260 :
261 6 : real(dp) :: T, Rho, eps, d_eps_dt, d_eps_dd
262 6 : T = exp10(logT); Rho = exp10(logRho)
263 :
264 6 : write (*, *)
265 6 : write (*, '(a40, f20.9)') 'log temp', logT
266 6 : write (*, '(a40, f20.9)') 'log rho', logRho
267 6 : eps = eps_nuc
268 6 : d_eps_dt = d_eps_nuc_dT
269 6 : d_eps_dd = d_eps_nuc_dRho
270 6 : write (*, *)
271 6 : write (*, '(a40, f20.9)') 'log(eps_nuc)', safe_log10(eps_nuc)
272 6 : write (*, *)
273 6 : write (*, '(a40, e20.9)') 'eps_nuc', eps_nuc
274 6 : write (*, *)
275 6 : write (*, '(a40, f20.9)') 'd_lneps_dlnT', d_eps_dt*T/eps
276 6 : write (*, '(a40, f20.9)') 'd_lneps_dlnRho', d_eps_dd*Rho/eps
277 :
278 6 : end subroutine show_summary_results
279 :
280 0 : subroutine show_all_rates( &
281 : g, num_reactions, n, logT, logRho, sorted)
282 : type(Net_General_Info), pointer :: g
283 : type(net_info) :: n
284 : integer, intent(in) :: num_reactions
285 : real(dp), intent(in) :: logT, logRho
286 : logical, intent(in) :: sorted
287 0 : real(dp), dimension(num_reactions) :: rfact
288 :
289 : integer :: i
290 0 : real(dp) :: T, Rho
291 0 : T = exp10(logT); Rho = exp10(logRho)
292 :
293 0 : write (*, *)
294 0 : write (*, *) 'summary of log raw rates'
295 0 : write (*, *)
296 0 : call show_log_rates(g, n%rate_raw, T, Rho, sorted)
297 0 : write (*, *)
298 0 : write (*, *) 'summary of screening factors'
299 0 : write (*, *)
300 0 : do i = 1, num_reactions
301 0 : if (n%rate_raw(i) > 1d-50) then
302 0 : rfact(i) = n%rate_screened(i)/n%rate_raw(i)
303 : else
304 0 : rfact(i) = 1
305 : end if
306 : end do
307 0 : call show_rates(g, rfact, T, Rho, sorted)
308 0 : write (*, *)
309 0 : write (*, *) 'summary of log screened rates (reactions/gm/sec)'
310 0 : write (*, *)
311 0 : call show_log_rates(g, rfact, T, Rho, sorted)
312 0 : write (*, *)
313 :
314 0 : end subroutine show_all_rates
315 :
316 6 : subroutine show_by_category( &
317 : g, num_reactions, eps_nuc_categories, &
318 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
319 : sorted)
320 : type(Net_General_Info), pointer :: g
321 : integer, intent(in) :: num_reactions
322 : real(dp), intent(in), dimension(num_categories) :: eps_nuc_categories
323 : real(dp), intent(in) :: &
324 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx(species)
325 : logical, intent(in) :: sorted
326 :
327 6 : real(dp) :: mx
328 : integer :: k, j, jmx
329 12 : logical :: flgs(rates_reaction_id_max)
330 :
331 6 : write (*, *)
332 6 : write (*, *) 'energy generation by category'
333 6 : write (*, *)
334 6 : write (*, '(a40, 3x, a20)') 'category', 'log rate (erg/g/sec)'
335 6 : write (*, *)
336 1935 : flgs = .false.
337 51 : do k = 1, num_categories
338 51 : if (.not. sorted) then
339 0 : jmx = k
340 : else
341 51 : mx = -99d99; jmx = -1
342 1275 : do j = 1, num_categories
343 1275 : if ((.not. flgs(j)) .and. eps_nuc_categories(j) > mx) then
344 186 : mx = eps_nuc_categories(j); jmx = j
345 : end if
346 : end do
347 51 : if (jmx <= 0) exit
348 51 : if (mx < 1) exit ! FOR TEST OUTPUT
349 45 : flgs(jmx) = .true.
350 : end if
351 : write (*, '(a40, 2x, f15.6, e15.6)') &
352 90 : trim(category_name(jmx)), safe_log10(eps_nuc_categories(jmx)), &
353 135 : eps_nuc_categories(jmx)
354 : end do
355 6 : write (*, *)
356 : !write(*, '(a40, 2x, f15.6, e15.6)') &
357 : ! 'log10(-photodisintegration)', safe_log10(-eps_nuc_categories(iphoto)), &
358 : ! -eps_nuc_categories(iphoto)
359 :
360 6 : write (*, *)
361 :
362 6 : end subroutine show_by_category
363 :
364 0 : subroutine show_rates(g, rts, T, Rho, sorted)
365 : type(Net_General_Info), pointer :: g
366 : real(dp), intent(in) :: rts(rates_reaction_id_max), T, Rho
367 : logical, intent(in) :: sorted
368 :
369 0 : logical :: flgs(rates_reaction_id_max)
370 0 : real(dp) :: mx
371 : integer :: k, j, jmx
372 :
373 0 : flgs = .false.
374 :
375 0 : do k = 1, g%num_reactions
376 0 : if (.not. sorted) then
377 0 : jmx = k; mx = rts(jmx)
378 : else
379 0 : mx = -99d99; jmx = -1
380 0 : do j = 1, g%num_reactions
381 0 : if ((.not. flgs(j)) .and. rts(j) > mx) then
382 0 : mx = rts(j); jmx = j
383 : end if
384 : end do
385 0 : if (jmx <= 0) exit
386 0 : if (mx < 1d-60) exit
387 0 : flgs(jmx) = .true.
388 : end if
389 0 : if (mx == 1) cycle
390 0 : write (*, '(a40, e20.9, 2e17.6)') trim(reaction_name(reaction_id(jmx))), mx
391 : end do
392 :
393 0 : end subroutine show_rates
394 :
395 0 : subroutine show_log_rates(g, rts, T, Rho, sorted)
396 : type(Net_General_Info), pointer :: g
397 : real(dp), intent(in) :: rts(:), T, Rho
398 : logical, intent(in) :: sorted
399 :
400 0 : logical :: flgs(rates_reaction_id_max)
401 0 : real(dp) :: mx
402 : integer :: k, j, jmx
403 :
404 0 : flgs = .false.
405 :
406 0 : do k = 1, g%num_reactions
407 0 : if (.not. sorted) then
408 0 : jmx = k; mx = rts(jmx)
409 : else
410 0 : mx = -99d99; jmx = -1
411 0 : do j = 1, g%num_reactions
412 0 : if ((.not. flgs(j)) .and. rts(j) > mx) then
413 0 : mx = rts(j); jmx = j
414 : end if
415 : end do
416 0 : if (jmx <= 0) exit
417 0 : if (mx < 1d-40) exit
418 0 : flgs(jmx) = .true.
419 : end if
420 0 : if (mx == 1) cycle
421 0 : write (*, '(a40, f20.9, 2e17.6)') trim(reaction_name(reaction_id(jmx))), safe_log10(mx)
422 : end do
423 :
424 0 : end subroutine show_log_rates
425 :
426 0 : subroutine show_dx_dt(g, species, xin, dxdt, sorted)
427 : type(Net_General_Info), pointer :: g
428 : integer, intent(in) :: species
429 : real(dp), intent(in) :: xin(species)
430 : real(dp), intent(in) :: dxdt(species)
431 : logical, intent(in) :: sorted
432 :
433 0 : write (*, *)
434 0 : write (*, *) 'summary of isotope mass abundance changes'
435 0 : write (*, *)
436 0 : write (*, '(a40, 2(a17))') 'isotope', 'x initial', 'dx_dt '
437 0 : call show_partials(g, species, xin, dxdt, .true., sorted)
438 :
439 0 : end subroutine show_dx_dt
440 :
441 0 : subroutine show_d_eps_nuc_dx(g, species, xin, d_eps_nuc_dx, sorted)
442 : type(Net_General_Info), pointer :: g
443 : integer, intent(in) :: species
444 : real(dp), intent(in) :: xin(species)
445 : real(dp), intent(in) :: d_eps_nuc_dx(species)
446 : logical, intent(in) :: sorted
447 :
448 0 : write (*, *)
449 0 : write (*, *) 'summary of d_eps_nuc_dx'
450 0 : write (*, *)
451 0 : write (*, '(a40, a17)') 'isotope', 'd_eps_nuc_dx'
452 0 : call show_partials(g, species, xin, d_eps_nuc_dx, .false., sorted)
453 :
454 0 : end subroutine show_d_eps_nuc_dx
455 :
456 0 : subroutine show_partials(g, species, xin, derivs, initX_flag, sorted)
457 : type(Net_General_Info), pointer :: g
458 : integer, intent(in) :: species
459 : real(dp), intent(in) :: xin(species)
460 : real(dp), intent(in) :: derivs(species)
461 : logical, intent(in) :: initX_flag, sorted
462 :
463 0 : real(dp) :: mx
464 : integer :: k, j, jmx
465 0 : integer, pointer :: chem_id(:)
466 0 : logical :: iflgs(species)
467 0 : chem_id => g%chem_id
468 0 : write (*, *)
469 0 : iflgs = .false.
470 0 : do k = 1, species
471 0 : if (.not. sorted) then
472 0 : jmx = k
473 : else
474 0 : mx = -99d99; jmx = -1
475 0 : do j = 1, species
476 0 : if ((.not. iflgs(j)) .and. abs(derivs(j)) > mx) then
477 0 : mx = abs(derivs(j)); jmx = j
478 : end if
479 : end do
480 0 : if (jmx <= 0) exit
481 0 : if (mx < 1d-40) exit
482 : end if
483 0 : if (initX_flag) then
484 0 : write (*, '(a40, 2e17.6)') trim(chem_isos%name(chem_id(jmx))), xin(jmx), derivs(jmx)
485 : else
486 0 : write (*, '(a40, e25.14)') trim(chem_isos%name(chem_id(jmx))), derivs(jmx)
487 : end if
488 0 : iflgs(jmx) = .true.
489 : end do
490 :
491 0 : end subroutine show_partials
492 :
493 : end module test_net_do_one
|
