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