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