Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 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 struct_burn_mix
21 :
22 : use star_private_def
23 : use const_def, only: dp, i8, ln10, secyer, lsun
24 : use utils_lib, only: is_bad
25 :
26 : implicit none
27 :
28 : private
29 : public :: do_struct_burn_mix
30 :
31 : contains
32 :
33 22 : integer function do_struct_burn_mix(s, skip_global_corr_coeff_limit)
34 : use mix_info, only: get_convection_sigmas
35 : use rates_def, only: num_rvs
36 : use hydro_vars, only: set_vars_if_needed
37 : use star_utils, only: start_time, update_time
38 :
39 : type (star_info), pointer :: s
40 : logical, intent(in) :: skip_global_corr_coeff_limit
41 :
42 : integer :: nz, nvar, species, ierr, k
43 : integer(i8) :: time0
44 : logical :: do_chem
45 : real(dp) :: dt, tol_correction_norm, tol_max_correction, total
46 :
47 : include 'formats'
48 :
49 11 : ierr = 0
50 11 : s% non_epsnuc_energy_change_from_split_burn = 0d0
51 11 : dt = s% dt
52 :
53 11 : if (s% rsp_flag) then
54 0 : do_struct_burn_mix = do_rsp_step(s,dt)
55 : s% total_num_solver_iterations = &
56 0 : s% total_num_solver_iterations + s% num_solver_iterations
57 0 : s% total_num_solver_calls_made = s% total_num_solver_calls_made + 1
58 0 : if (do_struct_burn_mix == keep_going) &
59 : s% total_num_solver_calls_converged = &
60 0 : s% total_num_solver_calls_converged + 1
61 0 : return
62 : end if
63 :
64 11 : if (s% use_other_before_struct_burn_mix) then
65 0 : call s% other_before_struct_burn_mix(s% id, dt, do_struct_burn_mix)
66 0 : if (do_struct_burn_mix /= keep_going) return
67 : end if
68 :
69 11 : if (s% doing_timing) call start_time(s, time0, total)
70 :
71 11 : s% doing_struct_burn_mix = .true.
72 11 : nz = s% nz
73 :
74 11 : species = s% species
75 11 : s% num_solver_iterations = 0
76 :
77 11 : do_struct_burn_mix = retry
78 :
79 11 : s% do_burn = (s% dxdt_nuc_factor > 0d0)
80 11 : s% do_mix = (s% mix_factor > 0d0)
81 :
82 11 : if (s% op_split_burn) then
83 0 : do k=1,nz
84 0 : s% burn_num_iters(k) = 0
85 0 : if (s% T(k) >= s% op_split_burn_min_T) then
86 0 : s% eps_nuc(k) = 0d0
87 0 : s% d_epsnuc_dlnd(k) = 0d0
88 0 : s% d_epsnuc_dlnT(k) = 0d0
89 0 : s% d_epsnuc_dx(:,k) = 0d0
90 0 : s% eps_nuc_categories(:,k) = 0d0
91 0 : s% dxdt_nuc(:,k) = 0d0
92 0 : s% d_dxdt_nuc_dRho(:,k) = 0d0
93 0 : s% d_dxdt_nuc_dT(:,k) = 0d0
94 0 : s% d_dxdt_nuc_dx(:,:,k) = 0d0
95 0 : s% eps_nuc_neu_total(k) = 0d0
96 : end if
97 : end do
98 : end if
99 :
100 11 : if (s% do_burn .and. s% op_split_burn) then
101 0 : total = 0
102 0 : do k=1,s% nz
103 0 : total = total - s% energy(k)*s% dm(k)
104 : end do
105 0 : if (s% trace_evolve) write(*,*) 'call do_burn'
106 0 : do_struct_burn_mix = do_burn(s, dt)
107 : if (do_struct_burn_mix /= keep_going) then
108 : write(*,2) 'failed in do_burn', s% model_number
109 : call mesa_error(__FILE__,__LINE__,'do_struct_burn_mix')
110 : return
111 : end if
112 0 : call set_vars_if_needed(s, s% dt, 'after do_burn', ierr)
113 0 : if (ierr /= 0) return
114 0 : do k=1,s% nz
115 0 : total = total + s% energy(k)*s% dm(k)
116 : end do
117 0 : s% non_epsnuc_energy_change_from_split_burn = total
118 0 : if (s% trace_evolve) write(*,*) 'done do_burn'
119 : end if
120 :
121 11 : if (s% doing_first_model_of_run) then
122 1 : if (s% i_lum /= 0) then
123 1 : s% L_phot_old = s% xh(s% i_lum,1)/Lsun
124 : else
125 0 : s% L_phot_old = 0
126 : end if
127 : end if
128 :
129 11 : if (s% do_mix) then
130 11 : call get_convection_sigmas(s, dt, ierr)
131 11 : if (ierr /= 0) then
132 0 : if (s% report_ierr) write(*,*) 'get_convection_sigmas failed'
133 0 : return
134 : end if
135 : else
136 0 : s% sig(1:s% nz) = 0
137 : end if
138 :
139 11 : do_chem = (s% do_burn .or. s% do_mix)
140 : if (do_chem) then ! include abundances
141 11 : nvar = s% nvar_total
142 : else ! no chem => just do structure
143 0 : nvar = s% nvar_hydro
144 : end if
145 :
146 : call set_tol_correction(s, maxval(s% T(1:s% nz)), &
147 13098 : tol_correction_norm, tol_max_correction)
148 11 : call set_surf_info(s, nvar)
149 :
150 11 : if (s% w_div_wc_flag) then
151 0 : s% xh(s% i_w_div_wc,:s% nz) = s% w_div_w_crit_roche(:s% nz)
152 : end if
153 :
154 11 : if (s% j_rot_flag) then
155 0 : s% xh(s% i_j_rot,:s% nz) = s% j_rot(:s% nz)
156 0 : s% j_rot_start(:s% nz) = s% j_rot(:s% nz)
157 0 : do k=1, s% nz
158 0 : s% i_rot_start(k) = s% i_rot(k)% val
159 : end do
160 0 : s% total_abs_angular_momentum = dot_product(abs(s% j_rot(:s% nz)),s% dm_bar(:s% nz))
161 : end if
162 :
163 11 : call save_start_values(s, ierr)
164 11 : if (ierr /= 0) then
165 0 : if (s% report_ierr) write(*,*) 'save_start_values failed'
166 0 : return
167 : end if
168 :
169 11 : if (s% trace_evolve) write(*,*) 'call solver'
170 : do_struct_burn_mix = do_solver_converge( &
171 : s, nvar, skip_global_corr_coeff_limit, &
172 11 : tol_correction_norm, tol_max_correction)
173 11 : if (s% trace_evolve) write(*,*) 'done solver'
174 :
175 : s% total_num_solver_iterations = &
176 11 : s% total_num_solver_iterations + s% num_solver_iterations
177 11 : s% total_num_solver_calls_made = s% total_num_solver_calls_made + 1
178 11 : if (do_struct_burn_mix == keep_going) &
179 : s% total_num_solver_calls_converged = &
180 11 : s% total_num_solver_calls_converged + 1
181 11 : if (s% doing_relax) then
182 : s% total_num_solver_relax_iterations = &
183 0 : s% total_num_solver_relax_iterations + s% num_solver_iterations
184 0 : s% total_num_solver_relax_calls_made = s% total_num_solver_relax_calls_made + 1
185 0 : if (do_struct_burn_mix == keep_going) &
186 : s% total_num_solver_relax_calls_converged = &
187 0 : s% total_num_solver_relax_calls_converged + 1
188 : end if
189 :
190 11 : if (do_struct_burn_mix /= keep_going) return
191 :
192 11 : if (s% rotation_flag) then
193 : ! store ST mixing info for time smoothing
194 0 : do k=1, s% nz
195 0 : s% D_ST_start(k) = s% D_ST(k)
196 0 : s% nu_ST_start(k) = s% nu_ST(k)
197 : end do
198 0 : s% have_ST_start_info = .true.
199 : end if
200 :
201 11 : if (.not. s% j_rot_flag) &
202 11 : do_struct_burn_mix = do_mix_omega(s,dt)
203 :
204 11 : if (s% use_other_after_struct_burn_mix) &
205 0 : call s% other_after_struct_burn_mix(s% id, dt, do_struct_burn_mix)
206 :
207 11 : s% solver_iter = 0 ! to indicate that no longer doing solver iterations
208 11 : s% doing_struct_burn_mix = .false.
209 11 : if (s% doing_timing) call update_time(s, time0, total, s% time_struct_burn_mix)
210 :
211 : contains
212 :
213 : subroutine test(str)
214 : use chem_def, only: category_name
215 : character (len=*), intent(in) :: str
216 : include 'formats'
217 : integer :: k, i
218 : k = s% nz
219 : i = maxloc(s% eps_nuc_categories(1:num_categories,k),dim=1)
220 : write(*,3) trim(str) // ' eps_nuc_cat ' // trim(category_name(i)), &
221 : k, s% model_number, s% eps_nuc_categories(i,k)
222 : end subroutine test
223 :
224 : end function do_struct_burn_mix
225 :
226 :
227 11 : integer function do_mix_omega(s, dt)
228 : use solve_omega_mix, only: do_solve_omega_mix
229 : use hydro_rotation, only: set_i_rot, set_omega
230 :
231 : type (star_info), pointer :: s
232 : real(dp), intent(in) :: dt
233 :
234 11 : do_mix_omega = keep_going
235 :
236 11 : if (s% rotation_flag) then
237 : ! After solver is done with the structure, recompute moments of inertia and
238 : ! omega before angular momentum mix
239 0 : call set_i_rot(s, .false.)
240 0 : call set_omega(s, 'struct_burn_mix')
241 0 : if (s% premix_omega) then
242 0 : do_mix_omega = do_solve_omega_mix(s, 0.5d0*dt)
243 : else
244 0 : do_mix_omega = do_solve_omega_mix(s, dt)
245 : end if
246 0 : if (do_mix_omega /= keep_going) return
247 : end if
248 :
249 : end function do_mix_omega
250 :
251 :
252 0 : integer function do_rsp_step(s,dt)
253 : ! return keep_going, retry, or terminate
254 : use rsp, only: rsp_one_step
255 : type (star_info), pointer :: s
256 : real(dp), intent(in) :: dt
257 : integer :: ierr
258 0 : do_rsp_step = keep_going
259 : ierr = 0
260 0 : call rsp_one_step(s,ierr)
261 0 : if (ierr /= 0) then
262 0 : if (s% report_ierr) write(*,*) 'ierr from rsp_one_step'
263 : do_rsp_step = retry
264 : end if
265 0 : end function do_rsp_step
266 :
267 :
268 11 : subroutine save_start_values(s, ierr)
269 : use hydro_rsp2, only: set_etrb_start_vars
270 : use star_utils, only: eval_total_energy_integrals, set_luminosity_by_category, get_Peos_face_val
271 : type (star_info), pointer :: s
272 : integer, intent(out) :: ierr
273 : integer :: k, j
274 : include 'formats'
275 11 : ierr = 0
276 :
277 11 : call set_luminosity_by_category(s)
278 :
279 13098 : do k=1,s% nz
280 117783 : do j=1,s% species
281 117783 : s% dxdt_nuc_start(j,k) = s% dxdt_nuc(j,k)
282 : end do
283 327186 : do j=1,num_categories
284 : s% luminosity_by_category_start(j,k) = &
285 327175 : s% luminosity_by_category(j,k)
286 : end do
287 : end do
288 :
289 13098 : do k=1,s% nz
290 : !s% lnT_start(k) = s% lnT(k)
291 : !s% T_start(k) set elsewhere
292 : !s% lnd_start(k) = s% lnd(k) set elsewhere
293 : !s% rho_start(k) = s% rho(k)
294 : !s% r_start(k) set elsewhere
295 : !s% rmid_start(k) set elsewhere
296 : !s% v_start(k) set elsewhere
297 : !s% csound_start(k) set elsewhere
298 13087 : s% lnPeos_start(k) = s% lnPeos(k)
299 13087 : s% Peos_start(k) = s% Peos(k)
300 13087 : s% Peos_face_start(k) = get_Peos_face_val(s,k)
301 13087 : s% lnPgas_start(k) = s% lnPgas(k)
302 13087 : s% energy_start(k) = s% energy(k)
303 13087 : s% lnR_start(k) = s% lnR(k)
304 13087 : s% u_start(k) = s% u(k)
305 13087 : s% u_face_start(k) = 0d0 ! s% u_face_ad(k)%val
306 13087 : s% P_face_start(k) = -1d0 ! mark as unset s% P_face_ad(k)%val
307 13087 : s% L_start(k) = s% L(k)
308 13087 : s% Y_face_start(k) = s% Y_face(k)
309 13087 : s% omega_start(k) = s% omega(k)
310 13087 : s% ye_start(k) = s% ye(k)
311 13087 : s% j_rot_start(k) = s% j_rot(k)
312 13087 : s% eps_nuc_start(k) = s% eps_nuc(k)
313 13087 : s% non_nuc_neu_start(k) = s% non_nuc_neu(k)
314 13087 : s% Pvsc_start(k) = -1d99
315 13087 : s% grada_start(k) = s% grada(k)
316 13087 : s% chiT_start(k) = s% chiT(k)
317 13087 : s% chiRho_start(k) = s% chiRho(k)
318 13087 : s% cp_start(k) = s% cp(k)
319 13087 : s% Cv_start(k) = s% Cv(k)
320 13087 : s% dE_dRho_start(k) = s% dE_dRho(k)
321 13087 : s% gam_start(k) = s% gam(k)
322 13087 : s% lnS_start(k) = s% lnS(k)
323 13087 : s% zbar_start(k) = s% zbar(k)
324 13087 : s% mu_start(k) = s% mu(k)
325 13087 : s% phase_start(k) = s% phase(k)
326 13087 : s% latent_ddlnT_start(k) = s% latent_ddlnT(k)
327 13087 : s% latent_ddlnRho_start(k) = s% latent_ddlnRho(k)
328 13087 : s% eps_nuc_start(k) = s% eps_nuc(k)
329 13087 : s% opacity_start(k) = s% opacity(k)
330 13098 : s% m_grav_start(k) = s% m_grav(k)
331 : end do
332 :
333 11 : if (s% RSP2_flag) then
334 0 : call set_etrb_start_vars(s,ierr)
335 : end if
336 :
337 13098 : do k=1,s% nz
338 65446 : do j=1,s% nvar_hydro
339 65435 : s% xh_start(j,k) = s% xh(j,k)
340 : end do
341 : end do
342 :
343 13098 : do k=1,s% nz
344 117794 : do j=1,s% species
345 117783 : s% xa_start(j,k) = s% xa(j,k)
346 : end do
347 : end do
348 :
349 : call eval_total_energy_integrals(s, &
350 : s% total_internal_energy_start, &
351 : s% total_gravitational_energy_start, &
352 : s% total_radial_kinetic_energy_start, &
353 : s% total_rotational_kinetic_energy_start, &
354 : s% total_turbulent_energy_start, &
355 11 : s% total_energy_start)
356 :
357 11 : end subroutine save_start_values
358 :
359 :
360 11 : integer function do_solver_converge( &
361 : s, nvar, skip_global_corr_coeff_limit, &
362 : tol_correction_norm, tol_max_correction)
363 : ! return keep_going, retry, or terminate
364 : use mtx_lib
365 : use mtx_def
366 : use num_def
367 : use star_utils, only: start_time, update_time
368 :
369 : type (star_info), pointer :: s
370 : integer, intent(in) :: nvar
371 : logical, intent(in) :: skip_global_corr_coeff_limit
372 : real(dp), intent(in) :: tol_correction_norm, tol_max_correction
373 :
374 : integer :: ierr, nz, n
375 : logical :: report
376 :
377 : include 'formats'
378 :
379 11 : if (s% dt <= 0d0) then
380 11 : do_solver_converge = keep_going
381 : return
382 : end if
383 :
384 11 : do_solver_converge = terminate
385 :
386 11 : ierr = 0
387 :
388 11 : nz = s% nz
389 11 : n = nz*nvar
390 :
391 11 : s% solver_call_number = s% solver_call_number + 1
392 :
393 : do_solver_converge = do_solver( &
394 : s, skip_global_corr_coeff_limit, &
395 : tol_correction_norm, tol_max_correction, &
396 11 : report, nz, nvar)
397 :
398 :
399 11 : end function do_solver_converge
400 :
401 :
402 11 : subroutine set_surf_info(s, nvar) ! set to values at start of step
403 : type (star_info), pointer :: s
404 : integer, intent(in) :: nvar
405 11 : s% surf_lnS = s% lnS(1)
406 11 : s% num_surf_revisions = 0
407 11 : end subroutine set_surf_info
408 :
409 :
410 11 : subroutine set_xh(s,nvar,ierr) ! set xh using current structure info
411 : use hydro_rsp2, only: RSP2_adjust_vars_before_call_solver
412 : type (star_info), pointer :: s
413 : integer, intent(in) :: nvar
414 : integer, intent(out) :: ierr
415 : integer :: j1, k, nz
416 : include 'formats'
417 11 : ierr = 0
418 11 : nz = s%nz
419 11 : if (s% RSP2_flag) then
420 0 : call RSP2_adjust_vars_before_call_solver(s, ierr)
421 0 : if (ierr /= 0) then
422 0 : if (s% report_ierr) write(*,*) 'failed in RSP2_adjust_vars_before_call_solver'
423 0 : return
424 : end if
425 : end if
426 55 : do j1 = 1, min(nvar,s% nvar_hydro)
427 55 : if (j1 == s% i_lnd .and. s% i_lnd <= nvar) then
428 13098 : do k = 1, nz
429 13098 : s% xh(j1,k) = s% lnd(k)
430 : end do
431 33 : else if (j1 == s% i_lnT .and. s% i_lnT <= nvar) then
432 13098 : do k = 1, nz
433 13098 : s% xh(j1,k) = s% lnT(k)
434 : end do
435 22 : else if (j1 == s% i_lnR .and. s% i_lnR <= nvar) then
436 13098 : do k = 1, nz
437 13098 : s% xh(j1,k) = s% lnR(k)
438 : end do
439 11 : else if (j1 == s% i_lum .and. s% i_lum <= nvar) then
440 13098 : do k = 1, nz
441 13098 : s% xh(j1,k) = s% L(k)
442 : end do
443 0 : else if (j1 == s% i_w .and. s% i_w <= nvar) then
444 0 : do k = 1, nz
445 0 : s% xh(j1,k) = s% w(k)
446 : end do
447 0 : else if (j1 == s% i_Hp .and. s% i_Hp <= nvar) then
448 0 : do k = 1, nz
449 0 : s% xh(j1,k) = s% Hp_face(k)
450 : end do
451 0 : else if (j1 == s% i_v .and. s% i_v <= nvar) then
452 0 : do k = 1, nz
453 0 : s% xh(j1,k) = s% v(k)
454 : end do
455 0 : else if (j1 == s% i_u .and. s% i_u <= nvar) then
456 0 : do k = 1, nz
457 0 : s% xh(j1,k) = s% u(k)
458 : end do
459 0 : else if (j1 == s% i_alpha_RTI .and. s% i_alpha_RTI <= nvar) then
460 0 : do k = 1, nz
461 0 : s% xh(j1,k) = s% alpha_RTI(k)
462 : end do
463 : end if
464 : end do
465 : end subroutine set_xh
466 :
467 :
468 11 : subroutine set_tol_correction( &
469 : s, T_max, tol_correction_norm, tol_max_correction)
470 : type (star_info), pointer :: s
471 : real(dp), intent(in) :: T_max
472 : real(dp), intent(out) :: tol_correction_norm, tol_max_correction
473 : include 'formats'
474 11 : if (T_max >= s% tol_correction_extreme_T_limit) then
475 0 : tol_correction_norm = s% tol_correction_norm_extreme_T
476 0 : tol_max_correction = s% tol_max_correction_extreme_T
477 11 : else if (T_max >= s% tol_correction_high_T_limit) then
478 0 : tol_correction_norm = s% tol_correction_norm_high_T
479 0 : tol_max_correction = s% tol_max_correction_high_T
480 : else
481 11 : tol_correction_norm = s% tol_correction_norm
482 11 : tol_max_correction = s% tol_max_correction
483 : end if
484 11 : end subroutine set_tol_correction
485 :
486 :
487 33 : integer function do_solver( &
488 : s, skip_global_corr_coeff_limit, &
489 : tol_correction_norm, tol_max_correction, &
490 : report, nz, nvar)
491 : ! return keep_going, retry, or terminate
492 :
493 : ! when using solver for hydro step,
494 : ! do not require that functions have been evaluated for starting configuration.
495 : ! when finish, will have functions evaluated for the final set of primary variables.
496 : ! for example, the reaction rates will have been computed, so they can be used
497 : ! as initial values in the following burn and mix.
498 :
499 : use num_def
500 : use alloc
501 :
502 : type (star_info), pointer :: s
503 : integer, intent(in) :: nz, nvar
504 : logical, intent(in) :: skip_global_corr_coeff_limit, report
505 : real(dp), intent(in) :: tol_correction_norm, tol_max_correction
506 :
507 : logical :: converged
508 : integer :: k, species, ierr, j1, j2, gold_tolerances_level
509 : real(dp) :: maxT
510 :
511 : include 'formats'
512 :
513 11 : species = s% species
514 11 : do_solver = keep_going
515 11 : s% using_gold_tolerances = .false.
516 11 : gold_tolerances_level = 0
517 :
518 11 : if ((s% use_gold2_tolerances .and. s% steps_before_use_gold2_tolerances < 0) .or. &
519 : (s% steps_before_use_gold2_tolerances >= 0 .and. &
520 : s% model_number > s% steps_before_use_gold2_tolerances + max(0,s% init_model_number))) then
521 0 : s% using_gold_tolerances = .true.
522 0 : gold_tolerances_level = 2
523 11 : else if ((s% use_gold_tolerances .and. s% steps_before_use_gold_tolerances < 0) .or. &
524 : (s% steps_before_use_gold_tolerances >= 0 .and. &
525 : s% model_number > s% steps_before_use_gold_tolerances + max(0,s% init_model_number))) then
526 11 : if (s% maxT_for_gold_tolerances > 0) then
527 13098 : maxT = maxval(s% T(1:nz))
528 : else
529 : maxT = -1d0
530 : end if
531 11 : if (maxT > s% maxT_for_gold_tolerances) then
532 : !write(*,2) 'exceed maxT_for_gold_tolerances', &
533 : ! s% model_number, maxT, s% maxT_for_gold_tolerances
534 : else ! okay for maxT, so check if also ok for eosPC_frac
535 11 : s% using_gold_tolerances = .true.
536 11 : gold_tolerances_level = 1
537 : end if
538 : end if
539 :
540 11 : call set_xh(s, nvar, ierr) ! set xh using current structure info
541 11 : if (ierr /= 0) then
542 0 : if (report) then
543 0 : write(*, *) 'set_xh returned ierr in struct_burn_mix', ierr
544 0 : write(*, *) 's% model_number', s% model_number
545 0 : write(*, *) 'nz', nz
546 0 : write(*, *) 's% num_retries', s% num_retries
547 0 : write(*, *)
548 : end if
549 0 : do_solver = retry
550 0 : s% result_reason = nonzero_ierr
551 : s% dt_why_retry_count(Tlim_solver) = &
552 0 : s% dt_why_retry_count(Tlim_solver) + 1
553 0 : return
554 : end if
555 :
556 13098 : do k = 1, nz
557 65446 : do j1 = 1, min(nvar, s% nvar_hydro)
558 65435 : s% solver_dx(j1,k) = s% xh(j1,k) - s% xh_start(j1,k)
559 : end do
560 : end do
561 :
562 11 : if (nvar >= s% nvar_hydro+1) then
563 13098 : do k = 1, nz
564 13087 : j2 = 1
565 117794 : do j1 = s% nvar_hydro+1, nvar
566 104696 : s% xa_sub_xa_start(j2,k) = s% xa(j2,k) - s% xa_start(j2,k)
567 104696 : s% solver_dx(j1,k) = s% xa_sub_xa_start(j2,k)
568 117783 : j2 = j2+1
569 : end do
570 : end do
571 : end if
572 :
573 : converged = .false.
574 : call hydro_solver_step( &
575 : s, nz, s% nvar_hydro, nvar, skip_global_corr_coeff_limit, &
576 : gold_tolerances_level, tol_max_correction, tol_correction_norm, &
577 11 : converged, ierr)
578 11 : if (ierr /= 0) then
579 0 : if (report) then
580 0 : write(*, *) 'hydro_solver_step returned ierr', ierr
581 0 : write(*, *) 's% model_number', s% model_number
582 0 : write(*, *) 'nz', nz
583 0 : write(*, *) 's% num_retries', s% num_retries
584 0 : write(*, *)
585 : end if
586 0 : do_solver = retry
587 0 : s% result_reason = nonzero_ierr
588 : s% dt_why_retry_count(Tlim_solver) = &
589 0 : s% dt_why_retry_count(Tlim_solver) + 1
590 0 : return
591 : end if
592 :
593 11 : if (converged) then ! sanity checks before accept it
594 11 : converged = check_after_converge(s, report, ierr)
595 11 : if (converged .and. s% RTI_flag) & ! special checks
596 0 : converged = RTI_check_after_converge(s, report, ierr)
597 : end if
598 :
599 11 : if (.not. converged) then
600 0 : do_solver = retry
601 0 : s% result_reason = hydro_failed_to_converge
602 : s% dt_why_retry_count(Tlim_solver) = &
603 0 : s% dt_why_retry_count(Tlim_solver) + 1
604 0 : if (report) then
605 0 : write(*,2) 'solver rejected trial model'
606 0 : write(*,2) 's% model_number', s% model_number
607 0 : write(*,2) 's% solver_call_number', s% solver_call_number
608 0 : write(*,2) 'nz', nz
609 0 : write(*,2) 's% num_retries', s% num_retries
610 0 : write(*,1) 'dt', s% dt
611 0 : write(*,1) 'log dt/secyer', log10(s% dt/secyer)
612 0 : write(*, *)
613 : end if
614 0 : return
615 : end if
616 :
617 : end function do_solver
618 :
619 :
620 0 : logical function RTI_check_after_converge(s, report, ierr) result(converged)
621 : use mesh_adjust, only: set_lnT_for_energy
622 : use micro, only: do_eos_for_cell
623 : use chem_def, only: ih1, ihe3, ihe4
624 : use star_utils, only: store_lnT_in_xh, get_T_and_lnT_from_xh
625 : type (star_info), pointer :: s
626 : logical, intent(in) :: report
627 : integer, intent(out) :: ierr
628 : integer :: k, nz
629 : real(dp) :: old_energy, old_IE, new_IE, old_KE, new_KE, new_u, new_v, &
630 : revised_energy, new_lnT
631 : include 'formats'
632 0 : ierr = 0
633 0 : nz = s% nz
634 0 : converged = .true.
635 : !return
636 0 : do k=1,nz
637 0 : if (k < nz .and. s% alpha_RTI(k) < 1d-10) cycle
638 0 : old_energy = s% energy(k)
639 0 : old_IE = old_energy*s% dm(k)
640 0 : if (s% energy(k) < s% RTI_energy_floor) then
641 : ! try to take from KE to give to IE
642 : ! else just bump energy and take hit to energy conservation
643 0 : s% energy(k) = s% RTI_energy_floor
644 0 : s% lnE(k) = log(s% energy(k))
645 : call set_lnT_for_energy(s, k, &
646 : s% net_iso(ih1), s% net_iso(ihe3), s% net_iso(ihe4), &
647 : s% species, s% xa(:,k), &
648 : s% rho(k), s% lnd(k)/ln10, s% energy(k), s% lnT(k), &
649 0 : new_lnT, revised_energy, ierr)
650 0 : if (ierr /= 0) return ! call mesa_error(__FILE__,__LINE__,'do_merge failed in set_lnT_for_energy')
651 0 : call store_lnT_in_xh(s, k, new_lnT)
652 0 : call get_T_and_lnT_from_xh(s, k, s% T(k), s% lnT(k))
653 : end if
654 0 : new_IE = s% energy(k)*s% dm(k)
655 0 : if (s% u_flag) then
656 0 : old_KE = 0.5d0*s% dm(k)*s% u(k)*s% u(k)
657 0 : new_KE = max(0d0, old_KE + old_IE - new_IE)
658 0 : new_u = sqrt(new_KE/(0.5d0*s% dm(k)))
659 0 : if (s% u(k) > 0d0) then
660 0 : s% u(k) = new_u
661 : else
662 0 : s% u(k) = -new_u
663 : end if
664 0 : s% xh(s% i_u, k) = s% u(k)
665 0 : else if (s% v_flag) then ! only rough approximation possible here
666 0 : old_KE = 0.5d0*s% dm_bar(k)*s% v(k)*s% v(k)
667 0 : new_KE = max(0d0, old_KE + old_IE - new_IE)
668 0 : new_v = sqrt(max(0d0,new_KE)/(0.5d0*s% dm_bar(k)))
669 0 : if (s% v(k) > 0d0) then
670 0 : s% v(k) = new_v
671 : else
672 0 : s% v(k) = -new_v
673 : end if
674 0 : s% xh(s% i_v, k) = s% v(k)
675 : end if
676 : end do
677 : end function RTI_check_after_converge
678 :
679 :
680 11 : logical function check_after_converge(s, report, ierr) result(converged)
681 : type (star_info), pointer :: s
682 : logical, intent(in) :: report
683 : integer, intent(out) :: ierr
684 : integer :: k, nz
685 : include 'formats'
686 11 : ierr = 0
687 11 : nz = s% nz
688 11 : converged = .true.
689 11 : if (s% R_center > 0) then
690 0 : if (s% R_center > exp(s% lnR(nz))) then
691 0 : if (report) &
692 0 : write(*,2) 'volume < 0 in cell nz', nz, &
693 0 : s% R_center - exp(s% lnR(nz)), s% R_center, exp(s% lnR(nz)), &
694 0 : s% dm(nz), s% rho(nz), s% dq(nz)
695 0 : converged = .false.
696 0 : return
697 : end if
698 : end if
699 13087 : do k=1,nz-1
700 13087 : if (s% lnR(k) <= s% lnR(k+1)) then
701 0 : if (report) write(*,2) 'after hydro, negative cell volume in cell k', &
702 0 : k, s% lnR(k) - s% lnR(k+1), s% lnR(k), s% lnR(k+1), &
703 0 : s% lnR_start(k) - s% lnR_start(k+1), s% lnR_start(k), s% lnR_start(k+1)
704 : converged = .false.; exit
705 : call mesa_error(__FILE__,__LINE__,'check_after_converge')
706 : else
707 13076 : if (s% lnT(k) > ln10*12) then
708 0 : if (report) write(*,2) 'after hydro, logT > 12 in cell k', k, s% lnT(k)
709 : converged = .false. !; exit
710 13076 : else if (s% lnT(k) < ln10) then
711 0 : if (report) write(*,*) 'after hydro, logT < 1 in cell k', k
712 : converged = .false. !; exit
713 13076 : else if (s% lnd(k) > ln10*12) then
714 0 : if (report) write(*,*) 'after hydro, logRho > 12 in cell k', k
715 : converged = .false. !; exit
716 13076 : else if (s% lnd(k) < -ln10*30) then
717 0 : if (report) write(*,*) 'after hydro, logRho < -30 in cell k', k
718 : converged = .false. !; exit
719 : end if
720 : end if
721 : end do
722 : end function check_after_converge
723 :
724 :
725 22 : subroutine hydro_solver_step( &
726 : s, nz, nvar_hydro, nvar, skip_global_corr_coeff_limit, &
727 : gold_tolerances_level, tol_max_correction, tol_correction_norm, &
728 : converged, ierr)
729 : use num_def
730 : use chem_def
731 : use mtx_lib
732 : use mtx_def
733 : use alloc
734 :
735 : type (star_info), pointer :: s
736 : integer, intent(in) :: nz, nvar_hydro, nvar
737 : logical, intent(in) :: skip_global_corr_coeff_limit
738 : real(dp), intent(in) :: tol_max_correction, tol_correction_norm
739 : integer, intent(in) :: gold_tolerances_level
740 : logical, intent(out) :: converged
741 : integer, intent(out) :: ierr
742 :
743 : integer :: k, j, neq
744 : logical :: failure
745 : logical, parameter :: dbg = .false.
746 :
747 : include 'formats'
748 :
749 : ierr = 0
750 :
751 11 : neq = nvar*nz
752 :
753 : if (dbg) write(*, *) 'enter hydro_solver_step'
754 :
755 11 : s% used_extra_iter_in_solver_for_accretion = .false.
756 :
757 11 : call check_sizes(s, ierr)
758 11 : if (ierr /= 0) then
759 0 : write(*,*) 'check_sizes failed'
760 : return
761 : end if
762 :
763 : if (dbg) write(*, *) 'call solver'
764 11 : call newt(ierr)
765 11 : if (ierr /= 0 .and. s% report_ierr) then
766 0 : write(*,*) 'solver failed for hydro'
767 : end if
768 :
769 11 : converged = (ierr == 0) .and. (.not. failure)
770 11 : if (converged) then
771 13098 : do k=1,nz
772 65446 : do j=1,min(nvar,nvar_hydro)
773 65435 : s% xh(j,k) = s% xh_start(j,k) + s% solver_dx(j,k)
774 : end do
775 : end do
776 : ! s% xa has already been updated by final call to set_solver_vars from solver
777 : end if
778 :
779 :
780 : contains
781 :
782 :
783 11 : subroutine newt(ierr)
784 : use star_solver, only: solver
785 : use rates_def, only: warn_rates_for_high_temp
786 : integer, intent(out) :: ierr
787 : logical :: save_warn_rates_flag
788 : include 'formats'
789 11 : s% doing_solver_iterations = .true.
790 11 : save_warn_rates_flag = warn_rates_for_high_temp
791 11 : warn_rates_for_high_temp = .false.
792 : call solver( &
793 : s, nvar, skip_global_corr_coeff_limit, &
794 : gold_tolerances_level, tol_max_correction, tol_correction_norm, &
795 11 : failure, ierr)
796 11 : s% doing_solver_iterations = .false.
797 11 : warn_rates_for_high_temp = save_warn_rates_flag
798 11 : end subroutine newt
799 :
800 :
801 : end subroutine hydro_solver_step
802 :
803 :
804 0 : integer function do_burn(s, dt)
805 : use star_utils, only: start_time, update_time
806 : use net, only: get_screening_mode
807 : use chem_def
808 : use micro, only: do_eos_for_cell
809 :
810 : type (star_info), pointer :: s
811 : real(dp), intent(in) :: dt
812 :
813 : integer :: &
814 : k_bad, ierr, max_num_iters_k, nz, op_err, &
815 : k, num_iters, species, max_num_iters_used, &
816 : screening_mode, kmin
817 : integer(i8) :: time0
818 : real(dp) :: total, avg_epsnuc, min_T_for_const_density_solver
819 : logical :: trace, dbg, skip_burn
820 : logical, parameter :: burn_dbg = .false.
821 :
822 : include 'formats'
823 :
824 0 : trace = .false.
825 :
826 0 : min_T_for_const_density_solver = s% op_split_burn_min_T_for_variable_T_solver
827 :
828 0 : do_burn = keep_going
829 : ierr = 0
830 0 : nz = s% nz
831 0 : species = s% species
832 :
833 0 : if (s% eps_nuc_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) then
834 0 : do k = 1, nz
835 0 : s% eps_nuc(k) = 0d0
836 0 : s% burn_num_iters(k) = 0
837 0 : s% burn_avg_epsnuc(k) = 0d0
838 0 : s% max_burn_correction(k) = 0d0
839 : end do
840 : return
841 : end if
842 :
843 0 : if (dt <= 0d0) return
844 :
845 0 : max_num_iters_used = 0
846 0 : max_num_iters_k = 0
847 0 : k_bad = 0
848 :
849 0 : screening_mode = get_screening_mode(s,ierr)
850 0 : if (ierr /= 0) then
851 0 : if (s% report_ierr) &
852 0 : write(*,*) 'unknown string for screening_mode: ' // trim(s% screening_mode)
853 0 : return
854 : call mesa_error(__FILE__,__LINE__,'do1_net')
855 : end if
856 :
857 0 : dbg = .false. ! (s% model_number == 1137)
858 :
859 0 : kmin = nz+1
860 0 : do k=1,nz
861 0 : if (s% T_start(k) < s% op_split_burn_min_T) then
862 : ! We get here if we have an off center ignition,
863 : ! the arrays wont have been initialised earlier as they stop at the
864 : ! first temperature that exceeds op_split_burn_min_T
865 0 : s% burn_num_iters(k) = 0
866 0 : s% burn_avg_epsnuc(k) = 0d0
867 : cycle
868 : end if
869 : kmin = k
870 0 : exit
871 : end do
872 :
873 0 : if (kmin > nz) return
874 :
875 : !skip_burn = s% fe_core_infall > s% op_split_burn_eps_nuc_infall_limit
876 0 : skip_burn = (minval(s% v_start(1:s% nz)) < -s% op_split_burn_eps_nuc_infall_limit)
877 :
878 0 : if (s% doing_timing) call start_time(s, time0, total)
879 :
880 0 : !$OMP PARALLEL DO PRIVATE(k,op_err,num_iters,avg_epsnuc) SCHEDULE(dynamic,2)
881 : do k = kmin, nz
882 : if (k_bad /= 0) cycle
883 : if (s% T_start(k) < s% op_split_burn_min_T) then
884 : ! We get here if we have an off center ignition,
885 : ! the arrays wont have been initialised earlier as they stop at the
886 : ! first temperature that exceeds op_split_burn_min_T
887 : s% burn_num_iters(k) = 0
888 : s% burn_avg_epsnuc(k) = 0d0
889 : cycle
890 : end if
891 : s% max_burn_correction(k) = 0d0
892 : op_err = 0
893 : call burn1_zone( &
894 : s, k, species, min_T_for_const_density_solver, skip_burn, &
895 : screening_mode, &
896 : dt, num_iters, avg_epsnuc, burn_dbg, op_err)
897 : if (op_err /= 0) then
898 : ierr = -1
899 : k_bad = k
900 : cycle
901 : end if
902 : call do_eos_for_cell(s,k,op_err)
903 : if (op_err /= 0) then
904 : write(*,2) 'do_burn failed in do_eos_for_cell', k
905 : ierr = -1
906 : k_bad = k
907 : cycle
908 : end if
909 : !write(*,3) 'num_iters', k, num_iters
910 : s% burn_num_iters(k) = num_iters
911 : s% burn_avg_epsnuc(k) = avg_epsnuc
912 : if (num_iters > max_num_iters_used) then
913 : max_num_iters_used = num_iters
914 : max_num_iters_k = k
915 : end if
916 : end do
917 : !$OMP END PARALLEL DO
918 :
919 0 : s% need_to_setvars = .true.
920 :
921 0 : if (s% doing_timing) &
922 0 : call update_time(s, time0, total, s% time_solve_burn)
923 :
924 0 : if (ierr /= 0) then
925 0 : if (s% report_ierr) write(*,2) 'do_burn failed', k_bad
926 0 : return
927 : call mesa_error(__FILE__,__LINE__,'do_burn')
928 :
929 :
930 : do_burn = retry
931 : if (trace .or. s% report_ierr) then
932 : write(*,*) 'do_burn ierr'
933 : end if
934 : call restore
935 : return
936 : end if
937 :
938 : if (dbg) write(*,2) 'done do_burn'
939 :
940 :
941 : contains
942 :
943 : subroutine restore
944 : integer :: j, k
945 : do k = 1, nz
946 : do j=1,species
947 : s% xa(j,k) = s% xa_start(j,k)
948 : end do
949 : end do
950 : end subroutine restore
951 :
952 : end function do_burn
953 :
954 :
955 0 : subroutine burn1_zone( &
956 : s, k, species, min_T_for_const_density_solver, skip_burn, &
957 : screening_mode, &
958 : dt, num_iters_out, avg_epsnuc, dbg_in, ierr)
959 : use net_lib, only: net_1_zone_burn_const_density, net_1_zone_burn, &
960 : show_net_reactions_and_info
961 : use rates_def, only: std_reaction_Qs, std_reaction_neuQs
962 : use chem_def, only: num_categories
963 : use net, only: do1_net
964 : use star_utils, only: store_lnT_in_xh, get_T_and_lnT_from_xh
965 : type (star_info), pointer :: s
966 : integer, intent(in) :: k, species, screening_mode
967 : real(dp), intent(in) :: dt, min_T_for_const_density_solver
968 : logical, intent(in) :: skip_burn, dbg_in
969 : real(dp), intent(out) :: avg_epsnuc
970 : integer, intent(out) :: num_iters_out, ierr
971 :
972 0 : real(dp), target :: xa_start_ary(species)
973 0 : real(dp), pointer :: xa_start(:)
974 :
975 : real(dp) :: stptry, eps, odescal, &
976 : starting_log10T, ending_log10T, ending_eps_neu_total, &
977 : Cv0, eta0, substep_start_time
978 : integer :: i, max_steps, nfcn, njac, ntry, naccpt, nrejct
979 : integer, parameter :: num_times = 1
980 : real(dp), target, dimension(4*num_times) :: log10Ts_ary, log10Rhos_ary, etas_ary
981 0 : real(dp), pointer, dimension(:) :: log10Ts_f1, log10Rhos_f1, etas_f1, &
982 0 : dxdt_source_term, times
983 : logical :: use_pivoting, trace, burn_dbg
984 :
985 : include 'formats'
986 :
987 0 : ierr = 0
988 0 : num_iters_out = 0
989 :
990 0 : if (skip_burn) then
991 0 : avg_epsnuc = 0d0
992 0 : s% eps_nuc(k) = 0d0
993 0 : s% d_epsnuc_dlnd(k) = 0d0
994 0 : s% d_epsnuc_dlnT(k) = 0d0
995 0 : s% d_epsnuc_dx(:,k) = 0d0
996 0 : s% dxdt_nuc(:,k) = 0d0
997 0 : s% eps_nuc_categories(:,k) = 0d0
998 0 : s% d_dxdt_nuc_dRho(:,k) = 0d0
999 0 : s% d_dxdt_nuc_dT(:,k) = 0d0
1000 0 : s% d_dxdt_nuc_dx(:,:,k) = 0d0
1001 0 : s% eps_nuc_neu_total(k) = 0d0
1002 0 : return
1003 : end if
1004 :
1005 0 : log10Ts_f1 => log10Ts_ary
1006 0 : log10Rhos_f1 => log10Rhos_ary
1007 0 : etas_f1 => etas_ary
1008 :
1009 0 : nullify(dxdt_source_term, times)
1010 :
1011 0 : xa_start => xa_start_ary
1012 :
1013 0 : stptry = 0d0
1014 0 : eps = s% op_split_burn_eps
1015 0 : odescal = s% op_split_burn_odescal
1016 0 : max_steps = s% burn_steps_hard_limit
1017 0 : use_pivoting = .false. ! .true.
1018 0 : trace = .false.
1019 0 : burn_dbg = .false.
1020 0 : starting_log10T = s% lnT(k)/ln10
1021 :
1022 0 : do i=1,species
1023 0 : xa_start(i) = s% xa(i,k)
1024 : end do
1025 :
1026 0 : substep_start_time = 0d0
1027 :
1028 0 : if (s% use_other_split_burn) then
1029 0 : log10Ts_f1 => log10Ts_ary
1030 0 : log10Rhos_f1 => log10Rhos_ary
1031 0 : etas_f1 => etas_ary
1032 : nullify(dxdt_source_term, times)
1033 0 : log10Ts_f1(1) = s% lnT(k)/ln10
1034 0 : log10Rhos_f1(1) = s% lnd(k)/ln10
1035 0 : etas_f1(1) = s% eta(k)
1036 : call s% other_split_burn( &
1037 : s%id, k, s% net_handle, s% eos_handle, species, s% num_reactions, 0d0, dt, xa_start, &
1038 : num_times, times, log10Ts_f1, log10Rhos_f1, etas_f1, dxdt_source_term, &
1039 : s% rate_factors, s% weak_rate_factor, &
1040 : std_reaction_Qs, std_reaction_neuQs, &
1041 : screening_mode, &
1042 : stptry, max_steps, eps, odescal, &
1043 : use_pivoting, trace, burn_dbg, burn_finish_substep, &
1044 : s% xa(1:species,k), &
1045 : s% eps_nuc_categories(:,k), &
1046 : avg_epsnuc, ending_eps_neu_total, &
1047 0 : nfcn, njac, ntry, naccpt, nrejct, ierr)
1048 0 : if (ierr /= 0) then
1049 0 : if (s% report_ierr) write(*,2) 'other_split_burn failed', k
1050 0 : return
1051 : call mesa_error(__FILE__,__LINE__,'burn1_zone')
1052 : end if
1053 :
1054 0 : else if (s% T(k) >= min_T_for_const_density_solver) then
1055 0 : Cv0 = s% Cv(k)
1056 0 : eta0 = s% eta(k)
1057 : call net_1_zone_burn_const_density( &
1058 : s% net_handle, s% eos_handle, species, s% num_reactions, 0d0, dt, &
1059 : xa_start, starting_log10T, s% lnd(k)/ln10, &
1060 : get_eos_info_for_burn_at_const_density, &
1061 : s% rate_factors, s% weak_rate_factor, &
1062 : std_reaction_Qs, std_reaction_neuQs, &
1063 : screening_mode, &
1064 : stptry, max_steps, eps, odescal, &
1065 : use_pivoting, trace, burn_dbg, burn_finish_substep, &
1066 : s% xa(1:species,k), &
1067 : s% eps_nuc_categories(:,k), &
1068 : ending_log10T, avg_epsnuc, ending_eps_neu_total, &
1069 0 : nfcn, njac, ntry, naccpt, nrejct, ierr)
1070 0 : if (ierr /= 0) then
1071 0 : if (s% report_ierr) write(*,2) 'net_1_zone_burn_const_density failed', k
1072 0 : return
1073 : call mesa_error(__FILE__,__LINE__,'burn1_zone')
1074 : end if
1075 : ! restore temperature
1076 0 : call store_lnT_in_xh(s, k, starting_log10T*ln10)
1077 0 : call get_T_and_lnT_from_xh(s, k, s% T(k), s% lnT(k))
1078 : else
1079 0 : log10Ts_f1 => log10Ts_ary
1080 0 : log10Rhos_f1 => log10Rhos_ary
1081 0 : etas_f1 => etas_ary
1082 : nullify(dxdt_source_term, times)
1083 0 : log10Ts_f1(1) = s% lnT(k)/ln10
1084 0 : log10Rhos_f1(1) = s% lnd(k)/ln10
1085 0 : etas_f1(1) = s% eta(k)
1086 : call net_1_zone_burn( &
1087 : s% net_handle, s% eos_handle, species, s% num_reactions, 0d0, dt, xa_start, &
1088 : num_times, times, log10Ts_f1, log10Rhos_f1, etas_f1, dxdt_source_term, &
1089 : s% rate_factors, s% weak_rate_factor, &
1090 : std_reaction_Qs, std_reaction_neuQs, &
1091 : screening_mode, &
1092 : stptry, max_steps, eps, odescal, &
1093 : use_pivoting, trace, burn_dbg, burn_finish_substep, &
1094 : s% xa(1:species,k), &
1095 : s% eps_nuc_categories(:,k), &
1096 : avg_epsnuc, ending_eps_neu_total, &
1097 0 : nfcn, njac, ntry, naccpt, nrejct, ierr)
1098 0 : if (ierr /= 0) then
1099 0 : if (s% report_ierr) write(*,2) 'net_1_zone_burn failed', k
1100 0 : return
1101 : call mesa_error(__FILE__,__LINE__,'burn1_zone')
1102 : end if
1103 : end if
1104 :
1105 0 : s% raw_rate(:,k) = 0d0
1106 0 : s% screened_rate(:,k) = 0d0
1107 0 : s% eps_nuc_rate(:,k) = 0d0
1108 0 : s% eps_neu_rate(:,k) = 0d0
1109 :
1110 0 : num_iters_out = naccpt
1111 :
1112 : ! make extra call to get eps_nuc_categories
1113 0 : call do1_net(s, k, s% species, s% num_reactions, .false., ierr)
1114 0 : if (ierr /= 0) then
1115 0 : if (s% report_ierr) &
1116 0 : write(*,2) 'net_1_zone_burn final call to do1_net failed', k
1117 0 : return
1118 : call mesa_error(__FILE__,__LINE__,'burn1_zone')
1119 : end if
1120 :
1121 0 : s% eps_nuc(k) = 0d0
1122 0 : s% d_epsnuc_dlnd(k) = 0d0
1123 0 : s% d_epsnuc_dlnT(k) = 0d0
1124 0 : s% d_epsnuc_dx(:,k) = 0d0
1125 0 : s% dxdt_nuc(:,k) = 0d0
1126 : !s% eps_nuc_categories(:,k) = 0d0
1127 0 : s% d_dxdt_nuc_dRho(:,k) = 0d0
1128 0 : s% d_dxdt_nuc_dT(:,k) = 0d0
1129 0 : s% d_dxdt_nuc_dx(:,:,k) = 0d0
1130 : ! below, restore eps_nuc_neu to op_split zones.
1131 0 : s% eps_nuc_neu_total(k) = ending_eps_neu_total
1132 :
1133 0 : do i=1,species ! for use by dX_nuc_drop timestep limiter
1134 0 : s% dxdt_nuc(i,k) = (s% xa(i,k)-xa_start(i))/dt
1135 : end do
1136 :
1137 : contains
1138 :
1139 0 : subroutine get_eos_info_for_burn_at_const_density( &
1140 0 : eos_handle, species, chem_id, net_iso, xa, &
1141 : Rho, logRho, T, logT, &
1142 : Cv, d_Cv_dlnT, eta, d_eta_dlnT, ierr)
1143 : use eos_lib, only: eosDT_get
1144 : use eos_def
1145 : integer, intent(in) :: eos_handle, species
1146 : integer, pointer :: chem_id(:) ! maps species to chem id
1147 : integer, pointer :: net_iso(:) ! maps chem id to species number
1148 : real(dp), intent(in) :: &
1149 : xa(:), rho, logRho, T, logT
1150 : real(dp), intent(out) :: &
1151 : Cv, d_Cv_dlnT, eta, d_eta_dlnT
1152 : integer, intent(out) :: ierr
1153 :
1154 : real(dp), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
1155 0 : real(dp) :: d_dxa(num_eos_d_dxa_results,species)
1156 :
1157 : include 'formats'
1158 : ierr = 0
1159 :
1160 : call eosDT_get( &
1161 : eos_handle, species, chem_id, net_iso, xa, &
1162 : Rho, logRho, T, logT, &
1163 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
1164 :
1165 0 : if (ierr /= 0) then
1166 0 : write(*,*) 'failed in eosDT_get'
1167 : return
1168 : end if
1169 :
1170 0 : Cv = res(i_cv)
1171 0 : d_Cv_dlnT = d_dlnT(i_cv)
1172 :
1173 0 : eta = res(i_eta)
1174 0 : d_eta_dlnT = d_dlnT(i_eta)
1175 :
1176 : end subroutine get_eos_info_for_burn_at_const_density
1177 :
1178 :
1179 0 : subroutine burn_finish_substep(nstp, time, y, ierr)
1180 : integer,intent(in) :: nstp
1181 : real(dp), intent(in) :: time, y(:)
1182 : integer, intent(out) :: ierr
1183 : !real(dp) :: frac, step_time
1184 : !integer :: j, i
1185 : include 'formats'
1186 0 : ierr = 0
1187 : ! This routine does nothing other than set ierr = 0,
1188 : ! but we need an empty routine here because
1189 : ! net_1_zone_burn_const_density
1190 : ! expects to be passed a routine burn_finish_substep,
1191 : ! and often that will be a routine that actually does something,
1192 : ! but here we don't want to do anything.
1193 :
1194 : !step_time = time - substep_start_time
1195 : !if (step_time <= 0d0) return
1196 : !frac = step_time/dt
1197 : !do j = 1, num_categories
1198 : ! s% eps_nuc_categories(j,k) = &
1199 : ! s% eps_nuc_categories(j,k) + frac*eps_nuc_cat(j)
1200 : !end do
1201 : !if (.false. .and. k == s% nz) then
1202 : ! i = maxloc(eps_nuc_cat(1:num_categories),dim=1)
1203 : ! write(*,3) 'frac time/dt eps_nuc_cat ' // trim(category_name(i)), &
1204 : ! i, k, frac, time/dt, eps_nuc_cat(i), s% eps_nuc_categories(i,k)
1205 : !end if
1206 : !substep_start_time = time
1207 0 : end subroutine burn_finish_substep
1208 :
1209 : end subroutine burn1_zone
1210 :
1211 :
1212 : end module struct_burn_mix
1213 :
1214 :
|
