Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 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 net_eval
21 :
22 : use const_def, only: dp, i8, Qconv, arg_not_provided
23 : use math_lib
24 : use chem_def
25 : use chem_lib, only: get_mass_excess
26 : use net_def, only: Net_General_Info, Net_Info
27 : use utils_lib, only: fill_with_NaNs
28 :
29 : implicit none
30 :
31 : contains
32 :
33 89108 : subroutine eval_net( &
34 : n, g, rates_only, just_dxdt, &
35 : num_isos, num_reactions, num_wk_reactions, &
36 89108 : x, temp, logtemp, rho, logrho, &
37 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
38 89108 : rate_factors, weak_rate_factor, &
39 : reaction_Qs, reaction_neuQs, &
40 89108 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
41 89108 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
42 : screening_mode, &
43 89108 : eps_nuc_categories, eps_neu_total, &
44 : actual_Qs, actual_neuQs, from_weaklib, symbolic, ierr)
45 : use net_initialize, only: &
46 : setup_net_info, set_ptrs_for_approx21
47 : use net_approx21, only: num_reactions_func => num_reactions
48 : use net_screen
49 : use net_derivs
50 : use net_def, only: &
51 : net_test_partials, &
52 : net_test_partials_val, net_test_partials_dval_dx, &
53 : net_test_partials_i, net_test_partials_iother
54 :
55 : type (Net_Info) :: n
56 : type (Net_General_Info), pointer :: g
57 : logical, intent(in) :: rates_only, just_dxdt
58 : integer, intent(in) :: num_isos
59 : integer, intent(in) :: num_reactions, num_wk_reactions
60 : real(dp), intent(in) :: x(:)
61 : real(dp), intent(in) :: temp, logtemp
62 : real(dp), intent(in) :: rho, logrho
63 : real(dp), intent(in) :: abar ! mean number of nucleons per nucleus
64 : real(dp), intent(in) :: zbar ! mean charge per nucleus
65 : real(dp), intent(in) :: z2bar ! mean charge squared per nucleus
66 : real(dp), intent(in) :: ye
67 : real(dp), intent(in) :: eta, d_eta_dlnT, d_eta_dlnRho ! eta and derivatives
68 : real(dp), intent(in) :: rate_factors(:)
69 : real(dp), intent(in) :: weak_rate_factor
70 : real(dp), pointer, intent(in) :: reaction_Qs(:) ! (rates_reaction_id_max)
71 : real(dp), pointer, intent(in) :: reaction_neuQs(:) ! (rates_reaction_id_max)
72 : real(dp), intent(out) :: eps_nuc ! ergs/gram/second from burning
73 : real(dp), intent(out) :: d_eps_nuc_dT
74 : real(dp), intent(out) :: d_eps_nuc_dRho
75 : real(dp), intent(inout) :: d_eps_nuc_dx(:)
76 : real(dp), intent(inout) :: dxdt(:)
77 : real(dp), intent(inout) :: d_dxdt_dRho(:)
78 : real(dp), intent(inout) :: d_dxdt_dT(:)
79 : real(dp), intent(inout) :: d_dxdt_dx(:,:)
80 : real(dp), intent(inout) :: eps_nuc_categories(:)
81 : real(dp), intent(out) :: eps_neu_total
82 : integer, intent(in) :: screening_mode
83 : real(dp), pointer, dimension(:) :: actual_Qs, actual_neuQs ! ignore if null
84 : logical, pointer :: from_weaklib(:) ! ignore if null
85 : logical, intent(in) :: symbolic
86 : integer, intent(out) :: ierr
87 :
88 : integer, parameter :: max_z_for_cache = 14
89 89108 : real(dp) :: enuc, T9, total, prev, curr, prev_T
90 : real(dp) :: eps_total, Ys, sum_dxdt, compare, Z_plus_N
91 356432 : real(qp) :: eps_nuc_MeV(num_rvs)
92 : integer :: ci, i, j, ir, weak_id, h1, iwork
93 : integer(i8) :: time0, time1
94 : logical :: doing_timing
95 :
96 : logical, parameter :: dbg = .false.
97 : !logical, parameter :: dbg = .true.
98 :
99 : include 'formats'
100 :
101 : if (dbg) write(*,*) 'enter eval_net'
102 :
103 89108 : doing_timing = g% doing_timing
104 89108 : if (doing_timing) then
105 0 : call system_clock(time0)
106 0 : g% doing_timing = .false.
107 : end if
108 :
109 89108 : if (.not. g% net_has_been_defined) then
110 0 : ierr = -1
111 : if (dbg) write(*,*) 'failed (.not. g% net_has_been_defined)'
112 0 : return
113 : end if
114 :
115 : if (temp == arg_not_provided .or. logtemp == arg_not_provided .or. &
116 89108 : rho == arg_not_provided .or. logrho == arg_not_provided) then
117 0 : write(*,*) "You must now eplxicity pass both the linear and log values of the temperature and density"
118 0 : ierr = -1
119 0 : return
120 : end if
121 :
122 89108 : ierr = 0
123 89108 : n% g => g
124 :
125 : if (dbg) write(*,*) 'call setup_net_info'
126 89108 : call setup_net_info(n)
127 :
128 89108 : n% reaction_Qs => reaction_Qs
129 89108 : n% reaction_neuQs => reaction_neuQs
130 89108 : n% weak_rate_factor = weak_rate_factor
131 89108 : n% logT = logtemp
132 89108 : n% temp = temp
133 89108 : n% logRho = logrho
134 89108 : n% rho = rho
135 89108 : n% screening_mode = screening_mode
136 1009500 : n% x = x
137 89108 : n% zbar = zbar
138 89108 : n% abar = abar
139 89108 : n% z2bar = z2bar
140 89108 : n% ye = ye
141 89108 : n% eta = eta
142 89108 : n% d_eta_dlnT = d_eta_dlnT
143 89108 : n% d_eta_dlnRho = d_eta_dlnRho
144 26383522 : n% rate_factors = rate_factors
145 :
146 89108 : if (n% logT < rattab_tlo) then ! clip to table so can eval beta decays
147 18826 : n% logT = rattab_tlo
148 18826 : n% temp = rattab_temp_lo
149 : end if
150 :
151 89108 : T9 = n% temp*1d-9
152 :
153 89108 : if (g% doing_approx21) then
154 9106 : call set_ptrs_for_approx21(n)
155 : end if
156 :
157 : if (dbg) write(*,*) 'call set_molar_abundances'
158 89108 : call set_molar_abundances(n, dbg, ierr)
159 89108 : if (ierr /= 0) then
160 : if (dbg) write(*,*) 'failed in set_molar_abundances'
161 : return
162 : end if
163 :
164 89108 : if (num_wk_reactions > 0) then
165 : if (dbg) write(*,*) 'call get_weaklib_rates'
166 89108 : call get_weaklib_rates(n, ierr)
167 89108 : if (ierr /= 0) then
168 : if (dbg) write(*,*) 'failed in get_weaklib_rates'
169 : return
170 : end if
171 : end if
172 :
173 89108 : if (associated(actual_Qs) .and. associated(actual_neuQs)) then
174 84224 : do i = 1, g% num_reactions
175 83328 : ir = g% reaction_id(i)
176 83328 : from_weaklib(i) = .false.
177 83328 : actual_Qs(i) = n% reaction_Qs(ir)
178 83328 : actual_neuQs(i) = n% reaction_neuQs(ir)
179 83328 : weak_id = g% weak_reaction_index(i)
180 84224 : if (weak_id > 0) then
181 2688 : if (g% weaklib_ids(weak_id) > 0) then
182 0 : from_weaklib(i) = .true.
183 0 : actual_Qs(i) = n% Q(weak_id)
184 0 : actual_neuQs(i) = n% Qneu(weak_id)
185 : end if
186 : end if
187 : end do
188 : end if
189 :
190 89108 : if (doing_timing) then
191 0 : call system_clock(time1)
192 0 : g% clock_net_eval = g% clock_net_eval + (time1 - time0)
193 0 : time0 = time1
194 : end if
195 :
196 : if (dbg) write(*,*) 'call get_rates_with_screening'
197 89108 : call get_rates_with_screening(n, ierr)
198 : if (dbg) write(*,*) 'done get_rates_with_screening'
199 89108 : if (ierr /= 0) then
200 : if (dbg) write(*,*) 'failed in get_rates_with_screening'
201 : return
202 : end if
203 :
204 89108 : if (rates_only) return
205 :
206 : ! n% d_eps_nuc_dT = 0
207 : ! n% d_eps_nuc_dRho = 0
208 : ! n% d_eps_nuc_dx(:) = 0
209 :
210 : ! n% dxdt(:) = 0
211 : ! n% d_dxdt_dRho(:) = 0
212 : ! n% d_dxdt_dT(:) = 0
213 : ! if (.not. just_dxdt) d_dxdt_dx(:,:) = 0
214 2227700 : n% eps_nuc_categories(:) = 0
215 89108 : n% eps_neu_total = 0
216 920392 : n% d_eps_nuc_dy = 0
217 :
218 89108 : if (g% doing_approx21) then
219 9106 : call eval_net_approx21_procs(n, just_dxdt, ierr)
220 9106 : if (ierr /= 0) return
221 :
222 9106 : if (net_test_partials) then
223 0 : net_test_partials_val = eps_nuc
224 0 : net_test_partials_dval_dx = d_eps_nuc_dx(net_test_partials_i)
225 0 : if (g% add_co56_to_approx21) then
226 0 : write(*,*) 'net: eval_net_approx21_plus_co56'
227 : else
228 0 : write(*,*) 'net: eval_net_approx21'
229 : end if
230 : end if
231 :
232 : call unpack_for_export(n, eps_nuc, d_eps_nuc_dT, d_eps_nuc_dRho, d_eps_nuc_dx, &
233 : eps_neu_total, &
234 : dxdt, d_dxdt_dT, d_dxdt_dRho, d_dxdt_dx, &
235 9106 : eps_nuc_categories)
236 :
237 9106 : return
238 : end if ! End of approx21
239 :
240 : if (dbg) write(*,*) 'call get_derivs'
241 : call get_derivs( &
242 : n, n% dydt, eps_nuc_MeV(1:num_rvs), n% eta, n% ye, &
243 : n% logT, n% temp, n% rho, n% abar, n% zbar, &
244 : num_reactions, n% rate_factors, &
245 80002 : symbolic, just_dxdt, ierr)
246 80002 : if (ierr /= 0) then
247 : if (dbg) write(*,*) 'failed in get_derivs'
248 : return
249 : end if
250 :
251 80002 : if (symbolic) then
252 0 : do j=1, num_isos
253 0 : do i=1, num_isos
254 0 : d_dxdt_dx(i,j) = n% d_dydt_dy(i,j)
255 : end do
256 : end do
257 : return
258 : end if
259 :
260 80002 : if (doing_timing) then
261 0 : call system_clock(time1)
262 0 : g% clock_net_derivs = g% clock_net_derivs + (time1 - time0)
263 0 : time0 = time1
264 : end if
265 :
266 : ! convert the eps_nuc_categories
267 2000050 : do i=1,num_categories
268 2000050 : n% eps_nuc_categories(i) = Qconv * n% eps_nuc_categories(i)
269 : end do
270 :
271 : ! store the results
272 720056 : do i=1,num_isos
273 640054 : ci = g% chem_id(i)
274 720056 : n% dxdt(i) = chem_isos% Z_plus_N(ci) * n% dydt(i_rate, i)
275 : end do
276 :
277 80002 : if (.not. just_dxdt) call store_partials(n)
278 :
279 80002 : n% eps_nuc = eps_nuc_MeV(i_rate)*Qconv
280 80002 : n% d_eps_nuc_dT = eps_nuc_MeV(i_rate_dT)*Qconv
281 80002 : n% d_eps_nuc_dRho = eps_nuc_MeV(i_rate_dRho)*Qconv
282 :
283 80002 : n% eps_neu_total = n% eps_neu_total * Qconv
284 :
285 80002 : if (doing_timing) then
286 0 : call system_clock(time1)
287 0 : g% clock_net_eval = g% clock_net_eval + (time1 - time0)
288 0 : g% doing_timing = .true.
289 : end if
290 :
291 80002 : if (net_test_partials) then
292 : !net_test_partials_val = eps_nuc
293 : !net_test_partials_dval_dx = d_eps_nuc_dx(net_test_partials_i)
294 : net_test_partials_val = &
295 : n% rate_screened(g% net_reaction(irn14_to_c12))/ &
296 0 : n% rate_raw(g% net_reaction(irn14_to_c12))
297 0 : net_test_partials_dval_dx = 0d0
298 0 : write(*,*) 'net_test_partials'
299 : end if
300 :
301 : call unpack_for_export(n, eps_nuc, d_eps_nuc_dT, d_eps_nuc_dRho, d_eps_nuc_dx, &
302 : eps_neu_total, &
303 : dxdt, d_dxdt_dT, d_dxdt_dRho, d_dxdt_dx,&
304 80002 : eps_nuc_categories)
305 :
306 89108 : end subroutine eval_net
307 :
308 89108 : subroutine unpack_for_export(n, eps_nuc, d_eps_nuc_dT, d_eps_nuc_dRho, d_eps_nuc_dx, &
309 : eps_neu_total, &
310 89108 : dxdt, d_dxdt_dT, d_dxdt_dRho, d_dxdt_dx,&
311 89108 : eps_nuc_categories)
312 : type(net_info) :: n
313 : real(dp), intent(out) :: eps_nuc ! ergs/gram/second from burning
314 : real(dp), intent(out) :: d_eps_nuc_dT
315 : real(dp), intent(out) :: d_eps_nuc_dRho
316 : real(dp), intent(out) :: d_eps_nuc_dx(:)
317 : real(dp), intent(out) :: dxdt(:)
318 : real(dp), intent(out) :: d_dxdt_dRho(:)
319 : real(dp), intent(out) :: d_dxdt_dT(:)
320 : real(dp), intent(out) :: d_dxdt_dx(:,:)
321 : real(dp), intent(out) :: eps_neu_total
322 : real(dp), intent(out) :: eps_nuc_categories(:)
323 :
324 89108 : eps_nuc = n% eps_nuc
325 89108 : d_eps_nuc_dT = n% d_eps_nuc_dT
326 89108 : d_eps_nuc_dRho = n% d_eps_nuc_dRho
327 920392 : d_eps_nuc_dx = n% d_eps_nuc_dx
328 :
329 920392 : dxdt = n% dxdt
330 920392 : d_dxdt_dT = n% d_dxdt_dT
331 920392 : d_dxdt_dRho = n% d_dxdt_dRho
332 10057424 : d_dxdt_dx = n% d_dxdt_dx
333 :
334 89108 : eps_neu_total = n% eps_neu_total
335 :
336 2227700 : eps_nuc_categories = n% eps_nuc_categories
337 :
338 89108 : end subroutine unpack_for_export
339 :
340 :
341 9106 : subroutine eval_net_approx21_procs(n,just_dxdt, ierr)
342 : use net_approx21
343 : use rates_def
344 : type(net_info) :: n
345 : type(net_general_info), pointer :: g=>null()
346 : logical,intent(in) :: just_dxdt
347 : integer :: ierr
348 :
349 : integer :: ci, i, j, num_isos
350 9106 : real(dp) :: Z_plus_N
351 :
352 9106 : ierr = 0
353 :
354 9106 : g => n% g
355 :
356 9106 : num_isos = g% num_isos
357 :
358 : call approx21_special_reactions( &
359 : n% temp, n% rho, n% abar, n% zbar, n% y, &
360 : g% use_3a_fl87, Qconv * n% reaction_Qs(ir_he4_he4_he4_to_c12), &
361 : n% rate_screened, n% rate_screened_dT, n% rate_screened_dRho, &
362 9106 : n% dratdumdy1, n% dratdumdy2, g% add_co56_to_approx21, ierr)
363 9106 : if (ierr /= 0) return
364 :
365 : call approx21_dydt( &
366 : n% y, n% rate_screened, n% rate_screened, &
367 : n% dydt1, .false., g% fe56ec_fake_factor, g% min_T_for_fe56ec_fake_factor, &
368 9106 : g% fe56ec_n_neut, n% temp, n% rho, g% add_co56_to_approx21, ierr)
369 9106 : if (ierr /= 0) return
370 :
371 9106 : n% fII = approx21_eval_PPII_fraction(n% y, n% rate_screened)
372 :
373 : call get_approx21_eps_info( n, &
374 : n% dydt1, n% rate_screened, .true., n% eps_total, n% eps_neu_total, &
375 9106 : g% add_co56_to_approx21, ierr)
376 :
377 9106 : if (ierr /= 0) return
378 9106 : n% eps_nuc = n% eps_total - n% eps_neu_total
379 :
380 200336 : do i=1, num_isos
381 200336 : n% dxdt(i) = chem_isos% Z_plus_N(g% chem_id(i)) * n% dydt1(i)
382 : end do
383 :
384 9106 : if (just_dxdt) return
385 :
386 : call approx21_dfdy( &
387 : n% y, n% dfdy, &
388 : g% fe56ec_fake_factor, g% min_T_for_fe56ec_fake_factor, g% fe56ec_n_neut, &
389 : n% rate_screened, n% rate_screened_dT, n% rate_screened_dRho, &
390 6044 : n% dratdumdy1, n% dratdumdy2, n% temp, g% add_co56_to_approx21, ierr)
391 6044 : if (ierr /= 0) return
392 :
393 : call approx21_dfdT_dfdRho( &
394 :
395 : ! NOTE: currently this gives d_eps_total_dy -- should fix to account for neutrinos too
396 :
397 : n% y, g% mion, n% dfdy, n% rate_screened, n% rate_screened_dT, n% rate_screened_dRho, &
398 : g% fe56ec_fake_factor, g% min_T_for_fe56ec_fake_factor, &
399 6044 : g% fe56ec_n_neut, n% temp, n% rho, n% dfdT, n% dfdRho, n% d_epsnuc_dy, g% add_co56_to_approx21, ierr)
400 6044 : if (ierr /= 0) return
401 :
402 : call get_approx21_eps_info( n, &
403 : n% dfdT, n% rate_screened_dT, .false., n% deps_total_dT, n% deps_neu_dT, &
404 6044 : g% add_co56_to_approx21, ierr)
405 :
406 6044 : if (ierr /= 0) return
407 6044 : n% d_eps_nuc_dT = n% deps_total_dT - n% deps_neu_dT
408 :
409 : call get_approx21_eps_info( n, &
410 : n% dfdRho, n% rate_screened_dRho, .false., n% deps_total_dRho, n% deps_neu_dRho, &
411 6044 : g% add_co56_to_approx21, ierr)
412 :
413 6044 : if (ierr /= 0) return
414 6044 : n% d_eps_nuc_dRho = n% deps_total_dRho - n% deps_neu_dRho
415 :
416 : call approx21_d_epsneu_dy( &
417 : n% y, n% rate_screened, &
418 : n% reaction_neuQs(irpp_to_he3), &
419 : n% reaction_neuQs(ir34_pp2), &
420 : n% reaction_neuQs(ir34_pp3), &
421 : n% reaction_neuQs(irc12_to_n14), &
422 : n% reaction_neuQs(irn14_to_c12), &
423 : n% reaction_neuQs(iro16_to_n14), &
424 : n% d_epsneu_dy, &
425 6044 : g% add_co56_to_approx21, ierr)
426 6044 : if (ierr /= 0) return
427 :
428 132972 : do i=1, n% g%num_isos
429 126928 : ci = g% chem_id(i)
430 126928 : Z_plus_N = dble(chem_isos% Z_plus_N(ci))
431 126928 : n% d_eps_nuc_dx(i) = (n% d_epsnuc_dy(i) - n% d_epsneu_dy(i))/Z_plus_N
432 126928 : n% d_dxdt_dRho(i) = Z_plus_N * n% dfdRho(i)
433 126928 : n% d_dxdt_dT(i) = Z_plus_N * n% dfdT(i)
434 2798548 : do j=1, num_isos
435 : n% d_dxdt_dx(i,j) = &
436 2792504 : n% dfdy(i,j)*Z_plus_N/chem_isos% Z_plus_N(g% chem_id(j))
437 : end do
438 : end do
439 :
440 9106 : end subroutine eval_net_approx21_procs
441 :
442 :
443 21194 : subroutine get_approx21_eps_info(n, &
444 21194 : dydt1, rate_screened, do_eps_nuc_categories, eps_total, eps_neu_total, plus_co56, ierr)
445 9106 : use net_approx21, only: approx21_eps_info
446 : use rates_def
447 : type(net_info) :: n
448 : type(net_general_info), pointer :: g=>null()
449 : real(dp), intent(in), dimension(:) :: dydt1, rate_screened
450 : logical, intent(in) :: do_eps_nuc_categories
451 : real(dp), intent(out) :: eps_total, eps_neu_total
452 : logical, intent(in) :: plus_co56
453 : integer, intent(out) :: ierr
454 : real(dp) :: Qtotal_rfe56ec, Qneu_rfe56ec
455 :
456 21194 : g => n% g
457 :
458 : ! Indexes into reaction_Qs and reaction_neuQs should be in terms of the
459 : ! normal rate ids not the approx21 rate ids (in net_approx21.f90)
460 :
461 21194 : call get_Qs_rfe56ec(n, Qtotal_rfe56ec, Qneu_rfe56ec)
462 :
463 : call approx21_eps_info( &
464 : n, n% y, g% mion, dydt1, rate_screened, n% fII, &
465 : n% reaction_Qs(irpp_to_he3), n% reaction_neuQs(irpp_to_he3), &
466 : n% reaction_Qs(ir_he3_he3_to_h1_h1_he4), &
467 : n% reaction_Qs(ir34_pp2), n% reaction_neuQs(ir34_pp2), &
468 : n% reaction_Qs(ir34_pp3), n% reaction_neuQs(ir34_pp3), &
469 : n% reaction_Qs(irc12_to_n14), n% reaction_neuQs(irc12_to_n14), &
470 : n% reaction_Qs(irn14_to_c12), n% reaction_neuQs(irn14_to_c12), &
471 : n% reaction_Qs(iro16_to_n14), n% reaction_neuQs(iro16_to_n14), &
472 : n% reaction_Qs(irn14_to_o16), &
473 :
474 : n% reaction_Qs(irprot_to_neut), n% reaction_neuQs(irprot_to_neut), &
475 : n% reaction_Qs(irneut_to_prot), n% reaction_neuQs(irneut_to_prot), &
476 : n% reaction_Qs(irni56ec_to_co56), n% reaction_neuQs(irni56ec_to_co56), &
477 : n% reaction_Qs(irco56ec_to_fe56), n% reaction_neuQs(irco56ec_to_fe56), &
478 : Qtotal_rfe56ec, Qneu_rfe56ec, &
479 :
480 : n% reaction_Qs(irn14ag_lite), &
481 : n% reaction_Qs(ir_he4_he4_he4_to_c12), &
482 : n% reaction_Qs(ir_c12_ag_o16), n% reaction_Qs(ir_o16_ag_ne20), &
483 : n% reaction_Qs(ir1212), &
484 : n% reaction_Qs(ir1216_to_mg24), n% reaction_Qs(ir1216_to_si28), &
485 : n% reaction_Qs(ir1616a), n% reaction_Qs(ir1616g), &
486 : n% reaction_Qs(ir_ne20_ag_mg24), &
487 : n% reaction_Qs(ir_mg24_ag_si28), &
488 : n% reaction_Qs(ir_si28_ag_s32), &
489 : n% reaction_Qs(ir_s32_ag_ar36), &
490 : n% reaction_Qs(ir_ar36_ag_ca40), &
491 : n% reaction_Qs(ir_ca40_ag_ti44), &
492 : n% reaction_Qs(ir_ti44_ag_cr48), &
493 : n% reaction_Qs(ir_cr48_ag_fe52), &
494 : n% reaction_Qs(ir_fe52_ag_ni56), &
495 : n% reaction_Qs(ir_fe52_ng_fe53), &
496 : n% reaction_Qs(ir_fe53_ng_fe54), &
497 : n% reaction_Qs(ir_fe54_ng_fe55), &
498 : n% reaction_Qs(ir_fe55_ng_fe56), &
499 : n% reaction_Qs(irfe52neut_to_fe54), &
500 : n% reaction_Qs(irfe52aprot_to_fe54), &
501 : n% reaction_Qs(irfe54ng_to_fe56), &
502 : n% reaction_Qs(irfe54aprot_to_fe56), &
503 : n% reaction_Qs(irfe52aprot_to_ni56), &
504 : n% reaction_Qs(irfe54prot_to_ni56), &
505 : n% reaction_Qs(irhe4_breakup), &
506 : n% reaction_Qs(irhe4_rebuild), &
507 : eps_total, eps_neu_total, & ! Dont use n% here as we call this for both eps_neu and eps_neu_dt and drho
508 : do_eps_nuc_categories, n% eps_nuc_categories, &
509 21194 : .false., plus_co56, ierr)
510 :
511 21194 : end subroutine get_approx21_eps_info
512 :
513 21194 : subroutine get_Qs_rfe56ec(n, Qtotal, Qneu)
514 21194 : use chem_def
515 : use rates_def
516 : type(net_info) :: n
517 : real(dp), intent(out) :: Qtotal, Qneu
518 : integer :: id, ir
519 : include 'formats'
520 21194 : id = n% g% approx21_ye_iso
521 21194 : if (id == imn56) then
522 : ir = irfe56ec_fake_to_mn56
523 21194 : else if (id == imn57) then
524 : ir = irfe56ec_fake_to_mn57
525 21194 : else if (id == icr56) then
526 : ir = irfe56ec_fake_to_cr56
527 6 : else if (id == icr57) then
528 : ir = irfe56ec_fake_to_cr57
529 6 : else if (id == icr58) then
530 : ir = irfe56ec_fake_to_cr58
531 6 : else if (id == icr59) then
532 : ir = irfe56ec_fake_to_cr59
533 6 : else if (id == icr60) then
534 : ir = irfe56ec_fake_to_cr60
535 0 : else if (id == icr61) then
536 : ir = irfe56ec_fake_to_cr61
537 0 : else if (id == icr62) then
538 : ir = irfe56ec_fake_to_cr62
539 0 : else if (id == icr63) then
540 : ir = irfe56ec_fake_to_cr63
541 0 : else if (id == icr64) then
542 : ir = irfe56ec_fake_to_cr64
543 0 : else if (id == icr65) then
544 : ir = irfe56ec_fake_to_cr65
545 0 : else if (id == icr66) then
546 : ir = irfe56ec_fake_to_cr66
547 : else
548 0 : ir = irco56ec_to_fe56
549 : end if
550 21194 : Qtotal = n% reaction_Qs(ir)
551 21194 : Qneu = n% reaction_neuQs(ir)
552 21194 : end subroutine get_Qs_rfe56ec
553 :
554 80002 : subroutine store_partials(n)
555 21194 : use rates_def, only: i_rate, i_rate_dT, i_rate_dRho
556 : type(net_info) :: n
557 : type(net_general_info), pointer :: g => null()
558 : integer :: i, j, ci
559 80002 : real(dp) :: Z_plus_N
560 : include 'formats'
561 80002 : g => n%g
562 720056 : do i=1, g% num_isos
563 640054 : ci = g% chem_id(i)
564 640054 : Z_plus_N = dble(chem_isos% Z_plus_N(ci))
565 640054 : n% d_eps_nuc_dx(i) = Qconv*n% d_eps_nuc_dy(i)/Z_plus_N
566 640054 : n% dxdt(i) = Z_plus_N * n% dydt(i_rate, i)
567 640054 : n% d_dxdt_dRho(i) = Z_plus_N * n% dydt(i_rate_dRho, i)
568 640054 : n% d_dxdt_dT(i) = Z_plus_N*n% dydt(i_rate_dT, i)
569 5841170 : do j=1, g% num_isos
570 : n% d_dxdt_dx(i,j) = &
571 5761168 : n% d_dydt_dy(i,j)*Z_plus_N/chem_isos% Z_plus_N(g% chem_id(j))
572 : end do
573 : end do
574 :
575 80002 : end subroutine store_partials
576 :
577 89108 : subroutine get_rates_with_screening(n, ierr)
578 80002 : use rates_def, only: reaction_inputs, nrattab
579 : use rates_lib, only: eval_using_rate_tables
580 : use net_approx21, only: num_reactions_func => num_reactions
581 : use net_screen
582 :
583 : type(net_info) :: n
584 : type(net_general_info),pointer :: g=> null()
585 :
586 : integer, intent(out) :: ierr
587 :
588 : logical, parameter :: dbg=.false.
589 : integer(i8) :: time0, time1
590 :
591 : integer :: i, num, num_reactions
592 : real(dp) :: f
593 : logical :: okay
594 :
595 : include 'formats'
596 :
597 89108 : g => n% g
598 :
599 :
600 3416130 : do i=1,g% num_reactions
601 3416130 : if (g% reaction_id(i) <= 0) then
602 0 : write(*,2) 'g% reaction_id(i)', i, g% reaction_id(i)
603 0 : call mesa_error(__FILE__,__LINE__,'get_rates_with_screening')
604 : end if
605 : end do
606 :
607 : if (dbg) write(*,*) 'call eval_using_rate_tables'
608 : call eval_using_rate_tables( &
609 : g% num_reactions, g% reaction_id, g% rate_table, g% rattab_f1, nrattab, &
610 : n% ye, n% logT, n% temp, n% rho, n% rate_factors, g% logttab, &
611 89108 : n% rate_raw, n% rate_raw_dT, n% rate_raw_dRho, ierr)
612 89108 : if (ierr /= 0) then
613 : if (dbg) write(*,*) 'ierr from eval_using_rate_tables'
614 : return
615 : end if
616 :
617 89108 : if (g% doing_timing) then
618 0 : call system_clock(time0)
619 : end if
620 :
621 89108 : if (g% doing_approx21) then
622 9106 : call approx21_rates(n, g% add_co56_to_approx21,ierr)
623 9106 : if (ierr /= 0) return
624 : end if
625 :
626 : ! get the reaction rates including screening factors
627 : if (dbg) write(*,*) 'call screen_net with init=.false.'
628 : call screen_net( &
629 : g, g% num_isos, n% y, n% temp, n% rho, n% logT, n% logrho, .false., &
630 : n% rate_raw, n% rate_raw_dT, n% rate_raw_dRho, &
631 : n% rate_screened, n% rate_screened_dT, n% rate_screened_dRho, &
632 : n% screening_mode, &
633 89108 : n% zbar, n% abar, n% z2bar, n% ye, ierr)
634 : if (dbg) write(*,*) 'done screen_net with init=.false.'
635 89108 : if (ierr /= 0) return
636 89108 : if (g% doing_approx21) then
637 9106 : num = num_reactions_func(g% add_co56_to_approx21)
638 218544 : do i=g% num_reactions+1,num
639 209438 : n% rate_screened(i) = n% rate_raw(i)
640 209438 : n% rate_screened_dT(i) = n% rate_raw_dT(i)
641 218544 : n% rate_screened_dRho(i) = n% rate_raw_dRho(i)
642 : end do
643 1065406 : do i=1,num
644 1056300 : n% dratdumdy1(i) = 0d0
645 1065406 : n% dratdumdy2(i) = 0d0
646 : end do
647 : end if
648 :
649 :
650 89108 : if (g% doing_timing) then
651 0 : call system_clock(time1)
652 0 : g% clock_net_screen = g% clock_net_screen + (time1 - time0)
653 : end if
654 :
655 89108 : end subroutine get_rates_with_screening
656 :
657 9106 : subroutine approx21_rates(n, plus_co56, ierr)
658 : use net_approx21, only: &
659 89108 : approx21_pa_pg_fractions, approx21_weak_rates
660 : type(net_info) :: n
661 : logical, intent(in) :: plus_co56
662 : integer, intent(out) :: ierr
663 : ierr = 0
664 : call approx21_pa_pg_fractions( &
665 9106 : n% rate_raw, n% rate_raw_dT, n% rate_raw_dRho, ierr)
666 9106 : if (ierr /= 0) return
667 : call approx21_weak_rates( &
668 : n% y, n% rate_raw, n% rate_raw_dT, n% rate_raw_dRho, &
669 : n% temp, n% rho, n% ye, n% eta, n% zbar, &
670 9106 : n% weak_rate_factor, plus_co56, ierr)
671 9106 : if (ierr /= 0) return
672 9106 : end subroutine approx21_rates
673 :
674 :
675 89108 : subroutine get_weaklib_rates(n, ierr)
676 9106 : use rates_def, only : Coulomb_Info
677 : use rates_lib, only: eval_weak_reaction_info, coulomb_set_context
678 : use net_def, only: other_kind
679 :
680 : type (net_info) :: n
681 : type(net_general_info), pointer :: g
682 :
683 : type (Coulomb_Info), target :: cc_info
684 : type (Coulomb_Info), pointer :: cc
685 :
686 : integer, intent(out) :: ierr
687 : integer :: i, j, id, ir
688 : integer(i8) :: time0, time1
689 :
690 : include 'formats'
691 :
692 : ! before getting the weaklib rates, the Coulomb_Info
693 : ! structure must be populated. the ecapture routines need
694 : ! to know some local quantities (functions of the density,
695 : ! temperature, and composition), to calculate Coulomb
696 : ! corrections to the rates
697 :
698 89108 : ierr = 0
699 89108 : cc => cc_info
700 :
701 89108 : g => n% g
702 :
703 89108 : if(g% doing_timing) then
704 0 : call system_clock(time0)
705 : end if
706 :
707 89108 : call coulomb_set_context(cc, n% temp, n% rho, n% logT, n% logRho, n% zbar, n% abar, n% z2bar)
708 :
709 : call eval_weak_reaction_info( &
710 : g% num_wk_reactions, &
711 : g% weaklib_ids(1:g% num_wk_reactions), &
712 : g% reaction_id_for_weak_reactions(1:g% num_wk_reactions), &
713 : cc, n% temp*1d-9, n% ye * n% rho, n% eta, n% d_eta_dlnT, n% d_eta_dlnRho, &
714 : n% lambda, n% dlambda_dlnT, n% dlambda_dlnRho, &
715 : n% Q, n% dQ_dlnT, n% dQ_dlnRho, &
716 : n% Qneu, n% dQneu_dlnT, n% dQneu_dlnRho, &
717 89108 : ierr)
718 89108 : if (n% weak_rate_factor < 1d0) then
719 0 : do i=1,g% num_wk_reactions
720 0 : n% lambda(i) = n% weak_rate_factor*n% lambda(i)
721 0 : n% dlambda_dlnT(i) = n% weak_rate_factor*n% dlambda_dlnT(i)
722 0 : n% dlambda_dlnRho(i) = n% weak_rate_factor*n% dlambda_dlnRho(i)
723 : end do
724 : end if
725 89108 : if (g% doing_timing) then
726 0 : call system_clock(time1)
727 0 : g% clock_net_weak_rates = g% clock_net_weak_rates + (time1 - time0)
728 : end if
729 :
730 89108 : end subroutine get_weaklib_rates
731 :
732 :
733 0 : subroutine get_T_limit_factor( &
734 : temp, lnT, T_lo, T_hi, lnT_lo, lnT_hi, &
735 : min_ln_factor, min_factor, &
736 : factor, d_factor_dT)
737 : real(dp), intent(in) :: &
738 : temp, lnT, T_lo, T_hi, lnT_lo, lnT_hi, &
739 : min_ln_factor, min_factor
740 : real(dp), intent(out) :: &
741 : factor, d_factor_dT
742 0 : real(dp) :: ln_factor, d_ln_factor_dlnT
743 0 : factor = 1d0
744 0 : d_factor_dT = 0d0
745 0 : if (temp <= T_lo) return
746 0 : if (temp >= T_hi) then
747 0 : factor = min_factor
748 0 : return
749 : end if
750 0 : ln_factor = min_ln_factor*(lnT - lnT_lo)/(lnT_hi - lnT_lo)
751 0 : d_ln_factor_dlnT = min_ln_factor/(lnT_hi - lnT_lo)
752 0 : factor = exp(ln_factor)
753 0 : d_factor_dT = d_ln_factor_dlnT*factor/temp
754 89108 : end subroutine get_T_limit_factor
755 :
756 :
757 89108 : subroutine set_molar_abundances(n, dbg, ierr)
758 : type (net_info) :: n
759 : type(net_general_info), pointer :: g
760 : logical, intent(in) :: dbg
761 : integer, intent(out) :: ierr
762 :
763 : real(dp) :: sum
764 : integer :: i, ci
765 : character (len=256) :: message
766 : include 'formats'
767 89108 : sum = 0
768 89108 : g => n% g
769 920392 : do i = 1, g% num_isos
770 831284 : sum = sum + n% x(i)
771 831284 : ci = g% chem_id(i)
772 831284 : if (ci <= 0) then
773 0 : write(*,*) 'problem in set_molar_abundances'
774 0 : write(*,*) 'i', i
775 0 : write(*,*) 'g% num_isos', g% num_isos
776 0 : write(*,*) 'g% chem_id(i)', g% chem_id(i)
777 0 : call mesa_error(__FILE__,__LINE__,'set_molar_abundances')
778 : end if
779 920392 : n% y(i) = min(1d0, max(n% x(i), 0d0)) / chem_isos% Z_plus_N(ci)
780 : end do
781 :
782 :
783 89108 : return ! let it go even with bad xsum.
784 :
785 :
786 :
787 : if (abs(sum - 1d0) > 1d-2) then
788 : ierr = -1
789 : if (dbg) then
790 : do i = 1, n% g% num_isos
791 : ci = g% chem_id(i)
792 : write(*,2) chem_isos% name(ci), i, n% x(i)
793 : end do
794 : write(*,1) 'abs(sum - 1d0)', abs(sum - 1d0)
795 : end if
796 : return
797 : end if
798 :
799 : end subroutine set_molar_abundances
800 :
801 10 : subroutine do_clean_up_fractions(nzlo, nzhi, species, nz, xa, max_sum_abs, xsum_tol, ierr)
802 : integer, intent(in) :: nzlo, nzhi, species, nz
803 : real(dp), intent(inout) :: xa(:,:) ! (species, nz)
804 : real(dp), intent(in) :: max_sum_abs, xsum_tol
805 : integer, intent(out) :: ierr
806 : integer :: k, op_err
807 10 : ierr = 0
808 10 : if (nzlo == nzhi) then
809 0 : call do_clean1(species, xa(1: species, nzlo), nzlo, max_sum_abs, xsum_tol, ierr)
810 0 : return
811 : end if
812 : !x$OMP PARALLEL DO PRIVATE(k, op_err)
813 12105 : do k = nzlo, nzhi
814 : op_err = 0
815 12095 : call do_clean1(species, xa(1: species, k), k, max_sum_abs, xsum_tol, op_err)
816 12105 : if (op_err /= 0) ierr = op_err
817 : end do
818 : !x$OMP END PARALLEL DO
819 : end subroutine do_clean_up_fractions
820 :
821 :
822 12095 : subroutine do_clean1(species, xa, k, max_sum_abs, xsum_tol, ierr)
823 : use utils_lib
824 : integer, intent(in) :: species, k
825 : real(dp), intent(inout) :: xa(:) ! (species)
826 : real(dp), intent(in) :: max_sum_abs, xsum_tol
827 : integer, intent(out) :: ierr
828 : integer :: j
829 : real(dp) :: xsum
830 : real(dp), parameter :: tiny_x = 1d-99
831 : character (len=256) :: message
832 12095 : if (max_sum_abs > 1) then ! check for crazy values
833 108855 : xsum = sum(abs(xa(1: species)))
834 12095 : if (is_bad(xsum) .or. xsum > max_sum_abs) then
835 0 : ierr = -1
836 0 : return
837 : end if
838 : end if
839 12095 : ierr = 0
840 108855 : do j = 1, species
841 96760 : if (xa(j) < tiny_x) xa(j) = tiny_x
842 108855 : if (xa(j) > 1) xa(j) = 1
843 : end do
844 108855 : xsum = sum(xa(1: species))
845 12095 : if (abs(xsum-1) > xsum_tol) then
846 0 : ierr = -1
847 0 : return
848 : end if
849 108855 : xa(1: species) = xa(1: species)/xsum
850 : end subroutine do_clean1
851 :
852 : end module net_eval
|
