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