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