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