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_burn
21 : use const_def, only: dp, Qconv
22 : use math_lib
23 : use chem_def
24 : use net_def
25 : use utils_lib, only: is_bad,fill_with_NaNs,fill_with_NaNs_2D
26 : use net_burn_support, only: netint
27 : use net_approx21, only : num_reactions_func => num_reactions
28 :
29 : implicit none
30 :
31 : !logical, parameter :: use_ludcmp = .true.
32 : logical, parameter :: use_ludcmp = .false.
33 :
34 : !logical, parameter :: show_mesa_rates = .true.
35 : logical, parameter :: show_mesa_rates = .false.
36 :
37 : !logical, parameter :: report_ierr = .true.
38 : logical, parameter :: report_ierr = .false.
39 :
40 : contains
41 :
42 2 : subroutine get_pointers( &
43 : g, species, num_reactions, &
44 : dratdumdy1, dratdumdy2, dens_dfdy, dmat, i, ierr)
45 : use net_approx21, only: approx21_nrat
46 : type (Net_General_Info), pointer :: g
47 : integer, intent(in) :: species, num_reactions
48 : real(dp), allocatable, dimension(:) :: dratdumdy1, dratdumdy2
49 : real(dp), allocatable, dimension(:,:) :: dens_dfdy, dmat
50 : integer, intent(inout) :: i
51 : integer, intent(out) :: ierr
52 : integer :: sz
53 : include 'formats'
54 2 : ierr = 0
55 :
56 2 : sz = num_reactions
57 2 : if (g% doing_approx21) sz = num_reactions_func(g% add_co56_to_approx21)
58 :
59 2 : allocate(dratdumdy1(1:sz))
60 2 : allocate(dratdumdy2(1:sz))
61 :
62 4 : allocate(dens_dfdy(1:species,1:species))
63 2 : allocate(dmat(1:species,1:species))
64 :
65 2 : if (g% fill_arrays_with_nans) then
66 0 : call fill_with_NaNs(dratdumdy1)
67 0 : call fill_with_NaNs(dratdumdy2)
68 0 : call fill_with_NaNs_2D(dens_dfdy)
69 0 : call fill_with_NaNs_2D(dmat)
70 : end if
71 :
72 2 : end subroutine get_pointers
73 :
74 :
75 2 : subroutine burn_1_zone( &
76 2 : net_handle, eos_handle, species, num_reactions, t_start, t_end, starting_x, &
77 : ntimes, times, log10Ts_f1, log10Rhos_f1, etas_f1, dxdt_source_term, &
78 : rate_factors, weak_rate_factor, &
79 : reaction_Qs, reaction_neuQs, &
80 : screening_mode, &
81 : stptry_in, max_steps, eps, odescal, &
82 : use_pivoting, trace, burn_dbg, burner_finish_substep, &
83 2 : ending_x, eps_nuc_categories, avg_eps_nuc, eps_neu_total, &
84 : nfcn, njac, nstep, naccpt, nrejct, ierr)
85 2 : use num_def
86 : use num_lib
87 : use mtx_lib
88 : use mtx_def
89 : use rates_def, only: rates_reaction_id_max, reaction_Name, reaction_categories
90 : use rates_lib, only: rates_reaction_id
91 : use net_initialize, only: setup_net_info
92 : use chem_lib, only: basic_composition_info, get_Q
93 : use net_approx21, only: approx21_nrat
94 :
95 : integer, intent(in) :: net_handle, eos_handle
96 : integer, intent(in) :: species
97 : integer, intent(in) :: num_reactions
98 : real(dp), intent(in) :: t_start, t_end, starting_x(:) ! (species)
99 : integer, intent(in) :: ntimes ! ending time is times(num_times); starting time is 0
100 : real(dp), pointer, intent(in) :: times(:) ! (num_times)
101 : real(dp), pointer, intent(in) :: log10Ts_f1(:)
102 : ! =(4,numtimes) interpolant for log10T(time)
103 : real(dp), pointer, intent(in) :: log10Rhos_f1(:)
104 : ! =(4,numtimes) interpolant for log10Rho(time)
105 : real(dp), pointer, intent(in) :: etas_f1(:)
106 : ! =(4,numtimes) interpolant for eta(time)
107 : real(dp), pointer, intent(in) :: dxdt_source_term(:)
108 : ! (species) or null if no source term.
109 : real(dp), intent(in), pointer :: rate_factors(:) ! (num_reactions)
110 : real(dp), intent(in) :: weak_rate_factor
111 : real(dp), pointer, intent(in) :: reaction_Qs(:) ! (rates_reaction_id_max)
112 : real(dp), pointer, intent(in) :: reaction_neuQs(:) ! (rates_reaction_id_max)
113 : integer, intent(in) :: screening_mode ! see screen_def
114 : real(dp), intent(in) :: stptry_in
115 : integer, intent(in) :: max_steps ! maximal number of allowed steps.
116 : real(dp), intent(in) :: eps, odescal ! tolerances. e.g., set both to 1d-6
117 : logical, intent(in) :: use_pivoting
118 : logical, intent(in) :: trace, burn_dbg
119 : interface
120 : include 'burner_finish_substep.inc'
121 : end interface
122 : real(dp), intent(inout) :: ending_x(:) ! (species)
123 : real(dp), intent(inout) :: eps_nuc_categories(:) ! (num_categories)
124 : real(dp), intent(out) :: avg_eps_nuc, eps_neu_total
125 : integer, intent(out) :: nfcn ! number of function evaluations
126 : integer, intent(out) :: njac ! number of jacobian evaluations
127 : integer, intent(out) :: nstep ! number of computed steps
128 : integer, intent(out) :: naccpt ! number of accepted steps
129 : integer, intent(out) :: nrejct ! number of rejected steps
130 : integer, intent(out) :: ierr
131 :
132 : type (Net_General_Info), pointer :: g
133 : integer :: ijac, lrd, lid, lout, i, j, ir, idid, sz
134 : logical :: okay, have_set_rate_screened
135 4 : real(dp) :: temp, rho, eta, lgT, lgRho, r, prev_lgRho, prev_lgT
136 :
137 : integer :: stpmax, imax_dydx, nstp
138 : real(dp) :: &
139 2 : h, start, stptry, stpmin, stopp, max_dydx, abs_max_dydx, &
140 44 : burn_ergs, dx
141 :
142 132 : real(dp), dimension(species), target :: starting_y_a, ending_y_a, save_x_a
143 2 : real(dp), dimension(:), pointer :: starting_y, ending_y, save_x
144 2 : real(dp), dimension(:), allocatable :: dratdumdy1, dratdumdy2
145 :
146 2 : real(dp), dimension(:,:), allocatable :: dens_dfdy, dmat
147 :
148 2 : real(dp) :: xh, xhe, z, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
149 44 : real(dp) :: aion(species)
150 :
151 : logical :: dbg
152 :
153 2 : type (Net_Info) :: n
154 :
155 : integer :: iwork, cid
156 :
157 : include 'formats'
158 :
159 : !dbg = .true.
160 2 : dbg = burn_dbg
161 :
162 2 : if (dbg) then
163 0 : do i=1,species
164 0 : write(*,2) 'starting_x', i, starting_x(i)
165 : end do
166 : end if
167 :
168 2 : starting_y => starting_y_a
169 2 : ending_y => ending_y_a
170 2 : save_x => save_x_a
171 :
172 2 : have_set_rate_screened = .false.
173 :
174 2 : lgT = log10Ts_f1(1)
175 2 : temp = exp10(lgT)
176 2 : lgRho = log10Rhos_f1(1)
177 2 : rho = exp10(lgRho)
178 2 : eta = etas_f1(1)
179 2 : prev_lgT = -1d99
180 2 : prev_lgRho = -1d99
181 :
182 : ierr = 0
183 2 : call get_net_ptr(net_handle, g, ierr)
184 2 : if (ierr /= 0) then
185 : if (report_ierr) write(*,*) 'invalid handle for burn_1_zone'
186 : return
187 : end if
188 :
189 2 : if (g% num_isos /= species) then
190 0 : write(*,*) 'invalid species', species
191 0 : return
192 : end if
193 :
194 2 : if (g% num_reactions /= num_reactions) then
195 0 : write(*,*) 'invalid num_reactions', num_reactions
196 0 : return
197 : end if
198 :
199 2 : nfcn = 0
200 2 : njac = 0
201 2 : nstep = 0
202 2 : naccpt = 0
203 2 : nrejct = 0
204 :
205 44 : do i=1,species
206 42 : cid = g% chem_id(i)
207 42 : if (cid < 0) cid = g% approx21_ye_iso
208 42 : aion(i) = chem_isos% Z_plus_N(cid)
209 42 : save_x(i) = starting_x(i)
210 42 : ending_x(i) = starting_x(i)
211 42 : starting_y(i) = starting_x(i)/aion(i)
212 44 : ending_y(i) = starting_y(i)
213 : end do
214 :
215 2 : start = t_start
216 2 : stptry = stptry_in
217 2 : if (stptry == 0d0) stptry = t_end
218 :
219 : !write(*,1) 'stptry', stptry
220 :
221 2 : stpmin = min(t_end*1d-20,stptry*1d-6)
222 2 : stopp = t_end
223 2 : stpmax = max_steps
224 :
225 2 : n% screening_mode = screening_mode
226 2 : n% g => g
227 :
228 2 : if (dbg) write(*,*) 'call setup_net_info'
229 2 : call setup_net_info(n)
230 2 : if (dbg) write(*,*) 'done setup_net_info'
231 :
232 0 : if (dbg) write(*,*) 'call get_pointers'
233 : call get_pointers( &
234 : g, species, num_reactions, &
235 2 : dratdumdy1, dratdumdy2, dens_dfdy, dmat, iwork, ierr)
236 2 : if (ierr /= 0) return
237 :
238 : call basic_composition_info( &
239 : species, g% chem_id, starting_x, xh, xhe, z, &
240 2 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
241 :
242 2 : stptry = max(start * 1.0d-10,1.0d-16)
243 2 : stpmin = stptry * 1.0d-12
244 :
245 2 : if (dbg) write(*,*) 'call netint'
246 : call netint( &
247 : start,stptry,stpmin,max_steps,stopp,ending_y, &
248 : eps,species,species,naccpt,nrejct,nstep,odescal,dens_dfdy,dmat, &
249 2 : burner_derivs,burner_jakob,burner_finish_substep,ierr)
250 2 : if (dbg) write(*,*) 'done netint'
251 2 : if (ierr /= 0) then
252 : return
253 : !write(*,*) 'netint ierr'
254 : !stop
255 : end if
256 :
257 44 : do i=1,species
258 44 : ending_x(i) = ending_y(i)*aion(i)
259 : end do
260 :
261 : ! set burn_ergs according to change in abundances
262 2 : burn_ergs = 0
263 44 : do i=1,species
264 42 : ending_x(i) = ending_y(i)*aion(i)
265 42 : dx = ending_x(i) - save_x(i)
266 : !write(*,2) 'dx aion end_x', i, dx, aion(i), ending_x(i)
267 42 : cid = g% chem_id(i)
268 : burn_ergs = burn_ergs + &
269 44 : (get_Q(chem_isos,cid))*dx/chem_isos% Z_plus_N(cid)
270 : end do
271 2 : burn_ergs = burn_ergs*Qconv
272 : !write(*,1) 'burn_ergs', burn_ergs
273 4 : avg_eps_nuc = burn_ergs/(t_end - t_start) - eps_neu_total
274 :
275 : contains
276 :
277 3062 : subroutine burner_derivs(x,y,f,species,ierr)
278 : integer, intent(in) :: species
279 : real(dp) :: x, y(:), f(:)
280 : integer, intent(out) :: ierr
281 : integer, parameter :: ld_dfdx = 0
282 3062 : real(dp), target :: dfdx_arry(ld_dfdx,species)
283 : real(dp), pointer :: dfdx(:,:)
284 : real(dp) :: dxdt_sum, dxdt_sum_approx21, &
285 : Z_plus_N, xsum, r, r1, r2
286 : integer :: i, ir, ci, j, k, ibad
287 : logical :: okay
288 :
289 : real(dp), target :: f21_a(species)
290 : real(dp), pointer :: f21(:)
291 :
292 : include 'formats'
293 :
294 : ierr = 0
295 3062 : nfcn = nfcn + 1
296 3062 : dfdx => dfdx_arry
297 3062 : call jakob_or_derivs(x,y,f,dfdx,ierr)
298 3062 : if (ierr /= 0) return
299 :
300 2 : end subroutine burner_derivs
301 :
302 54 : subroutine burner_jakob(x,y,dfdy,species,ierr)
303 : integer, intent(in) :: species
304 : real(dp) :: x, y(:)
305 : real(dp) :: dfdy(:,:)
306 : integer, intent(out) :: ierr
307 : real(dp), target :: f_arry(0)
308 : real(dp), pointer :: f(:)
309 :
310 : real(dp) :: Z_plus_N, df_t, df_m
311 : integer :: i, ci, j, cj
312 : logical :: okay
313 : include 'formats'
314 :
315 : ierr = 0
316 54 : njac = njac + 1
317 54 : f => f_arry
318 :
319 54 : call jakob_or_derivs(x,y,f,dfdy,ierr)
320 54 : if (ierr /= 0) return
321 :
322 : end subroutine burner_jakob
323 :
324 3116 : subroutine jakob_or_derivs(time,y,f,dfdy,ierr)
325 : use chem_lib, only: basic_composition_info
326 : use chem_def, only: chem_isos, num_categories, category_name, ih1
327 : use net_eval, only: eval_net
328 : use rates_def, only: rates_reaction_id_max, i_rate, i_rate_dT, i_rate_dRho
329 : use interp_1d_lib, only: interp_value
330 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results, i_eta
331 : use eos_lib, only: eosDT_get
332 :
333 : real(dp) :: time, y(:), f(:)
334 : real(dp) :: dfdy(:,:)
335 : integer, intent(out) :: ierr
336 :
337 3116 : real(dp) :: rho, lgRho, T, lgT, rate_limit, rat, dratdt, dratdd
338 3116 : real(dp) :: eta, d_eta_dlnT, d_eta_dlnRho
339 3116 : real(dp) :: eps_nuc
340 3116 : real(dp) :: d_eps_nuc_dT
341 3116 : real(dp) :: d_eps_nuc_dRho
342 68552 : real(dp) :: d_eps_nuc_dx(species)
343 68552 : real(dp) :: dxdt(species)
344 68552 : real(dp) :: d_dxdt_dRho(species)
345 68552 : real(dp) :: d_dxdt_dT(species)
346 1442708 : real(dp) :: d_dxdt_dx(species, species)
347 :
348 : logical :: rates_only, dxdt_only, okay
349 : integer :: i, j, k, ir
350 :
351 3116 : real(dp), pointer, dimension(:) :: actual_Qs, actual_neuQs
352 3116 : real(dp) :: xsum
353 3116 : logical, pointer :: from_weaklib(:)
354 :
355 1511260 : real(dp), target :: x_a(species), dfdx_a(species,species)
356 3116 : real(dp), pointer :: x(:), dfdx(:,:)
357 :
358 249280 : real(dp), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
359 199424 : real(dp) :: d_dxa(num_eos_d_dxa_results,species)
360 :
361 : include 'formats'
362 :
363 3116 : ierr = 0
364 :
365 3116 : x => x_a
366 3116 : dfdx => dfdx_a
367 :
368 3116 : actual_Qs => null()
369 3116 : actual_neuQs => null()
370 3116 : from_weaklib => null()
371 :
372 3116 : if (ntimes == 1) then
373 :
374 0 : lgT = log10Ts_f1(1)
375 0 : lgRho = log10Rhos_f1(1)
376 :
377 : else
378 :
379 3116 : call interp_value(times, ntimes, log10Ts_f1, time, lgT, ierr)
380 3116 : if (ierr /= 0) then
381 : if (report_ierr) &
382 : write(*,1) 'interp_value for lgT failed in jakob for 1 zone burn', time
383 : return
384 : end if
385 :
386 3116 : call interp_value(times, ntimes, log10Rhos_f1, time, lgRho, ierr)
387 3116 : if (ierr /= 0) then
388 : if (report_ierr) &
389 : write(*,1) 'interp_value for lgRho failed in jakob for 1 zone burn', time
390 : return
391 : end if
392 :
393 : end if
394 :
395 3116 : xsum = 0
396 68552 : do i=1,species
397 65436 : if (is_bad(y(i))) then
398 0 : ierr = -1
399 : if (report_ierr) &
400 : write(*,2) 'net_burn failed in jakob_or_derivs: bad y(i) lgT lgRho', i, y(i), lgT, lgRho
401 0 : return
402 : stop
403 : end if
404 65436 : y(i) = min(1.0d0, max(y(i),1.0d-30))
405 65436 : x(i) = y(i)*aion(i)
406 68552 : xsum = xsum + x(i)
407 : end do
408 3116 : if (trace .and. xsum > 2) write(*,*) 'sum_x, time', xsum, time
409 :
410 3116 : rho = exp10(lgRho)
411 3116 : T = exp10(lgT)
412 :
413 : call basic_composition_info( &
414 : species, g% chem_id, x, xh, xhe, z, &
415 3116 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
416 :
417 :
418 : call eosDT_get( &
419 : eos_handle, species, g% chem_id, g% net_iso, x, &
420 : Rho, lgRho, T, lgT, &
421 3116 : res, d_dlnd, d_dlnT, d_dxa, ierr)
422 3116 : if (ierr /= 0) then
423 : if (report_ierr) write(*,*) 'failed in eosDT_get'
424 : return
425 : end if
426 3116 : eta = res(i_eta)
427 3116 : d_eta_dlnT = d_dlnT(i_eta)
428 3116 : d_eta_dlnRho = d_dlnd(i_eta)
429 :
430 :
431 3116 : rates_only = .false.
432 3116 : dxdt_only = (size(dfdy,dim=1) == 0)
433 :
434 : call eval_net( &
435 : n, g, rates_only, dxdt_only, &
436 : species, num_reactions, g% num_wk_reactions, &
437 : x, T, lgT, rho, lgRho, &
438 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
439 : rate_factors, weak_rate_factor, &
440 : reaction_Qs, reaction_neuQs, &
441 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
442 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
443 : screening_mode, &
444 : eps_nuc_categories, eps_neu_total, &
445 3116 : actual_Qs, actual_neuQs, from_weaklib, .false., ierr)
446 :
447 3116 : if (size(f,dim=1) > 0) then
448 67364 : do j = 1, species
449 64302 : f(j) = dxdt(j)/aion(j)
450 3062 : if (.false. .and. is_bad(f(j))) then
451 : write(*,1) 'x', x
452 : write(*,2) 'f(j)', j, f(j)
453 : call mesa_error(__FILE__,__LINE__,'jakob_or_derivs')
454 : end if
455 : end do
456 : end if
457 :
458 3116 : if (.not. dxdt_only) then
459 1188 : do j = 1, species
460 25002 : do i = 1, species
461 23814 : dfdy(i,j) = d_dxdt_dx(i,j)*aion(j)/aion(i)
462 1134 : if (.false. .and. is_bad(dfdy(i,j))) then
463 : write(*,1) 'x', x
464 : write(*,3) 'dfdy(i,j)', i, j, dfdy(i,j)
465 : call mesa_error(__FILE__,__LINE__,'jakob_or_derivs')
466 : end if
467 : end do
468 : end do
469 : end if
470 :
471 6232 : end subroutine jakob_or_derivs
472 :
473 : end subroutine burn_1_zone
474 :
475 : end module net_burn
|
