Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012-2019 Bill Paxton & 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 mod_one_zone_support
21 :
22 : use chem_def
23 : use chem_lib
24 : use math_lib
25 : use net_def
26 : use net_lib
27 : use const_def, only: dp, i8, Qconv, secyer, kerg, avo, ln10
28 : use rates_def
29 : use utils_lib, only: mesa_error
30 :
31 : implicit none
32 :
33 : character(len=256) :: net_name
34 :
35 : character(len=256):: final_abundances_filename
36 : logical :: save_final_abundances
37 :
38 : character(len=256):: initial_abundances_filename
39 : logical :: read_initial_abundances
40 :
41 : character(len=256):: T_Rho_history_filename
42 : logical :: read_T_Rho_history
43 :
44 : character(len=256):: burn_filename
45 : real(dp) :: burn_tend, burn_rho, burn_temp, &
46 : burn_rtol, burn_atol, burn_P, burn_xmin, burn_xmax, &
47 : eps, odescal, stptry
48 : logical :: trace, burn_dbg, use_pivoting
49 : real(dp) :: min_for_show_peak_abundances
50 : integer :: max_num_for_show_peak_abundances
51 :
52 : integer, parameter :: max_num_burn_isos_to_show = 1000
53 : character(len=iso_name_length) :: names_of_isos_to_show(max_num_burn_isos_to_show)
54 : integer :: num_names_of_isos_to_show
55 :
56 : integer, parameter :: max_num_isos_for_Xinit = 1000
57 : character(len=iso_name_length) :: names_of_isos_for_Xinit(max_num_isos_for_Xinit)
58 : real(dp) :: values_for_Xinit(max_num_isos_for_Xinit)
59 : integer :: num_isos_for_Xinit
60 : logical :: uniform_Xinit
61 :
62 : integer :: screening_mode
63 :
64 : integer, parameter :: io_out = 35
65 : real(dp) :: data_output_min_t
66 :
67 : character(len=32) :: which_solver
68 : integer :: decsol_switch
69 : character(len=32) :: small_mtx_decsol, large_mtx_decsol
70 :
71 : logical :: show_net_reactions_info
72 :
73 : real(dp) :: rattab_logT_lower_bound, rattab_logT_upper_bound
74 :
75 : character(len=256):: data_filename, data_heading_line
76 : character(len=64) :: net_file, cache_suffix
77 :
78 : integer :: handle, eos_handle, net_handle
79 : type(Net_General_Info), pointer :: g
80 : integer :: species, num_reactions
81 : integer, pointer :: reaction_id(:)
82 : integer, dimension(:), pointer :: net_iso, chem_id
83 :
84 : real(dp) :: z, abar, zbar, z2bar, z53bar, ye, eps_neu_total
85 : real(dp) :: eta, d_eta_dlnT, d_eta_dlnRho
86 : real(dp), dimension(:), pointer :: &
87 : xin, xin_copy, d_eps_nuc_dx, dxdt, d_dxdt_dRho, d_dxdt_dT
88 : real(dp), pointer :: d_dxdt_dx(:, :)
89 :
90 : real(dp) :: weak_rate_factor
91 :
92 : integer :: max_steps ! maximal number of allowed steps.
93 :
94 : real(dp), pointer :: x_initial(:) ! (species)
95 :
96 : real(dp) :: burn_lnE, burn_lnS
97 : real(dp) :: burn_logT, burn_logRho, &
98 : burn_eta, burn_deta_dlnT, burn_Cv, burn_d_Cv_dlnT
99 :
100 : real(dp) :: T_prev, time_prev, eps_nuc_prev, eps_neu_prev, cp_prev
101 : real(dp), pointer :: x_previous(:) ! (species)
102 :
103 : real(dp), dimension(:), pointer :: peak_abundance, peak_time
104 :
105 : integer :: num_times_for_burn
106 : integer, parameter :: max_num_times = 10000
107 : real(dp), dimension(max_num_times) :: &
108 : times_for_burn, log10Ts_for_burn, log10Rhos_for_burn, &
109 : etas_for_burn, log10Ps_for_burn
110 :
111 : logical :: burn_at_constant_P, clip, show_peak_x_and_time, &
112 : burn_at_constant_density
113 : real(dp) :: starting_temp, pressure
114 :
115 : real(dp) :: max_step_size ! maximal step size.
116 :
117 : real(dp), pointer :: rate_factors(:) ! (num_reactions)
118 : integer, pointer :: net_reaction_ptr(:)
119 :
120 : integer, parameter :: max_num_reactions_to_track = 100
121 : integer :: num_reactions_to_track
122 : character(len=maxlen_reaction_Name) :: &
123 : reaction_to_track(max_num_reactions_to_track)
124 : integer :: index_for_reaction_to_track(max_num_reactions_to_track)
125 :
126 : integer, parameter :: max_num_special_rate_factors = 100
127 : integer :: num_special_rate_factors
128 : real(dp) :: special_rate_factor(max_num_special_rate_factors)
129 : character(len=maxlen_reaction_Name) :: &
130 : reaction_for_special_factor(max_num_special_rate_factors)
131 :
132 : character(len=16) :: set_rate_c12ag, set_rate_n14pg, set_rate_3a, &
133 : set_rate_1212
134 :
135 : logical :: show_Qs, quiet, complete_silence_please, &
136 : show_ye_stuff
137 :
138 : real(dp) :: starting_logT
139 :
140 : logical, parameter :: dbg = .false.
141 :
142 : contains
143 :
144 6 : integer function burn_isos_for_Xinit(i)
145 : integer, intent(in) :: i
146 6 : burn_isos_for_Xinit = chem_get_iso_id(names_of_isos_for_Xinit(i))
147 6 : end function burn_isos_for_Xinit
148 :
149 27000 : integer function burn_isos_to_show(i)
150 : integer, intent(in) :: i
151 27000 : burn_isos_to_show = chem_get_iso_id(names_of_isos_to_show(i))
152 27000 : end function burn_isos_to_show
153 :
154 4 : subroutine Do_One_Zone_Burn(net_file_in)
155 : use num_lib, only: solver_option
156 : use mtx_lib, only: decsol_option
157 : use eos_lib
158 : use mtx_def
159 : use interp_1d_lib, only: interp_pm
160 : use interp_1d_def, only: pm_work_size
161 : use test_net_support, only: Setup_eos
162 : use net_lib, only: get_net_reaction_table_ptr
163 : use rates_lib, only: rates_reaction_id
164 : use utils_lib, only: set_nan
165 :
166 : character(len=*), intent(in) :: net_file_in
167 :
168 : character(len=256) :: net_file
169 4 : real(dp) :: logRho, logT, Rho, T, xsum, &
170 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT
171 100 : real(dp), target :: eps_nuc_categories(num_categories)
172 :
173 : integer :: i, j, ierr, iounit, decsol_choice, solver_choice, num_times
174 :
175 4 : real(dp), dimension(:), pointer :: times, dxdt_source_term
176 : real(dp), dimension(:), pointer :: &
177 4 : log10Ts_f1, log10Rhos_f1, etas_f1, log10Ps_f1
178 : real(dp), dimension(:, :), pointer :: &
179 4 : log10Ts_f, log10Rhos_f, etas_f, log10Ps_f
180 :
181 : ! args to control the solver -- see num/public/num_isolve.dek
182 4 : real(dp) :: h
183 : ! absolute and relative error tolerances
184 4 : real(dp), pointer :: rtol(:) ! relative error tolerance (species)
185 4 : real(dp), pointer :: atol(:) ! absolute error tolerance (species)
186 : integer :: itol ! switch for rtol and atol
187 :
188 4 : real(dp), pointer :: ending_x(:) ! (species)
189 : integer :: nfcn ! number of function evaluations
190 : integer :: njac ! number of jacobian evaluations
191 : integer :: nstep ! number of computed steps
192 : integer :: naccpt ! number of accepted steps
193 : integer :: nrejct ! number of rejected steps
194 : integer :: max_order_used
195 :
196 : integer :: iout, caller_id, cid, ir
197 : integer(i8) :: time0, time1, clock_rate
198 :
199 8 : real(dp) :: ending_temp, ending_rho, ending_lnS, initial_rho, initial_lnS, dt
200 8 : real(dp) :: ending_log10T, starting_log10T, avg_eps_nuc, ending_eps_neu_total
201 4 : real(dp) :: time_doing_net, time_doing_eos, told
202 8 : real(dp) :: xh, xhe, z, abar, zbar, z2bar, z53bar, ye, mass_correction
203 :
204 : integer, parameter :: nwork = pm_work_size
205 4 : real(dp), pointer :: pm_work(:)
206 :
207 : type(Net_Info) :: n
208 : type(Net_General_Info), pointer :: g
209 :
210 : include 'formats'
211 :
212 4 : ierr = 0
213 4 : told = 0
214 :
215 4 : net_file = net_file_in
216 :
217 4 : call test_net_setup(net_file)
218 :
219 4 : net_handle = handle
220 4 : call get_net_ptr(net_handle, g, ierr)
221 4 : if (ierr /= 0) return
222 :
223 4 : call Setup_eos(eos_handle)
224 : ! g% max_rate_times_dt = max_rate_times_dt
225 :
226 4 : logT = burn_logT
227 4 : T = burn_temp
228 4 : logRho = burn_logRho
229 4 : Rho = burn_rho
230 :
231 4 : if (read_T_Rho_history) then
232 0 : num_times = max_num_times
233 4 : else if (num_times_for_burn <= 0) then
234 2 : num_times = 1
235 : else
236 2 : num_times = num_times_for_burn
237 : end if
238 :
239 4 : if (num_names_of_isos_to_show < 0) num_names_of_isos_to_show = species
240 :
241 : allocate ( &
242 : rate_factors(num_reactions), rtol(species), atol(species), &
243 : x_initial(species), x_previous(species), ending_x(species), times(num_times), &
244 : log10Ts_f1(4*num_times), log10Rhos_f1(4*num_times), &
245 : etas_f1(4*num_times), log10Ps_f1(4*num_times), &
246 : peak_abundance(species), peak_time(species), &
247 4 : stat=ierr)
248 4 : if (ierr /= 0) then
249 0 : write (*, *) 'allocate failed for Do_One_Zone_Burn'
250 0 : call mesa_error(__FILE__, __LINE__)
251 : end if
252 :
253 4 : call get_net_reaction_table_ptr(net_handle, net_reaction_ptr, ierr)
254 4 : if (ierr /= 0) then
255 0 : write (*, *) 'bad net? get_net_reaction_table_ptr failed'
256 0 : return
257 : end if
258 :
259 376 : rate_factors(:) = 1
260 :
261 4 : if (num_special_rate_factors > 0) then
262 0 : do i = 1, num_special_rate_factors
263 0 : if (len_trim(reaction_for_special_factor(i)) == 0) cycle
264 0 : ir = rates_reaction_id(reaction_for_special_factor(i))
265 0 : j = 0
266 0 : if (ir > 0) j = net_reaction_ptr(ir)
267 0 : if (j <= 0) cycle
268 0 : rate_factors(j) = special_rate_factor(i)
269 : write (*, 1) 'set special rate factor for '// &
270 0 : trim(reaction_for_special_factor(i)), special_rate_factor(i)
271 : end do
272 : end if
273 :
274 4 : if (num_reactions_to_track > 0) then
275 0 : do i = 1, num_reactions_to_track
276 0 : index_for_reaction_to_track(i) = 0
277 0 : if (len_trim(reaction_to_track(i)) == 0) cycle
278 0 : ir = rates_reaction_id(reaction_to_track(i))
279 0 : j = 0
280 0 : if (ir > 0) j = net_reaction_ptr(ir)
281 0 : if (j <= 0) cycle
282 0 : index_for_reaction_to_track(i) = j
283 0 : write (*, 1) 'track rate '//trim(reaction_to_track(i))
284 : end do
285 : end if
286 :
287 4 : log10Ts_f(1:4, 1:num_times) => log10Ts_f1(1:4*num_times)
288 4 : log10Rhos_f(1:4, 1:num_times) => log10Rhos_f1(1:4*num_times)
289 4 : etas_f(1:4, 1:num_times) => etas_f1(1:4*num_times)
290 4 : log10Ps_f(1:4, 1:num_times) => log10Ps_f1(1:4*num_times)
291 :
292 88 : peak_abundance(:) = 0
293 :
294 88 : xin = 0
295 4 : eta = 0
296 :
297 4 : iout = 1
298 4 : itol = 0
299 :
300 88 : rtol(:) = burn_rtol
301 88 : atol(:) = burn_atol
302 :
303 88 : xin = 0
304 4 : if (read_initial_abundances) then
305 0 : call read_X(ierr)
306 0 : if (ierr /= 0) return
307 4 : else if (uniform_Xinit) then
308 0 : xin(:) = 0.5d0/(species - 1)
309 0 : j = net_iso(ih1)
310 0 : if (j <= 0) call mesa_error(__FILE__, __LINE__, 'where is the h?')
311 0 : xin(j) = 0.5d0
312 : else
313 10 : do i = 1, num_isos_for_Xinit
314 6 : cid = burn_isos_for_Xinit(i)
315 6 : j = net_iso(cid)
316 6 : if (j <= 0 .or. j > species) then
317 : write (*, *) 'bad names_of_isos_for_Xinit '// &
318 0 : trim(names_of_isos_for_Xinit(i))
319 0 : call mesa_error(__FILE__, __LINE__)
320 : end if
321 10 : xin(j) = values_for_Xinit(i)
322 : end do
323 : end if
324 :
325 : !xin(:) = xin(:)/sum(xin(:))
326 :
327 4 : if (read_T_Rho_history) then
328 0 : call do_read_T_Rho_history(ierr)
329 0 : if (ierr /= 0) return
330 : end if
331 :
332 4 : if (num_times_for_burn <= 0) then
333 2 : times(1) = burn_tend
334 2 : log10Ts_f(1, 1) = logT
335 8 : log10Ts_f(2:4, 1) = 0
336 2 : log10Rhos_f(1, 1) = logRho
337 8 : log10Rhos_f(2:4, 1) = 0
338 4 : etas_f(1:1, 1) = burn_eta
339 8 : etas_f(2:4, 1) = 0
340 : else
341 8 : times(1:num_times) = times_for_burn(1:num_times)
342 8 : log10Ts_f(1, 1:num_times) = log10Ts_for_burn(1:num_times)
343 26 : log10Ts_f(2:4, 1:num_times) = 0
344 8 : log10Rhos_f(1, 1:num_times) = log10Rhos_for_burn(1:num_times)
345 26 : log10Rhos_f(2:4, 1:num_times) = 0
346 8 : etas_f(1, 1:num_times) = etas_for_burn(1:num_times)
347 26 : etas_f(2:4, 1:num_times) = 0
348 2 : logRho = log10Rhos_for_burn(1)
349 2 : burn_rho = exp10(logRho)
350 2 : logT = log10Ts_for_burn(1)
351 2 : burn_temp = exp10(logT)
352 2 : burn_tend = times_for_burn(num_times)
353 : if (.false.) then
354 : write (*, '(A)')
355 : do i = 1, num_times
356 : write (*, 2) 'history', i, times_for_burn(i), log10Ts_for_burn(i), &
357 : log10Rhos_for_burn(i), etas_for_burn(i)
358 : end do
359 : write (*, '(A)')
360 : end if
361 : end if
362 4 : starting_logT = logT
363 :
364 4 : h = 1d-2*burn_tend ! 1d-14
365 : !write(*,1) 'h', h
366 : !stop
367 :
368 172 : x_initial(1:species) = xin(1:species)
369 172 : x_previous(1:species) = xin(1:species)
370 4 : caller_id = 0
371 4 : dxdt_source_term => null()
372 :
373 4 : if (.not. quiet) then
374 0 : write (*, '(A)')
375 0 : write (*, '(A)')
376 0 : write (*, 1) 'one zone burn '//trim(net_file)
377 0 : write (*, '(A)')
378 0 : write (*, 2) 'number of species', species
379 0 : write (*, '(A)')
380 0 : write (*, 1) 'temp', burn_temp
381 0 : write (*, 1) 'logT', burn_logT
382 0 : write (*, 1) 'rho', burn_rho
383 0 : write (*, 1) 'logRho', burn_logRho
384 0 : write (*, 1) 'eta', burn_eta
385 0 : write (*, '(A)')
386 0 : write (*, 1) 'tend', burn_tend
387 0 : write (*, 1) 'tend/secyer', burn_tend/secyer
388 0 : write (*, '(A)')
389 0 : write (*, 1) 'initial abundances'
390 0 : call show_X(xin, .false., .false.)
391 : end if
392 :
393 : ! data_heading_line was not set and writing out nulls. change it - fxt
394 : ! write(io_out,'(a)') trim(data_heading_line)
395 4 : write (data_heading_line, '(99(a,1pe14.6))') 'temp =', burn_temp, ' rho =', burn_rho
396 4 : write (io_out, '(a)') trim(data_heading_line)
397 :
398 : write (io_out, '(a7,99(a26,1x))', advance='no') &
399 4 : 'i', &
400 4 : 'signed_log_avg_eps_nuc', &
401 4 : 'avg_eps_nuc', &
402 4 : 'signed_log_eps_nuc', &
403 4 : 'eps_nuc', &
404 4 : 'eps_neu', &
405 4 : 'delta_logT', &
406 4 : 'logT', &
407 4 : 'logRho', &
408 4 : 'logPgas', &
409 4 : 'logP', &
410 4 : 'abar', &
411 4 : 'zbar', &
412 4 : 'entropy', &
413 4 : 'logE', &
414 4 : 'time', &
415 4 : 'lg_time', &
416 4 : 'lg_yrs', &
417 4 : 'dt', &
418 4 : 'lg_dt', &
419 4 : 'ye', &
420 8 : 'xsum_sub_1'
421 :
422 64 : do i = 1, num_names_of_isos_to_show
423 60 : if (num_names_of_isos_to_show < species) then
424 60 : cid = burn_isos_to_show(i)
425 : else
426 0 : if (i > species) exit
427 0 : cid = chem_id(i)
428 : end if
429 60 : j = net_iso(cid)
430 60 : if (j == 0) cycle
431 64 : write (io_out, '(a26,1x)', advance='no') 'lg_'//trim(chem_isos%name(cid))
432 : end do
433 64 : do i = 1, num_names_of_isos_to_show
434 60 : if (num_names_of_isos_to_show < species) then
435 60 : cid = burn_isos_to_show(i)
436 : else
437 0 : if (i > species) exit
438 0 : cid = chem_id(i)
439 : end if
440 60 : j = net_iso(cid)
441 60 : if (j == 0) cycle
442 64 : write (io_out, '(a26,1x)', advance='no') trim(chem_isos%name(cid))
443 : end do
444 4 : do i = 1, num_reactions_to_track
445 0 : j = index_for_reaction_to_track(i)
446 0 : if (j == 0) cycle
447 0 : write (io_out, '(a26,1x)', advance='no') 'eps_'//trim(reaction_to_track(i))
448 0 : write (io_out, '(a26,1x)', advance='no') 'raw_'//trim(reaction_to_track(i))
449 4 : write (io_out, '(a26,1x)', advance='no') 'scrn_'//trim(reaction_to_track(i))
450 : end do
451 4 : write (io_out, *)
452 :
453 4 : if (show_net_reactions_info) then
454 0 : write (*, '(a)') ' species'
455 0 : do j = 1, species
456 0 : cid = chem_id(j)
457 0 : write (*, '(a)') chem_isos%name(cid)
458 : end do
459 0 : do j = 1, species
460 0 : cid = chem_id(j)
461 0 : write (*, '(3i5,3x,a)') j, &
462 0 : chem_isos%Z(cid), &
463 0 : chem_isos%N(cid), &
464 0 : chem_isos%name(cid)
465 : end do
466 0 : write (*, '(A)')
467 :
468 0 : call show_net_reactions_and_info(handle, 6, ierr)
469 0 : if (ierr /= 0) then
470 0 : write (*, *) 'failed in show_net_reactions'
471 0 : call mesa_error(__FILE__, __LINE__)
472 : end if
473 0 : write (*, '(A)')
474 : end if
475 :
476 4 : if (.not. quiet) then
477 0 : write (*, 1) 'h', h
478 0 : write (*, 1) 'max_step_size', max_step_size
479 0 : write (*, 2) 'max_steps', max_steps
480 0 : write (*, 2) 'screening_mode', screening_mode
481 0 : write (*, '(A)')
482 : end if
483 :
484 4 : if (species >= decsol_switch) then
485 0 : decsol_choice = decsol_option(large_mtx_decsol, ierr)
486 0 : if (ierr /= 0) then
487 0 : write (*, *) 'ERROR: unknown large_mtx_decsol '//trim(large_mtx_decsol)
488 0 : return
489 : end if
490 : else
491 4 : decsol_choice = decsol_option(small_mtx_decsol, ierr)
492 4 : if (ierr /= 0) then
493 0 : write (*, *) 'ERROR: unknown small_mtx_decsol '//trim(small_mtx_decsol)
494 0 : return
495 : end if
496 : end if
497 :
498 4 : solver_choice = solver_option(which_solver, ierr)
499 4 : if (ierr /= 0) then
500 0 : write (*, *) 'ERROR: unknown value for which_solver '//trim(which_solver)
501 0 : return
502 : end if
503 :
504 4 : nullify (pm_work)
505 :
506 4 : call system_clock(time0, clock_rate)
507 4 : time_doing_net = -1
508 4 : time_doing_eos = -1
509 :
510 4 : if (burn_at_constant_density) then
511 :
512 0 : starting_log10T = burn_logT
513 0 : logT = burn_logT
514 0 : logRho = burn_logRho
515 :
516 : call net_1_zone_burn_const_density( &
517 : net_handle, eos_handle, species, num_reactions, &
518 : 0d0, burn_tend, xin, logT, logRho, &
519 : get_eos_info_for_burn_at_const_density, &
520 : rate_factors, weak_rate_factor, &
521 : std_reaction_Qs, std_reaction_neuQs, &
522 : screening_mode, &
523 : stptry, max_steps, eps, odescal, &
524 : use_pivoting, trace, burn_dbg, burn_finish_substep, &
525 : ending_x, eps_nuc_categories, ending_log10T, &
526 : avg_eps_nuc, ending_eps_neu_total, &
527 0 : nfcn, njac, nstep, naccpt, nrejct, ierr)
528 0 : if (ierr == 0 .and. .not. quiet) then
529 0 : write (*, '(A)')
530 0 : write (*, 1) 'constant log10Rho', burn_logRho
531 0 : write (*, 1) 'starting log10T', starting_log10T
532 0 : write (*, 1) 'ending log10T', ending_log10T
533 0 : write (*, 1) 'change in log10T', ending_log10T - starting_log10T
534 0 : write (*, 1) 'change in T', exp10(ending_log10T/starting_log10T)
535 0 : write (*, 1) 'avg_eps_nuc', avg_eps_nuc
536 0 : write (*, 1) 'ending_eps_neu_total', ending_eps_neu_total
537 : end if
538 :
539 4 : else if (burn_at_constant_P) then
540 2 : if (num_times_for_burn <= 0) then
541 2 : log10Ps_f(1, 1) = log10(pressure)
542 10 : log10Ps_f(2:4, 1:num_times) = 0
543 : else
544 0 : log10Ps_f(1, 1:num_times) = log10Ps_for_burn(1:num_times)
545 0 : log10Ps_f(2:4, 1:num_times) = 0
546 0 : pressure = log10Ps_f(1, 1)
547 0 : starting_temp = log10Ts_f(1, 1)
548 : end if
549 2 : if (.not. quiet) then
550 0 : write (*, 1) 'pressure', pressure
551 0 : write (*, 1) 'starting_temp', starting_temp
552 : end if
553 2 : if (num_times > 1) then ! create interpolant
554 0 : allocate (pm_work(num_times*nwork))
555 0 : call interp_pm(times, num_times, log10Ps_f1, nwork, pm_work, 'net_1_zone_burn', ierr)
556 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'failed in interp for logTs')
557 : end if
558 :
559 : call net_1_zone_burn_const_P( &
560 : net_handle, eos_handle, species, num_reactions, &
561 : solver_choice, starting_temp, xin, clip, &
562 : num_times, times, log10Ps_f1, &
563 : rate_factors, weak_rate_factor, &
564 : std_reaction_Qs, std_reaction_neuQs, screening_mode, &
565 : h, max_step_size, max_steps, rtol, atol, itol, burn_xmin, burn_xmax, &
566 : decsol_choice, caller_id, burn_solout, iout, &
567 : ending_x, ending_temp, ending_rho, ending_lnS, initial_rho, initial_lnS, &
568 2 : nfcn, njac, nstep, naccpt, nrejct, time_doing_net, time_doing_eos, ierr)
569 2 : if (ierr == 0 .and. .not. quiet) then
570 0 : write (*, '(A)')
571 0 : write (*, 1) 'pressure', pressure
572 0 : write (*, 1) 'starting_temp', starting_temp
573 0 : write (*, 1) 'ending_temp', ending_temp
574 0 : write (*, 1) 'initial_rho', initial_rho
575 0 : write (*, 1) 'ending_rho', ending_rho
576 0 : write (*, 1) 'initial_entropy', exp(initial_lnS)/(avo*kerg)
577 0 : write (*, 1) 'ending_entropy', exp(ending_lnS)/(avo*kerg)
578 : end if
579 : else
580 2 : if (num_times > 1) then ! create interpolants
581 2 : allocate (pm_work(num_times*nwork))
582 2 : call interp_pm(times, num_times, log10Ts_f1, nwork, pm_work, 'net_1_zone_burn', ierr)
583 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'failed in interp for logTs')
584 2 : call interp_pm(times, num_times, log10Rhos_f1, nwork, pm_work, 'net_1_zone_burn', ierr)
585 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'failed in interp for logRhos')
586 2 : call interp_pm(times, num_times, etas_f1, nwork, pm_work, 'net_1_zone_burn', ierr)
587 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'failed in interp for etas')
588 : end if
589 :
590 : call net_1_zone_burn( &
591 : net_handle, eos_handle, species, num_reactions, 0d0, burn_tend, xin, &
592 : num_times, times, log10Ts_f1, log10Rhos_f1, etas_f1, dxdt_source_term, &
593 : rate_factors, weak_rate_factor, &
594 : std_reaction_Qs, std_reaction_neuQs, &
595 : screening_mode, &
596 : stptry, max_steps, eps, odescal, &
597 : use_pivoting, trace, burn_dbg, burn_finish_substep, &
598 : ending_x, eps_nuc_categories, &
599 : avg_eps_nuc, eps_neu_total, &
600 2 : nfcn, njac, nstep, naccpt, nrejct, ierr)
601 :
602 : end if
603 4 : call system_clock(time1, clock_rate)
604 4 : dt = dble(time1 - time0)/clock_rate
605 :
606 4 : if (ierr /= 0) then
607 0 : write (*, *) 'net_1_zone_burn returned ierr', ierr
608 0 : call mesa_error(__FILE__, __LINE__)
609 : end if
610 :
611 4 : if (.not. quiet) then
612 0 : write (*, '(A)')
613 0 : write (*, '(A)')
614 : end if
615 :
616 4 : if (.not. complete_silence_please) then
617 2 : if (nstep >= max_steps) then
618 0 : write (*, 2) 'hit max number of steps', nstep
619 0 : call mesa_error(__FILE__, __LINE__, 'burn')
620 : end if
621 2 : write (*, 2) 'number of species', species
622 2 : write (*, 1) 'large final abundances', min_for_show_peak_abundances
623 2 : call show_X(ending_x(:), show_peak_x_and_time, .true.)
624 2 : if (show_ye_stuff) then
625 : call basic_composition_info( &
626 : species, chem_id, ending_x, xh, xhe, z, abar, zbar, z2bar, z53bar, ye, &
627 0 : mass_correction, xsum)
628 0 : write (*, 1) 'changes in abundances'
629 0 : do i = 1, species
630 0 : x_previous(i) = ending_x(i) - x_initial(i)
631 : end do
632 0 : min_for_show_peak_abundances = -1d99
633 0 : call show_X(x_previous(:), .false., .true.)
634 0 : write (*, '(A)')
635 0 : write (*, 1) 'ye', ye
636 0 : write (*, '(A)')
637 0 : if (burn_at_constant_density) then
638 0 : write (*, '(A)')
639 0 : write (*, 1) 'constant log10Rho', burn_logRho
640 0 : write (*, 1) 'starting log10T', starting_log10T
641 0 : write (*, 1) 'ending log10T', ending_log10T
642 0 : write (*, 1) 'change in log10T', ending_log10T - starting_log10T
643 0 : write (*, 1) 'change in T', exp10(ending_log10T) - exp10(starting_log10T)
644 0 : write (*, 1) 'ending_eps_neu_total', ending_eps_neu_total
645 0 : write (*, 1) 'avg_eps_nuc', avg_eps_nuc
646 0 : write (*, '(A)')
647 0 : write (*, 1) 'tolerance eps', eps
648 0 : write (*, 1) 'odescal', odescal
649 0 : else if (.not. burn_at_constant_P) then
650 0 : write (*, 1) 'ending eps_neu_total', eps_neu_total
651 0 : write (*, 1) 'avg_eps_nuc', avg_eps_nuc
652 0 : write (*, '(A)')
653 0 : write (*, 1) 'logT', logT
654 0 : write (*, 1) 'logRho', logRho
655 0 : write (*, '(A)')
656 0 : write (*, 1) 'burn_temp', burn_temp
657 0 : write (*, 1) 'burn_rho', burn_rho
658 0 : write (*, '(A)')
659 0 : write (*, 1) 'burn_rtol', burn_rtol
660 0 : write (*, 1) 'burn_atol', burn_atol
661 : else
662 0 : write (*, 1) 'burn_temp', burn_temp
663 0 : write (*, 1) 'burn_rho', burn_rho
664 0 : write (*, '(A)')
665 0 : write (*, 1) 'burn_rtol', burn_rtol
666 0 : write (*, 1) 'burn_atol', burn_atol
667 : end if
668 0 : write (*, '(A)')
669 0 : write (*, 1) 'burn_tend (seconds)', burn_tend
670 0 : write (*, '(A)')
671 0 : write (*, 1) trim(net_name)
672 : end if
673 2 : write (*, '(A)')
674 : end if
675 :
676 4 : if (.not. quiet) then
677 0 : write (*, 2) 'nfcn', nfcn
678 0 : write (*, 2) 'njac', njac
679 0 : write (*, 2) 'naccpt', naccpt
680 0 : write (*, 2) 'nrejct', nrejct
681 0 : write (*, '(A)')
682 0 : write (*, 2) 'number of steps', nstep
683 0 : write (*, '(A)')
684 0 : write (*, '(A)')
685 0 : write (*, 1) 'output file '//trim(data_filename)
686 0 : write (*, '(A)')
687 0 : write (*, '(/,a30,99f18.3,/)') 'runtime (seconds)', dt
688 0 : write (*, '(A)')
689 : end if
690 :
691 4 : if (save_final_abundances) then
692 0 : if (.not. quiet) write (*, *) 'save final abundances to '//trim(final_abundances_filename)
693 0 : open (newunit=iounit, file=trim(final_abundances_filename), iostat=ierr)
694 0 : if (ierr /= 0) then
695 0 : write (*, *) 'failed to open final_abundances_filename '//trim(final_abundances_filename)
696 : else
697 0 : write (iounit, 2) 'species', species
698 0 : do j = 1, species
699 0 : cid = chem_id(j)
700 0 : write (iounit, 3) trim(chem_isos%name(cid)), &
701 0 : chem_isos%Z(cid), &
702 0 : chem_isos%N(cid), &
703 0 : max(0d0, ending_x(j))
704 : end do
705 0 : close (iounit)
706 : end if
707 : end if
708 :
709 4 : if (associated(pm_work)) deallocate (pm_work)
710 0 : deallocate ( &
711 0 : rate_factors, rtol, atol, &
712 0 : ending_x, x_initial, x_previous, times, &
713 0 : log10Ts_f1, log10Rhos_f1, &
714 0 : etas_f1, log10Ps_f1, &
715 20 : peak_abundance, peak_time)
716 :
717 : contains
718 :
719 900 : subroutine get_eos_info_for_burn_at_const_density( &
720 0 : eos_handle, species, chem_id, net_iso, xa, &
721 : Rho, logRho, T, logT, &
722 : Cv, d_Cv_dlnT, eta, d_eta_dlnT, ierr)
723 4 : use eos_lib, only: eosDT_get
724 : use eos_def
725 : integer, intent(in) :: eos_handle, species
726 : integer, pointer :: chem_id(:) ! maps species to chem id
727 : integer, pointer :: net_iso(:) ! maps chem id to species number
728 : real(dp), intent(in) :: &
729 : xa(:), rho, logRho, T, logT
730 : real(dp), intent(out) :: &
731 : Cv, d_Cv_dlnT, eta, d_eta_dlnT
732 : integer, intent(out) :: ierr
733 :
734 0 : real(dp), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
735 0 : real(dp) :: d_dxa(num_eos_d_dxa_results, species)
736 :
737 : include 'formats'
738 0 : ierr = 0
739 :
740 0 : Cv = burn_Cv
741 0 : d_Cv_dlnT = burn_d_Cv_dlnT
742 0 : eta = burn_eta
743 0 : d_eta_dlnT = burn_deta_dlnT
744 :
745 0 : if (burn_Cv <= 0d0 .or. burn_eta <= 0d0) then
746 : call eosDT_get( &
747 0 : eos_handle, species, chem_id, net_iso, xa, &
748 : Rho, logRho, T, logT, &
749 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
750 :
751 0 : if (ierr /= 0) then
752 0 : write (*, *) 'failed in eosDT_get'
753 0 : return
754 : end if
755 0 : if (burn_Cv <= 0d0) then
756 0 : Cv = res(i_cv)
757 0 : d_Cv_dlnT = d_dlnT(i_cv)
758 : end if
759 0 : if (burn_eta <= 0d0) then
760 0 : eta = res(i_eta)
761 0 : d_eta_dlnT = d_dlnT(i_eta)
762 : end if
763 : end if
764 :
765 0 : end subroutine get_eos_info_for_burn_at_const_density
766 :
767 56 : subroutine burn_finish_substep(step, time, y, ierr)
768 : integer, intent(in) :: step
769 : real(dp), intent(in) :: time, y(:)
770 1232 : real(dp), dimension(size(y)) :: x
771 : integer, intent(out) :: ierr
772 : include 'formats'
773 56 : ierr = 0
774 : !if (burn_at_constant_density) then
775 : !write(*,2) 'finish_substep time xh logT T', &
776 : ! step, time, y(1), y(species+1)/ln10, exp(y(species+1))
777 : !return
778 : !end if
779 :
780 1232 : do i = 1, species
781 1176 : cid = chem_id(i)
782 1232 : x(i) = y(i)*dble(chem_isos%Z_plus_N(cid))
783 : end do
784 :
785 56 : if (burn_at_constant_density) then
786 0 : x(species + 1) = y(species + 1)
787 0 : logT = x(species + 1)/ln10
788 : end if
789 : call burn_solout1( &
790 56 : step, told, time, logT, logRho, species, x, ierr)
791 56 : told = time
792 0 : end subroutine burn_finish_substep
793 :
794 : real(dp) function interp_y(i, s, rwork_y, iwork_y, ierr)
795 : use const_def, only: dp
796 : integer, intent(in) :: i ! result is interpolated approximation of y(i) at x=s.
797 : real(dp), intent(in) :: s ! interpolation x value (between xold and x).
798 : real(dp), intent(inout), target :: rwork_y(*)
799 : integer, intent(inout), target :: iwork_y(*)
800 : integer, intent(out) :: ierr
801 : ierr = 0
802 : interp_y = 0
803 : end function interp_y
804 :
805 0 : subroutine do_read_T_Rho_history(ierr)
806 : integer, intent(out) :: ierr
807 : character(len=256) :: buffer, string
808 : integer :: i, n, iounit, t, num_isos, id, k
809 :
810 : include 'formats'
811 :
812 0 : ierr = 0
813 0 : write (*, *) 'read T Rho history from '//trim(T_Rho_history_filename)
814 : open (newunit=iounit, file=trim(T_Rho_history_filename), &
815 0 : action='read', status='old', iostat=ierr)
816 0 : if (ierr /= 0) then
817 0 : write (*, *) 'failed to open'
818 0 : return
819 : end if
820 :
821 0 : read (iounit, *, iostat=ierr) num_times
822 0 : if (ierr /= 0) then
823 0 : write (*, *) 'first line should have num_times'
824 0 : close (iounit)
825 0 : return
826 : end if
827 :
828 0 : if (num_times > max_num_times) then
829 0 : write (*, 3) 'num_times > max_num_times', num_times, max_num_times
830 0 : close (iounit)
831 0 : return
832 : end if
833 :
834 : if (.false.) write (*, 2) 'num_times', num_times
835 :
836 0 : do i = 1, num_times
837 0 : read (iounit, *, iostat=ierr) times_for_burn(i), log10Ts_for_burn(i), &
838 0 : log10Rhos_for_burn(i), etas_for_burn(i)
839 0 : if (ierr /= 0) then
840 0 : write (*, 2) 'failed reading line', i + 1
841 0 : write (*, '(a)') 'line should have values for time, logT, logRho, and eta'
842 0 : exit
843 : end if
844 : if (.false.) write (*, 2) 'history', i, times_for_burn(i), log10Ts_for_burn(i), &
845 0 : log10Rhos_for_burn(i), etas_for_burn(i)
846 : end do
847 :
848 0 : close (iounit)
849 :
850 0 : num_times_for_burn = num_times
851 :
852 : end subroutine do_read_T_Rho_history
853 :
854 0 : subroutine read_X(ierr)
855 : use utils_def
856 : use utils_lib
857 : integer, intent(out) :: ierr
858 : character(len=256) :: buffer, string
859 : integer :: i, n, iounit, t, num_isos, id, k
860 :
861 : include 'formats'
862 :
863 0 : write (*, *) 'read initial abundances from '//trim(initial_abundances_filename)
864 : open (newunit=iounit, file=trim(initial_abundances_filename), &
865 0 : action='read', status='old', iostat=ierr)
866 0 : if (ierr /= 0) then
867 0 : write (*, *) 'failed to open'
868 0 : return
869 : end if
870 :
871 0 : n = 0
872 0 : i = 0
873 0 : t = token(iounit, n, i, buffer, string)
874 0 : if (t /= name_token .or. string /= 'species') then
875 0 : write (*, *) 'expect to find specification of number of species at start of file'
876 0 : ierr = -1; return
877 : end if
878 0 : t = token(iounit, n, i, buffer, string)
879 0 : read (string, fmt=*, iostat=ierr) num_isos
880 0 : if (t /= name_token .or. ierr /= 0) then
881 0 : write (*, *) 'expect to find specification of number of species at start of file'
882 0 : ierr = -1; return
883 : end if
884 0 : if (num_isos /= species) then
885 0 : write (*, 2) 'expect to find number of species equal to those in current net', species
886 0 : ierr = -1; return
887 : end if
888 0 : do k = 1, species
889 0 : t = token(iounit, n, i, buffer, string)
890 0 : if (t /= name_token) then
891 0 : write (*, *) 'failed to find iso name at start of line: '//trim(string)
892 0 : ierr = -1; return
893 : end if
894 0 : id = get_nuclide_index(string)
895 0 : if (id <= 0) then
896 0 : write (*, *) 'failed to recognize name of iso '//trim(string)
897 0 : ierr = -1
898 0 : return
899 : end if
900 0 : j = net_iso(id)
901 0 : if (j <= 0 .or. j > species) then
902 0 : write (*, *) 'iso not in current net '//trim(string)
903 0 : ierr = 1
904 0 : return
905 : end if
906 0 : t = token(iounit, n, i, buffer, string)
907 0 : if (t /= name_token) then
908 : write (*, *) 'failed to read iso abundance: '// &
909 0 : trim(chem_isos%name(id))//' '//trim(string)
910 0 : ierr = -1; return
911 : end if
912 0 : read (string, fmt=*, iostat=ierr) xin(j)
913 0 : if (ierr /= 0) then
914 : write (*, *) 'failed to read iso abundance: ' &
915 0 : //trim(chem_isos%name(id))//' '//trim(string)
916 : end if
917 : end do
918 0 : close (iounit)
919 : end subroutine read_X
920 :
921 2 : subroutine show_X(X, show_peak, do_sort)
922 : use num_lib, only: qsort
923 : real(dp) :: X(:)
924 : logical, intent(in) :: show_peak, do_sort
925 44 : real(dp), target :: v_t(species)
926 4 : integer, target :: index_t(species)
927 2 : real(dp), pointer :: v(:)
928 2 : integer, pointer :: index(:)
929 : integer :: j
930 2 : real(dp) :: xsum
931 : include 'formats'
932 2 : v => v_t
933 2 : index => index_t
934 :
935 2 : if (do_sort) then
936 :
937 44 : v(1:species) = abs(x(1:species))
938 2 : call qsort(index, species, v)
939 44 : do i = 1, species
940 42 : if (i > max_num_for_show_peak_abundances) exit
941 42 : j = index(species + 1 - i)
942 42 : if (x(j) > min_for_show_peak_abundances) &
943 7 : write (*, 2) trim(chem_isos%name(chem_id(j))), i, x(j)
944 44 : if (x(j) > 1.1d0 .or. x(j) < -0.1d0) then
945 0 : write (*, 1) 'bad x for '//trim(chem_isos%name(chem_id(j))), x(j)
946 0 : call mesa_error(__FILE__, __LINE__)
947 : end if
948 : end do
949 :
950 : else
951 :
952 0 : do j = 1, species
953 0 : if (x(j) > min_for_show_peak_abundances) &
954 0 : write (*, 2) trim(chem_isos%name(chem_id(j))), j, x(j)
955 0 : if (x(j) > 1.1d0 .or. x(j) < -0.1d0) then
956 0 : write (*, 1) 'bad x for '//trim(chem_isos%name(chem_id(j))), x(j)
957 0 : call mesa_error(__FILE__, __LINE__)
958 : end if
959 :
960 : !write(*,1) 'xin(net_iso(i' // trim(chem_isos% name(chem_id(j))) // ')=', x(j)
961 :
962 : end do
963 :
964 : end if
965 : !stop
966 :
967 2 : write (*, '(A)')
968 44 : xsum = sum(x(1:species))
969 2 : write (*, 1) 'xsum', xsum
970 2 : write (*, '(A)')
971 2 : if (.not. show_peak) return
972 0 : write (*, '(A)')
973 0 : write (*, 1) 'peak x and time'
974 0 : do j = 1, species
975 0 : if (peak_abundance(j) >= min_for_show_peak_abundances) &
976 0 : write (*, 1) trim(chem_isos%name(chem_id(j))), &
977 0 : peak_abundance(j), peak_time(j)
978 : end do
979 0 : write (*, '(A)')
980 4 : end subroutine show_X
981 :
982 844 : subroutine burn_solout( &
983 844 : step, told, time, n, x, rwork_y, iwork_y, interp_y, &
984 : lrpar, rpar, lipar, ipar, irtrn)
985 2 : use const_def
986 : use eos_def
987 : use eos_lib
988 : integer, intent(in) :: step, n, lrpar, lipar
989 : real(dp), intent(in) :: told, time
990 : real(dp), intent(inout) :: x(:)
991 : real(dp), intent(inout), target :: rwork_y(*)
992 : integer, intent(inout), target :: iwork_y(*)
993 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
994 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
995 844 : real(dp) :: lgt, lgrho
996 : integer :: i, cid
997 : interface
998 : include 'num_interp_y.dek'
999 : end interface
1000 : integer, intent(out) :: irtrn ! < 0 causes solver to return to calling program.
1001 :
1002 844 : if (size(rpar) == burn_lrpar) then
1003 0 : lgt = rpar(r_burn_prev_lgT)
1004 0 : lgrho = rpar(r_burn_prev_lgRho)
1005 844 : else if (size(rpar) == burn_const_P_lrpar) then
1006 844 : lgt = log10(rpar(r_burn_const_P_temperature))
1007 844 : lgrho = log10(rpar(r_burn_const_P_rho))
1008 : end if
1009 :
1010 : call burn_solout1( &
1011 844 : step, told, time, lgt, lgrho, n, x, irtrn)
1012 844 : end subroutine burn_solout
1013 :
1014 900 : subroutine burn_solout1( &
1015 900 : step, told, time, logT_in, logRho_in, n, x, irtrn)
1016 844 : use const_def
1017 : use eos_def
1018 : use eos_lib
1019 : use chem_lib, only: get_Q
1020 : integer, intent(in) :: step, n
1021 : real(dp), intent(in) :: told, time, logT_in, logRho_in
1022 : real(dp), intent(in) :: x(:)
1023 : integer, intent(out) :: irtrn ! < 0 causes solver to return to calling program.
1024 :
1025 1800 : real(dp) :: logT, logRho, lgPgas, Pgas, Prad, lgP, avg_eps_nuc
1026 1800 : real(dp) :: eps_neu_total, eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, &
1027 22500 : eps_nuc_categories(num_categories)
1028 :
1029 900 : real(dp) :: dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas
1030 : real(dp) :: dlnRho_dlnT_const_P, d_epsnuc_dlnT_const_P, d_Cv_dlnT
1031 24300 : real(dp) :: res(num_eos_basic_results)
1032 24300 : real(dp) :: d_dlnRho_const_T(num_eos_basic_results)
1033 24300 : real(dp) :: d_dlnT_const_Rho(num_eos_basic_results)
1034 : real(dp) :: d_dabar_const_TRho(num_eos_basic_results)
1035 : real(dp) :: d_dzbar_const_TRho(num_eos_basic_results)
1036 57600 : real(dp) :: d_dxa_const_TRho(num_eos_d_dxa_results, species)
1037 :
1038 19800 : real(dp) :: Rho, T, xsum, d_eps_nuc_dx(species), dx, enuc, &
1039 1800 : dt, energy, entropy, burn_ergs, &
1040 900 : xh, xhe, Z, mass_correction
1041 :
1042 : integer :: i, j, adjustment_iso, cid, ierr, max_j
1043 57600 : real(dp), dimension(species) :: dabar_dx, dzbar_dx, eps_nuc_dx, dmc_dx
1044 900 : real(dp), pointer :: actual_Qs(:), actual_neuQs(:)
1045 900 : logical, pointer :: from_weaklib(:)
1046 : type(Net_General_Info), pointer :: g
1047 900 : type(net_info) :: netinfo
1048 : logical :: skip_jacobian
1049 :
1050 : include 'formats'
1051 :
1052 900 : irtrn = 0
1053 900 : if (time == 0) return
1054 :
1055 896 : logT = logT_in
1056 896 : logRho = logRho_in
1057 :
1058 896 : if ((.not. quiet) .and. step > 1 .and. mod(step, 50) == 0) then
1059 0 : max_j = maxloc(x(1:species), dim=1)
1060 0 : write (*, 2) 'step, time, logT, logRho, '//trim(chem_isos%name(chem_id(max_j))), &
1061 0 : step, time, logT, logRho, x(max_j)
1062 : if (.false.) then
1063 : do j = 1, species
1064 : write (*, 2) trim(chem_isos%name(chem_id(j))), j, x(j)
1065 : end do
1066 : write (*, 1) 'sum(x)', sum(x(1:species))
1067 : write (*, '(A)')
1068 : end if
1069 : end if
1070 :
1071 896 : ierr = 0
1072 :
1073 : allocate ( &
1074 : actual_Qs(num_reactions), actual_neuQs(num_reactions), &
1075 : from_weaklib(num_reactions), &
1076 896 : stat=ierr)
1077 896 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
1078 :
1079 19712 : xin(1:species) = x(1:species)
1080 :
1081 : call composition_info( &
1082 0 : species, chem_id, xin(1:species), xh, xhe, z, abar, zbar, z2bar, z53bar, ye, &
1083 896 : mass_correction, xsum, dabar_dx, dzbar_dx, dmc_dx)
1084 896 : Z = max(0d0, min(1d0, 1d0 - (xh + xhe)))
1085 :
1086 896 : if (burn_at_constant_P) then
1087 :
1088 842 : logT = x(n)/ln10
1089 842 : T = exp10(logT)
1090 842 : Prad = Radiation_Pressure(T)
1091 842 : Pgas = pressure - Prad
1092 842 : lgPgas = log10(Pgas)
1093 :
1094 : call eosPT_get( &
1095 : eos_handle, &
1096 0 : species, chem_id, net_iso, x, &
1097 : Pgas, lgPgas, T, logT, &
1098 : Rho, logRho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas, &
1099 : res, d_dlnRho_const_T, d_dlnT_const_Rho, &
1100 842 : d_dxa_const_TRho, ierr)
1101 :
1102 842 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
1103 :
1104 : else ! this is okay for burn_at_constant_density as well as constant T and Rho
1105 :
1106 : ! logT = rpar(r_burn_prev_lgT)
1107 : ! logRho = rpar(r_burn_prev_lgRho)
1108 54 : T = exp10(logT)
1109 54 : Rho = exp10(logRho)
1110 :
1111 : call eosDT_get( &
1112 : eos_handle, &
1113 0 : species, chem_id, net_iso, x, &
1114 : Rho, logRho, T, logT, &
1115 : res, d_dlnRho_const_T, d_dlnT_const_Rho, &
1116 54 : d_dxa_const_TRho, ierr)
1117 : !Pgas, Prad, energy, entropy, ierr)
1118 54 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
1119 :
1120 : end if
1121 :
1122 896 : lgPgas = res(i_lnPgas)/ln10
1123 896 : Pgas = exp10(lgPgas)
1124 896 : Prad = Radiation_Pressure(T)
1125 896 : lgP = log10(Pgas + Prad)
1126 896 : skip_jacobian = .false.
1127 896 : eta = res(i_eta)
1128 896 : d_eta_dlnT = 0d0
1129 896 : d_eta_dlnRho = 0d0
1130 :
1131 896 : call get_net_ptr(handle, g, ierr)
1132 896 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
1133 :
1134 896 : netinfo%g => g
1135 :
1136 : call net_get_with_Qs( &
1137 : handle, skip_jacobian, netinfo, species, num_reactions, &
1138 0 : xin(1:species), T, logT, Rho, logRho, &
1139 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
1140 : rate_factors, weak_rate_factor, &
1141 : std_reaction_Qs, std_reaction_neuQs, &
1142 0 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
1143 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
1144 : screening_mode, &
1145 : eps_nuc_categories, eps_neu_total, &
1146 896 : actual_Qs, actual_neuQs, from_weaklib, ierr)
1147 896 : if (ierr /= 0) then
1148 0 : write (*, '(A)')
1149 0 : write (*, 1) 'logT', logT
1150 0 : write (*, 1) 'logRho', logRho
1151 0 : write (*, *) 'bad return from net_get'
1152 0 : call mesa_error(__FILE__, __LINE__)
1153 : end if
1154 :
1155 896 : dt = time - told
1156 :
1157 : ! set burn_ergs according to change from initial abundances
1158 896 : eps_nuc = 0
1159 896 : xsum = 0
1160 896 : burn_ergs = 0
1161 19712 : do i = 1, species
1162 18816 : cid = chem_id(i)
1163 :
1164 18816 : dx = x(i) - x_initial(i)
1165 18816 : xsum = xsum + x(i)
1166 : burn_ergs = burn_ergs + &
1167 18816 : (get_Q(chem_isos, cid))*dx/chem_isos%Z_plus_N(cid)
1168 :
1169 18816 : dx = x(i) - x_previous(i)
1170 : eps_nuc = eps_nuc + &
1171 19712 : (get_Q(chem_isos, cid))*dx/chem_isos%Z_plus_N(cid)
1172 :
1173 : end do
1174 896 : avg_eps_nuc = burn_ergs*Qconv/time - eps_neu_total
1175 896 : eps_nuc = eps_nuc*Qconv/dt - eps_neu_total
1176 896 : burn_logT = logT
1177 896 : burn_logRho = logRho
1178 896 : burn_lnS = res(i_lnS)
1179 896 : burn_lnE = res(i_lnE)
1180 :
1181 19712 : x_previous(1:species) = x(1:species)
1182 :
1183 896 : if (time >= data_output_min_t) then
1184 :
1185 : write (io_out, '(i7,99(1pe26.16,1x))', advance='no') &
1186 896 : step, &
1187 896 : sign(1d0, avg_eps_nuc)*log10(max(1d0, abs(avg_eps_nuc))), &
1188 896 : avg_eps_nuc, &
1189 896 : sign(1d0, eps_nuc)*log10(max(1d0, abs(eps_nuc))), &
1190 896 : eps_nuc, &
1191 896 : eps_neu_total, &
1192 896 : burn_logT - starting_logT, &
1193 896 : burn_logT, &
1194 896 : burn_logRho, &
1195 896 : lgPgas, &
1196 896 : lgP, &
1197 896 : abar, &
1198 896 : zbar, &
1199 896 : exp(burn_lnS)/(avo*kerg), &
1200 896 : burn_lnE/ln10, &
1201 896 : time, safe_log10(time), safe_log10(time/secyer), &
1202 1792 : time - told, safe_log10(time - told), ye, xsum - 1
1203 14336 : do i = 1, num_names_of_isos_to_show
1204 13440 : if (num_names_of_isos_to_show < species) then
1205 13440 : cid = burn_isos_to_show(i)
1206 : else
1207 0 : if (i > species) exit
1208 0 : cid = chem_id(i)
1209 : end if
1210 13440 : j = net_iso(cid)
1211 13440 : if (j == 0) cycle
1212 :
1213 : ! output mass fractions, not abundances - fxt
1214 : ! write(io_out,'(1pe26.16,1x)',advance='no') safe_log10(x(j))
1215 14336 : write (io_out, '(1pe26.16,1x)', advance='no') safe_log10(x(j)*chem_isos%Z_plus_N(cid))
1216 : end do
1217 14336 : do i = 1, num_names_of_isos_to_show
1218 13440 : if (num_names_of_isos_to_show < species) then
1219 13440 : cid = burn_isos_to_show(i)
1220 : else
1221 0 : if (i > species) exit
1222 0 : cid = chem_id(i)
1223 : end if
1224 13440 : j = net_iso(cid)
1225 13440 : if (j == 0) cycle
1226 :
1227 : ! output mass fractions, not abundances - fxt
1228 : ! write(io_out,'(1pe26.16,1x)',advance='no') x(j)
1229 14336 : write (io_out, '(1pe26.16,1x)', advance='no') x(j)*chem_isos%Z_plus_N(cid)
1230 : end do
1231 896 : do i = 1, num_reactions_to_track
1232 0 : j = index_for_reaction_to_track(i)
1233 0 : if (j == 0) cycle
1234 0 : write (io_out, '(1pe26.16,1x)', advance='no') netinfo%rate_raw(j)
1235 896 : write (io_out, '(1pe26.16,1x)', advance='no') netinfo%rate_screened(j)
1236 : end do
1237 896 : write (io_out, *)
1238 19712 : do j = 1, species
1239 19712 : if (x(j) > peak_abundance(j)) then
1240 13276 : peak_abundance(j) = x(j)
1241 13276 : peak_time(j) = time
1242 : end if
1243 : end do
1244 : end if
1245 :
1246 896 : if (show_Qs) then
1247 0 : write (*, '(A)')
1248 0 : write (*, 1) 'logT', logT
1249 0 : write (*, 1) 'logRho', logRho
1250 0 : write (*, '(A)')
1251 0 : write (*, '(30x,4a20)') 'Q total', 'Q neutrino', 'Q total-neutrino'
1252 0 : do i = 1, num_reactions
1253 0 : if (from_weaklib(i)) then
1254 0 : write (*, '(a30,99f20.10)') 'weaklib '//trim(reaction_Name(reaction_id(i))), &
1255 0 : actual_Qs(i), actual_neuQs(i), actual_Qs(i) - actual_neuQs(i)
1256 : else
1257 0 : write (*, '(a30,99f20.10)') trim(reaction_Name(reaction_id(i))), &
1258 0 : actual_Qs(i), actual_neuQs(i), actual_Qs(i) - actual_neuQs(i)
1259 : end if
1260 : end do
1261 0 : write (*, '(A)')
1262 0 : call mesa_error(__FILE__, __LINE__, 'show_Qs')
1263 : end if
1264 :
1265 896 : deallocate (actual_Qs, actual_neuQs, from_weaklib)
1266 :
1267 2696 : end subroutine burn_solout1
1268 :
1269 : end subroutine Do_One_Zone_Burn
1270 :
1271 4 : subroutine test_net_setup(net_file_in)
1272 : character(len=*), intent(in) :: net_file_in
1273 : integer :: ierr, i
1274 :
1275 : include 'formats'
1276 :
1277 4 : net_file = net_file_in
1278 :
1279 4 : call net_init(ierr)
1280 4 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
1281 :
1282 4 : handle = alloc_net_handle(ierr)
1283 4 : if (ierr /= 0) then
1284 0 : write (*, *) 'alloc_net_handle failed'
1285 0 : call mesa_error(__FILE__, __LINE__)
1286 : end if
1287 :
1288 4 : call net_start_def(handle, ierr)
1289 4 : if (ierr /= 0) then
1290 0 : write (*, *) 'net_start_def failed'
1291 0 : call mesa_error(__FILE__, __LINE__)
1292 : end if
1293 :
1294 4 : call read_net_file(net_file, handle, ierr)
1295 4 : if (ierr /= 0) then
1296 0 : write (*, *) 'read_net_file failed ', trim(net_file)
1297 0 : call mesa_error(__FILE__, __LINE__)
1298 : end if
1299 :
1300 4 : call net_finish_def(handle, ierr)
1301 4 : if (ierr /= 0) then
1302 0 : write (*, *) 'net_finish_def failed'
1303 0 : call mesa_error(__FILE__, __LINE__)
1304 : end if
1305 :
1306 4 : call net_ptr(handle, g, ierr)
1307 4 : if (ierr /= 0) then
1308 0 : write (*, *) 'net_ptr failed'
1309 0 : call mesa_error(__FILE__, __LINE__)
1310 : end if
1311 :
1312 4 : species = g%num_isos
1313 4 : num_reactions = g%num_reactions
1314 :
1315 4 : allocate (reaction_id(num_reactions))
1316 :
1317 4 : call net_setup_tables(handle, cache_suffix, ierr)
1318 4 : if (ierr /= 0) then
1319 0 : write (*, *) 'net_setup_tables failed'
1320 0 : call mesa_error(__FILE__, __LINE__)
1321 : end if
1322 :
1323 4 : call get_chem_id_table(handle, species, chem_id, ierr)
1324 4 : if (ierr /= 0) then
1325 0 : write (*, *) 'get_chem_id_table failed'
1326 0 : call mesa_error(__FILE__, __LINE__)
1327 : end if
1328 :
1329 4 : call get_net_iso_table(handle, net_iso, ierr)
1330 4 : if (ierr /= 0) then
1331 0 : write (*, *) 'get_net_iso_table failed'
1332 0 : call mesa_error(__FILE__, __LINE__)
1333 : end if
1334 :
1335 4 : call get_reaction_id_table(handle, num_reactions, reaction_id, ierr)
1336 4 : if (ierr /= 0) then
1337 0 : write (*, *) 'get_reaction_id_table failed'
1338 0 : call mesa_error(__FILE__, __LINE__)
1339 : end if
1340 :
1341 : allocate ( &
1342 : xin(species), xin_copy(species), d_eps_nuc_dx(species), &
1343 4 : dxdt(species), d_dxdt_dRho(species), d_dxdt_dT(species), d_dxdt_dx(species, species))
1344 :
1345 28 : end subroutine test_net_setup
1346 :
1347 : end module mod_one_zone_support
1348 :
1349 : module mod_one_zone_burn
1350 : use chem_lib
1351 : use net_def
1352 : use net_lib
1353 : use rates_lib, only: rates_init
1354 : use rates_def
1355 : use const_def
1356 : use utils_lib
1357 : use mtx_def
1358 :
1359 : use mod_one_zone_support
1360 :
1361 : implicit none
1362 :
1363 : integer :: ierr, unit
1364 :
1365 : namelist /one_zone/ &
1366 : mesa_dir, net_name, quiet, show_ye_stuff, num_names_of_isos_to_show, names_of_isos_to_show, &
1367 : final_abundances_filename, save_final_abundances, show_peak_x_and_time, &
1368 : initial_abundances_filename, read_initial_abundances, &
1369 : read_T_Rho_history, T_Rho_history_filename, &
1370 : num_isos_for_Xinit, names_of_isos_for_Xinit, values_for_Xinit, uniform_Xinit, &
1371 : burn_tend, burn_rho, burn_temp, burn_logRho, burn_logT, burn_eta, burn_deta_dlnT, &
1372 : burn_Cv, burn_d_Cv_dlnT, &
1373 : eps, odescal, stptry, trace, burn_dbg, use_pivoting, &
1374 : burn_rtol, burn_atol, burn_xmin, burn_xmax, weak_rate_factor, &
1375 : min_for_show_peak_abundances, max_num_for_show_peak_abundances, &
1376 : data_output_min_t, data_filename, &
1377 : which_solver, screening_mode, &
1378 : data_heading_line, show_net_reactions_info, &
1379 : rattab_logT_lower_bound, rattab_logT_upper_bound, max_steps, max_step_size, &
1380 : decsol_switch, small_mtx_decsol, large_mtx_decsol, &
1381 : burn_at_constant_P, burn_at_constant_density, &
1382 : starting_temp, pressure, cache_suffix, &
1383 : num_times_for_burn, times_for_burn, log10Ts_for_burn, &
1384 : log10Rhos_for_burn, etas_for_burn, log10Ps_for_burn, &
1385 : show_Qs, num_reactions_to_track, reaction_to_track, &
1386 : num_special_rate_factors, reaction_for_special_factor, special_rate_factor
1387 :
1388 : contains
1389 :
1390 4 : subroutine do_one_burn(filename, qt)
1391 : character(len=*) :: filename
1392 : logical, intent(in) :: qt
1393 :
1394 : include 'formats'
1395 :
1396 : ! set defaults
1397 :
1398 4 : mesa_dir = '../..'
1399 4 : net_name = 'test.net'
1400 4 : quiet = .false.
1401 4 : show_ye_stuff = .false.
1402 4 : cache_suffix = '0'
1403 4 : final_abundances_filename = ''
1404 4 : save_final_abundances = .false.
1405 4 : initial_abundances_filename = ''
1406 4 : read_initial_abundances = .false.
1407 4 : read_T_Rho_history = .false.
1408 4 : T_Rho_history_filename = ''
1409 4 : burn_tend = 10 ! seconds
1410 4 : burn_rho = -1d99
1411 4 : burn_temp = -1d99
1412 4 : burn_logT = -1d99
1413 4 : burn_logRho = -1d99
1414 4 : burn_eta = -1d99
1415 4 : burn_deta_dlnT = -1d99
1416 4 : burn_Cv = -1d99
1417 4 : burn_d_Cv_dlnT = -1d99
1418 4 : burn_rtol = 1d-8
1419 4 : burn_atol = 1d-9
1420 4 : burn_xmin = -1d-10
1421 4 : burn_xmax = 1d0 + 1d-10
1422 4 : eps = 1d-8
1423 4 : odescal = 1d-12
1424 4 : stptry = 0d0
1425 4 : trace = .false.
1426 4 : burn_dbg = .false.
1427 4 : use_pivoting = .true.
1428 4 : show_net_reactions_info = .false.
1429 4 : show_Qs = .false.
1430 4 : decsol_switch = 50
1431 : ! if current number of species <= switch,
1432 : ! then use small_mtx_decsol,
1433 : ! else use large_mtx_decsol.
1434 4 : small_mtx_decsol = 'lapack'
1435 4 : large_mtx_decsol = ''
1436 4 : which_solver = 'rodas4_solver'
1437 4 : rattab_logT_lower_bound = -1
1438 4 : rattab_logT_upper_bound = -1
1439 4 : max_steps = 50000
1440 4 : uniform_Xinit = .false.
1441 4 : burn_at_constant_P = .false.
1442 4 : burn_at_constant_density = .false.
1443 4 : starting_temp = -1
1444 4 : pressure = -1
1445 4 : num_times_for_burn = 0 ! <= 0 means don't use the arrays
1446 4 : times_for_burn = 0
1447 4 : log10Ts_for_burn = 0
1448 4 : log10Rhos_for_burn = 0
1449 4 : etas_for_burn = 0
1450 4 : log10Ps_for_burn = 0
1451 4 : max_step_size = 0
1452 :
1453 4 : min_for_show_peak_abundances = 1d-3 ! show if peak is > this
1454 4 : max_num_for_show_peak_abundances = 21
1455 4 : show_peak_x_and_time = .true.
1456 :
1457 4 : data_filename = 'one_zone_burn.data'
1458 4 : data_output_min_t = -99
1459 :
1460 4 : num_names_of_isos_to_show = -1
1461 :
1462 4 : num_isos_for_Xinit = 4
1463 0 : names_of_isos_for_Xinit(1:num_isos_for_Xinit) = [ &
1464 20 : 'he4', 'c12', 'n14', 'o16']
1465 0 : values_for_Xinit(1:num_isos_for_Xinit) = [ &
1466 20 : 0.95d0, 0.005d0, 0.035d0, 0.010d0]
1467 :
1468 4 : screening_mode = extended_screening
1469 :
1470 4 : num_special_rate_factors = 0 ! must be <= max_num_special_rate_factors
1471 404 : reaction_for_special_factor(:) = ''
1472 404 : special_rate_factor(:) = 1
1473 :
1474 4 : num_reactions_to_track = 0
1475 404 : reaction_to_track(:) = ''
1476 :
1477 4 : weak_rate_factor = 1
1478 :
1479 : ! read inlist
1480 :
1481 4 : open (newunit=unit, file=trim(filename), action='read', delim='quote', iostat=ierr)
1482 4 : if (ierr /= 0) then
1483 0 : write (*, *) 'Failed to open control namelist file ', trim(filename)
1484 0 : call mesa_error(__FILE__, __LINE__)
1485 : else
1486 4 : read (unit, nml=one_zone, iostat=ierr)
1487 4 : close (unit)
1488 4 : if (ierr /= 0) then
1489 0 : write (*, *) 'Failed while trying to read control namelist file ', trim(filename)
1490 : write (*, '(a)') &
1491 0 : 'The following runtime error message might help you find the problem'
1492 0 : write (*, *)
1493 0 : open (newunit=unit, file=trim(filename), action='read', delim='quote', status='old', iostat=ierr)
1494 0 : read (unit, nml=one_zone)
1495 0 : call mesa_error(__FILE__, __LINE__)
1496 : end if
1497 : end if
1498 :
1499 : ! do initialization
1500 :
1501 4 : if (burn_temp < 0.d0 .and. burn_logT < 0.d0) then
1502 0 : call mesa_error(__FILE__, __LINE__, "Must set either burn_temp or burn_logT")
1503 0 : stop
1504 : end if
1505 4 : if (burn_temp < 0) burn_temp = exp10(burn_logT)
1506 4 : if (burn_logT < 0) burn_logT = log10(burn_temp)
1507 4 : if (burn_rho < 0.d0 .and. burn_logRho < 0.d0) then
1508 0 : call mesa_error(__FILE__, __LINE__, "Must set either burn_rho or burn_logRho")
1509 0 : stop
1510 : end if
1511 4 : if (burn_rho < 0) burn_rho = exp10(burn_logRho)
1512 4 : if (burn_logRho < 0) burn_logRho = log10(burn_rho)
1513 :
1514 4 : starting_temp = burn_temp
1515 :
1516 4 : allocate (net_iso(num_chem_isos), chem_id(num_chem_isos))
1517 :
1518 : !reaclib_filename = 'jina_reaclib_results_20130213default2'
1519 : !write(*,*) 'changing reaclib_filename'
1520 :
1521 4 : open (unit=io_out, file=trim(data_filename), action='write', iostat=ierr)
1522 4 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
1523 :
1524 4 : complete_silence_please = qt
1525 4 : call Do_One_Zone_Burn(net_name)
1526 :
1527 4 : open (unit=io_out, file=trim(data_filename), action='write', iostat=ierr)
1528 4 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
1529 :
1530 4 : end subroutine do_one_burn
1531 :
1532 : end module mod_one_zone_burn
|
