Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 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 evolve
21 :
22 : use star_private_def
23 : use const_def, only: dp, i8, pi, pi4, msun, lsun, crad, clight, secyer, secday
24 : use utils_lib, only: is_bad
25 : use star_utils
26 : use auto_diff_support
27 :
28 : implicit none
29 :
30 : private
31 : public :: do_evolve_step_part1
32 : public :: do_evolve_step_part2
33 : public :: pick_next_timestep
34 : public :: prepare_to_redo
35 : public :: prepare_to_retry
36 : public :: finish_step
37 : public :: set_age
38 :
39 : contains
40 :
41 11 : integer function do_evolve_step_part1(id, first_try)
42 : use alloc, only: fill_star_info_arrays_with_NaNs, &
43 : do_fill_arrays_with_NaNs, star_info_old_arrays
44 : logical, intent(in) :: first_try
45 : integer, intent(in) :: id
46 : type (star_info), pointer :: s
47 : integer :: ierr
48 : include 'formats'
49 :
50 11 : do_evolve_step_part1 = terminate
51 : ierr = 0
52 11 : call star_ptr(id, s, ierr)
53 11 : if (ierr /= 0) return
54 :
55 11 : if (s% trace_evolve) write(*,'(/,a)') 'start evolve step'
56 :
57 11 : if (is_bad(s% dt)) then
58 0 : write(*,1) 's% dt', s% dt
59 0 : call mesa_error(__FILE__,__LINE__,'do_evolve_step_part1')
60 : end if
61 :
62 11 : if (first_try .and. s% fill_arrays_with_NaNs .and. .not. s% RSP_flag) then
63 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
64 0 : write(*,*) 'fill_arrays_with_NaNs at start of step'
65 0 : call test_set_undefined
66 0 : call fill_star_info_arrays_with_NaNs(s, ierr)
67 0 : if (ierr /= 0) return
68 0 : call star_info_old_arrays(s, do_fill_arrays_with_NaNs, ierr)
69 0 : if (ierr /= 0) return
70 : end if
71 11 : do_evolve_step_part1 = do_step_part1(id, first_try)
72 11 : s% total_step_attempts = s% total_step_attempts + 1
73 11 : if (s% doing_relax) &
74 0 : s% total_relax_step_attempts = s% total_relax_step_attempts + 1
75 22 : if (do_evolve_step_part1 == redo) then
76 0 : s% total_step_redos = s% total_step_redos + 1
77 0 : if (s% doing_relax) &
78 0 : s% total_relax_step_redos = s% total_relax_step_redos + 1
79 11 : else if (do_evolve_step_part1 == retry) then
80 0 : s% total_step_retries = s% total_step_retries + 1
81 0 : if (s% doing_relax) &
82 0 : s% total_relax_step_retries = s% total_relax_step_retries + 1
83 : end if
84 :
85 : contains
86 :
87 0 : subroutine test_set_undefined
88 : ! should include everything in star_data_step_work.inc
89 : ! may be missing some. if so, please add them.
90 11 : use utils_lib, only: set_to_NaN
91 : integer :: j
92 : include 'formats'
93 :
94 0 : s% nz_old = -999
95 0 : s% model_number_old = -999
96 0 : s% prev_mesh_nz = -999
97 0 : s% num_conv_boundaries = -999
98 0 : s% num_mix_boundaries = -999
99 0 : s% num_mix_regions = -999
100 0 : s% num_mixing_regions = -999
101 0 : s% n_conv_regions = -999
102 0 : s% atm_structure_num_pts = -999
103 0 : s% generations = -999
104 0 : s% num_solver_iterations = -999
105 0 : s% num_diffusion_solver_iters = -999
106 0 : s% num_diffusion_solver_steps = -999
107 0 : s% num_rotation_solver_steps = -999
108 0 : s% result_reason = -999
109 0 : s% burner_num_threads = -999
110 0 : s% k_below_const_q = -999
111 0 : s% k_const_mass = -999
112 0 : s% k_for_test_CpT_absMdot_div_L = -999
113 0 : s% k_below_just_added = -999
114 0 : s% termination_code = -999
115 0 : s% dX_nuc_drop_max_k = -999
116 0 : s% dX_nuc_drop_max_j = -999
117 0 : s% solver_test_partials_var = -999
118 0 : s% solver_test_partials_dx_sink = -999
119 0 : s% n_conv_regions_old = -999
120 0 : do j=1,max_num_mixing_regions
121 0 : call set_to_NaN(s% cz_top_mass_old(j))
122 0 : call set_to_NaN(s% cz_bot_mass_old(j))
123 : end do
124 :
125 0 : do j=1,num_categories
126 0 : call set_to_NaN(s% L_by_category(j))
127 : end do
128 :
129 : ! sorted to help remove duplicates
130 0 : call set_to_NaN(s% L_center_old)
131 0 : call set_to_NaN(s% L_for_BB_outer_BC)
132 0 : call set_to_NaN(s% L_nuc_burn_total)
133 0 : call set_to_NaN(s% L_phot_old)
134 0 : call set_to_NaN(s% L_surf)
135 0 : call set_to_NaN(s% L_surf_old)
136 0 : call set_to_NaN(s% M_center_old)
137 0 : call set_to_NaN(s% R_center_old)
138 0 : call set_to_NaN(s% Teff_old)
139 0 : call set_to_NaN(s% acoustic_cutoff)
140 0 : call set_to_NaN(s% adjust_J_q)
141 0 : call set_to_NaN(s% adjust_mass_inner_frac_sub1)
142 0 : call set_to_NaN(s% adjust_mass_mid_frac_sub1)
143 0 : call set_to_NaN(s% adjust_mass_outer_frac_sub1)
144 0 : call set_to_NaN(s% angular_momentum_added)
145 0 : call set_to_NaN(s% angular_momentum_removed)
146 0 : call set_to_NaN(s% conv_mx1_bot)
147 0 : call set_to_NaN(s% conv_mx1_bot_r)
148 0 : call set_to_NaN(s% conv_mx1_top)
149 0 : call set_to_NaN(s% conv_mx1_top_r)
150 0 : call set_to_NaN(s% conv_mx2_bot)
151 0 : call set_to_NaN(s% conv_mx2_bot_r)
152 0 : call set_to_NaN(s% conv_mx2_top)
153 0 : call set_to_NaN(s% conv_mx2_top_r)
154 0 : call set_to_NaN(s% cumulative_energy_error_old)
155 0 : call set_to_NaN(s% cumulative_extra_heating_old)
156 0 : call set_to_NaN(s% d_vc_dv)
157 0 : call set_to_NaN(s% delta_Pg)
158 0 : call set_to_NaN(s% delta_nu)
159 0 : call set_to_NaN(s% dt_days)
160 0 : call set_to_NaN(s% dt_limit_ratio_old)
161 0 : call set_to_NaN(s% dt_old)
162 0 : call set_to_NaN(s% dt_start)
163 0 : call set_to_NaN(s% dt_years)
164 0 : call set_to_NaN(s% energy_change_from_do_adjust_mass_and_calculate_eps_mdot)
165 0 : call set_to_NaN(s% error_in_energy_conservation)
166 0 : call set_to_NaN(s% explicit_mstar_dot)
167 0 : call set_to_NaN(s% gradT_excess_max_lambda)
168 0 : call set_to_NaN(s% gradT_excess_min_beta)
169 0 : call set_to_NaN(s% h1_czb_mass)
170 0 : call set_to_NaN(s% initial_L_center)
171 0 : call set_to_NaN(s% initial_R_center)
172 0 : call set_to_NaN(s% initial_timestep)
173 0 : call set_to_NaN(s% initial_v_center)
174 0 : call set_to_NaN(s% max_conv_time_scale)
175 0 : call set_to_NaN(s% max_residual)
176 0 : call set_to_NaN(s% mdot_acoustic_surface)
177 0 : call set_to_NaN(s% mdot_adiabatic_surface)
178 0 : call set_to_NaN(s% mesh_adjust_IE_conservation)
179 0 : call set_to_NaN(s% mesh_adjust_KE_conservation)
180 0 : call set_to_NaN(s% mesh_adjust_PE_conservation)
181 0 : call set_to_NaN(s% min_conv_time_scale)
182 0 : call set_to_NaN(s% mstar_dot_old)
183 0 : call set_to_NaN(s% mstar_old)
184 0 : call set_to_NaN(s% mx1_bot)
185 0 : call set_to_NaN(s% mx1_bot_r)
186 0 : call set_to_NaN(s% mx1_top)
187 0 : call set_to_NaN(s% mx1_top_r)
188 0 : call set_to_NaN(s% mx2_bot)
189 0 : call set_to_NaN(s% mx2_bot_r)
190 0 : call set_to_NaN(s% mx2_top)
191 0 : call set_to_NaN(s% mx2_top_r)
192 0 : call set_to_NaN(s% non_epsnuc_energy_change_from_split_burn)
193 0 : call set_to_NaN(s% nu_max)
194 0 : call set_to_NaN(s% power_neutrinos)
195 0 : call set_to_NaN(s% power_nonnuc_neutrinos)
196 0 : call set_to_NaN(s% power_nuc_burn_old)
197 0 : call set_to_NaN(s% power_nuc_neutrinos)
198 0 : call set_to_NaN(s% prev_Ledd)
199 0 : call set_to_NaN(s% residual_norm)
200 0 : call set_to_NaN(s% rotational_mdot_boost)
201 0 : call set_to_NaN(s% shock_mass_start)
202 0 : call set_to_NaN(s% solver_test_partials_dval_dx)
203 0 : call set_to_NaN(s% solver_test_partials_val)
204 0 : call set_to_NaN(s% star_age)
205 0 : call set_to_NaN(s% starting_T_center)
206 0 : call set_to_NaN(s% super_eddington_wind_mdot)
207 0 : call set_to_NaN(s% surf_csound)
208 0 : call set_to_NaN(s% surf_r_equatorial)
209 0 : call set_to_NaN(s% surf_rho)
210 0 : call set_to_NaN(s% surface_cell_specific_total_energy_old)
211 0 : call set_to_NaN(s% tau_center)
212 0 : call set_to_NaN(s% time_days)
213 0 : call set_to_NaN(s% time_old)
214 0 : call set_to_NaN(s% time_step)
215 0 : call set_to_NaN(s% time_years)
216 0 : call set_to_NaN(s% total_WD_sedimentation_heating)
217 0 : call set_to_NaN(s% total_abs_angular_momentum)
218 0 : call set_to_NaN(s% total_angular_momentum)
219 0 : call set_to_NaN(s% total_angular_momentum_old)
220 0 : call set_to_NaN(s% total_energy)
221 0 : call set_to_NaN(s% total_energy_after_adjust_mass)
222 0 : call set_to_NaN(s% total_energy_before_adjust_mass)
223 0 : call set_to_NaN(s% total_energy_change_from_mdot)
224 0 : call set_to_NaN(s% total_energy_end)
225 0 : call set_to_NaN(s% total_energy_from_diffusion)
226 0 : call set_to_NaN(s% total_energy_from_phase_separation)
227 0 : call set_to_NaN(s% total_energy_old)
228 0 : call set_to_NaN(s% total_energy_sources_and_sinks)
229 0 : call set_to_NaN(s% total_energy_start)
230 0 : call set_to_NaN(s% total_eps_grav)
231 0 : call set_to_NaN(s% total_eps_mdot)
232 0 : call set_to_NaN(s% total_extra_heating)
233 0 : call set_to_NaN(s% total_gravitational_energy)
234 0 : call set_to_NaN(s% total_gravitational_energy_after_adjust_mass)
235 0 : call set_to_NaN(s% total_gravitational_energy_before_adjust_mass)
236 0 : call set_to_NaN(s% total_gravitational_energy_end)
237 0 : call set_to_NaN(s% total_gravitational_energy_initial)
238 0 : call set_to_NaN(s% total_gravitational_energy_old)
239 0 : call set_to_NaN(s% total_gravitational_energy_start)
240 0 : call set_to_NaN(s% total_internal_energy)
241 0 : call set_to_NaN(s% total_internal_energy_after_adjust_mass)
242 0 : call set_to_NaN(s% total_internal_energy_before_adjust_mass)
243 0 : call set_to_NaN(s% total_internal_energy_end)
244 0 : call set_to_NaN(s% total_internal_energy_initial)
245 0 : call set_to_NaN(s% total_internal_energy_old)
246 0 : call set_to_NaN(s% total_internal_energy_start)
247 0 : call set_to_NaN(s% total_irradiation_heating)
248 0 : call set_to_NaN(s% total_non_nuc_neu_cooling)
249 0 : call set_to_NaN(s% total_nuclear_heating)
250 0 : call set_to_NaN(s% total_radial_kinetic_energy)
251 0 : call set_to_NaN(s% total_radial_kinetic_energy_after_adjust_mass)
252 0 : call set_to_NaN(s% total_radial_kinetic_energy_before_adjust_mass)
253 0 : call set_to_NaN(s% total_radial_kinetic_energy_end)
254 0 : call set_to_NaN(s% total_radial_kinetic_energy_initial)
255 0 : call set_to_NaN(s% total_radial_kinetic_energy_old)
256 0 : call set_to_NaN(s% total_radial_kinetic_energy_start)
257 0 : call set_to_NaN(s% total_rotational_kinetic_energy)
258 0 : call set_to_NaN(s% total_rotational_kinetic_energy_after_adjust_mass)
259 0 : call set_to_NaN(s% total_rotational_kinetic_energy_before_adjust_mass)
260 0 : call set_to_NaN(s% total_rotational_kinetic_energy_end)
261 0 : call set_to_NaN(s% total_rotational_kinetic_energy_initial)
262 0 : call set_to_NaN(s% total_rotational_kinetic_energy_old)
263 0 : call set_to_NaN(s% total_rotational_kinetic_energy_start)
264 0 : call set_to_NaN(s% total_turbulent_energy)
265 0 : call set_to_NaN(s% total_turbulent_energy_after_adjust_mass)
266 0 : call set_to_NaN(s% total_turbulent_energy_before_adjust_mass)
267 0 : call set_to_NaN(s% total_turbulent_energy_end)
268 0 : call set_to_NaN(s% total_turbulent_energy_initial)
269 0 : call set_to_NaN(s% total_turbulent_energy_old)
270 0 : call set_to_NaN(s% total_turbulent_energy_start)
271 0 : call set_to_NaN(s% v_center_old)
272 0 : call set_to_NaN(s% virial_thm_P_avg)
273 0 : call set_to_NaN(s% work_inward_at_center)
274 0 : call set_to_NaN(s% work_outward_at_surface)
275 0 : call set_to_NaN(s% xmstar_old)
276 :
277 0 : end subroutine test_set_undefined
278 :
279 : end function do_evolve_step_part1
280 :
281 :
282 11 : integer function do_step_part1(id, first_try)
283 : use evolve_support, only: set_current_to_old
284 : use hydro_vars, only: set_vars
285 : use winds, only: set_mdot
286 : use alloc, only: check_sizes, fill_star_info_arrays_with_NaNs
287 : use do_one_utils, only: write_terminal_header
288 : use hydro_vars, only: set_vars_if_needed, set_vars, set_cgrav
289 : use mix_info, only: set_cz_bdy_mass
290 : use hydro_rotation, only: use_xh_to_update_i_rot, set_rotation_info
291 : use star_utils
292 : use report, only: do_report
293 : use rsp, only: rsp_total_energy_integrals
294 : logical, intent(in) :: first_try
295 : integer, intent(in) :: id
296 :
297 : type (star_info), pointer :: s
298 : integer :: ierr, k, nz
299 : integer(i8) :: time0, clock_rate
300 11 : real(dp) :: total_radiation
301 :
302 : logical, parameter :: dbg = .false.
303 :
304 : include 'formats'
305 :
306 11 : do_step_part1 = terminate
307 11 : ierr = 0
308 : call get_star_ptr(id, s, ierr)
309 11 : if (ierr /= 0) return
310 :
311 11 : s% termination_code = 0
312 11 : s% retry_message = ''
313 11 : s% doing_solver_iterations = .false.
314 11 : s% num_rotation_solver_steps = 0
315 11 : s% have_mixing_info = .false.
316 11 : s% rotational_mdot_boost = 0d0
317 11 : s% phase_sep_mixing_mass = -1d0 ! must be negative at start of step
318 11 : s% L_for_BB_outer_BC = -1 ! mark as not set
319 11 : s% need_to_setvars = .true. ! always start fresh
320 11 : s% okay_to_set_mixing_info = .true. ! set false by element diffusion
321 11 : s% okay_to_set_mlt_vc = .false. ! don't change mlt_vc until have set mlt_vc_old
322 :
323 11 : if (s% timestep_hold > s% model_number + 10000) then
324 0 : write(*,3) 'ERROR: s% timestep_hold', s% timestep_hold, s% model_number
325 0 : call mesa_error(__FILE__,__LINE__,'do_step_part1')
326 : end if
327 :
328 11 : if (s% u_flag .and. s% v_flag) then
329 0 : write(*,*) 'must not have both u_flag and v_flag at the same time'
330 0 : return
331 : end if
332 :
333 11 : call reset_starting_vectors(s)
334 11 : nz = s% nz
335 : call set_qs(s, nz, s% q, s% dq, ierr)
336 11 : if (failed('set_qs')) return
337 11 : call set_m_and_dm(s)
338 11 : call set_dm_bar(s, nz, s% dm, s% dm_bar)
339 :
340 11 : if (s% rotation_flag) then
341 : call set_cgrav(s, ierr)
342 0 : if (failed('set_cgrav')) return
343 0 : call use_xh_to_update_i_rot(s)
344 0 : s% total_angular_momentum = total_angular_momentum(s)
345 : ! set r and rmid from xh
346 0 : do k=1,nz
347 0 : call get_r_and_lnR_from_xh(s, k, s% r(k), s% lnR(k))
348 : end do
349 0 : call set_rmid(s, 1, nz, ierr)
350 0 : if (failed('set_rmid')) return
351 0 : call set_rotation_info(s, .true., ierr)
352 0 : if (failed('set_rotation_info')) return
353 : end if
354 :
355 11 : if (s% doing_first_model_of_run) then
356 1 : if (s% do_history_file) then
357 1 : if (first_try) then
358 1 : call write_terminal_header(s)
359 : else
360 0 : write(*,1) '1st model retry log10(dt/yr)', log10(s% dt/secyer)
361 : end if
362 : end if
363 1 : call system_clock(time0,clock_rate)
364 1 : s% starting_system_clock_time = time0
365 1 : s% system_clock_rate = clock_rate
366 1 : s% initial_timestep = s% dt_next
367 1 : s% initial_L_center = s% L_center
368 1 : s% initial_R_center = s% R_center
369 1 : s% initial_v_center = s% v_center
370 1 : s% timestep_hold = -111
371 1 : if (first_try) s% model_number_old = s% model_number
372 : end if
373 :
374 11 : call system_clock(s% system_clock_at_start_of_step, clock_rate)
375 :
376 11 : if (first_try) then ! i.e., not a redo or retry
377 11 : s% have_new_generation = .false.
378 11 : do_step_part1 = prepare_for_new_step(s)
379 11 : if (do_step_part1 /= keep_going) then
380 0 : if (s% report_ierr) &
381 0 : write(*,*) 'do_step_part1: prepare_for_new_step'
382 0 : return
383 : end if
384 11 : s% have_new_generation = .true.
385 11 : s% have_new_cz_bdy_info = .false.
386 11 : if ((s% steps_before_use_velocity_time_centering == 0) .or. &
387 : (s% steps_before_use_velocity_time_centering > 0 .and. &
388 : s% model_number >= s% steps_before_use_velocity_time_centering)) &
389 0 : s% using_velocity_time_centering = .true.
390 11 : if (.not. s% doing_relax .and. s% steps_before_use_TDC > 0) then
391 0 : if (s% model_number >= s% steps_before_use_TDC) then
392 0 : s% MLT_option = 'TDC'
393 : else
394 0 : s% MLT_option = 'Cox'
395 : end if
396 : end if
397 : end if
398 :
399 11 : call reset_epsnuc_vectors(s)
400 :
401 11 : call set_current_to_old(s)
402 11 : do_step_part1 = prepare_for_new_try(s)
403 11 : if (do_step_part1 /= keep_going) then
404 0 : if (s% report_ierr) &
405 0 : write(*,*) 'do_step_part1: prepare_for_new_try'
406 0 : return
407 : end if
408 :
409 : call set_start_of_step_info(s, 'after prepare_for_new_try', ierr) ! does set_vars_if_needed
410 11 : if (failed('set_start_of_step_info')) return
411 :
412 11 : if (.not. s% have_new_cz_bdy_info) then
413 : call set_cz_bdy_mass(s, ierr)
414 0 : if (failed('set_cz_bdy_mass')) return
415 : end if
416 :
417 11 : if (s% RSP_flag) then
418 0 : s% mstar_dot = 0
419 : call rsp_total_energy_integrals(s, &
420 : s% total_internal_energy_old, &
421 : s% total_gravitational_energy_old, &
422 : s% total_radial_kinetic_energy_old, &
423 : s% total_rotational_kinetic_energy_old, &
424 : s% total_turbulent_energy_old, &
425 0 : s% total_energy_old, total_radiation)
426 : else
427 11 : call set_mdot(s, s% L_phot*Lsun, s% mstar, s% Teff, ierr)
428 11 : if (failed('set_mdot')) return
429 : ! set energy info for new mesh
430 : call eval_total_energy_integrals(s, &
431 : s% total_internal_energy_old, &
432 : s% total_gravitational_energy_old, &
433 : s% total_radial_kinetic_energy_old, &
434 : s% total_rotational_kinetic_energy_old, &
435 : s% total_turbulent_energy_old, &
436 11 : s% total_energy_old)
437 : end if
438 :
439 11 : s% surface_cell_specific_total_energy_old = cell_specific_total_energy(s,1)
440 :
441 45 : if (.not. s% have_initial_energy_integrals) then
442 : s% total_internal_energy_initial = &
443 1 : s% total_internal_energy_old
444 : s% total_gravitational_energy_initial = &
445 1 : s% total_gravitational_energy_old
446 : s% total_radial_kinetic_energy_initial = &
447 1 : s% total_radial_kinetic_energy_old
448 : s% total_rotational_kinetic_energy_initial = &
449 1 : s% total_rotational_kinetic_energy_old
450 1 : s% total_turbulent_energy_initial = s% total_turbulent_energy_old
451 1 : s% total_energy_initial = s% total_energy_old
452 1 : s% have_initial_energy_integrals = .true.
453 : end if
454 :
455 : contains
456 :
457 33 : logical function failed(str)
458 : character (len=*), intent(in) :: str
459 33 : if (ierr == 0) then
460 33 : failed = .false.
461 : return
462 : end if
463 0 : failed = .true.
464 0 : do_step_part1 = retry
465 0 : if (s% report_ierr) write(*, *) 'do_step_part1 ' // trim(str)
466 0 : s% result_reason = nonzero_ierr
467 11 : end function failed
468 :
469 : end function do_step_part1
470 :
471 :
472 11 : integer function do_evolve_step_part2(id, first_try)
473 : logical, intent(in) :: first_try
474 : integer, intent(in) :: id
475 : type (star_info), pointer :: s
476 : integer :: ierr
477 11 : do_evolve_step_part2 = terminate
478 : ierr = 0
479 11 : call star_ptr(id, s, ierr)
480 11 : if (ierr /= 0) return
481 11 : do_evolve_step_part2 = do_step_part2(id, first_try)
482 11 : if (do_evolve_step_part2 == redo) then
483 0 : s% total_step_redos = s% total_step_redos + 1
484 0 : if (s% doing_relax) &
485 0 : s% total_relax_step_redos = s% total_relax_step_redos + 1
486 11 : else if (do_evolve_step_part2 == retry) then
487 0 : s% total_step_retries = s% total_step_retries + 1
488 0 : if (s% doing_relax) &
489 0 : s% total_relax_step_retries = s% total_relax_step_retries + 1
490 : else ! keep_going or terminate both count as finished
491 11 : s% total_steps_finished = s% total_steps_finished + 1
492 11 : if (s% doing_relax) &
493 0 : s% total_relax_steps_finished = s% total_relax_steps_finished + 1
494 : end if
495 11 : if (s% trace_evolve) write(*,'(/,a)') 'done evolve step'
496 : end function do_evolve_step_part2
497 :
498 :
499 11 : integer function do_step_part2(id, first_try)
500 : use num_def
501 : use chem_def
502 : use report, only: do_report, set_power_info
503 : use adjust_mass, only: do_adjust_mass
504 : use element_diffusion, only: do_element_diffusion, finish_element_diffusion
505 : use phase_separation, only: do_phase_separation
506 : use conv_premix, only: do_conv_premix
507 : use evolve_support, only: set_current_to_old
508 : use eps_mdot, only: calculate_eps_mdot
509 : use struct_burn_mix, only: do_struct_burn_mix
510 : use hydro_vars, only: set_vars_if_needed, set_vars, set_final_vars, set_cgrav
511 : use star_utils, only: start_time, update_time, get_phot_info
512 : use solve_omega_mix, only: do_solve_omega_mix
513 : use mix_info, only: set_cz_bdy_mass, set_mixing_info
514 : use hydro_rotation, only: set_rotation_info, set_i_rot
515 : use winds, only: set_mdot
516 : use star_utils, only: &
517 : eval_integrated_total_energy_profile, eval_deltaM_total_energy_integrals, &
518 : set_phase_of_evolution, set_luminosity_by_category
519 : use profile
520 :
521 : logical, intent(in) :: first_try
522 : integer, intent(in) :: id
523 :
524 : type (star_info), pointer :: s
525 : integer :: ierr, &
526 : k, mdot_redo_cnt, max_mdot_redo_cnt, nz
527 : integer(i8) :: time0, clock_rate
528 : logical :: trace, skip_global_corr_coeff_limit, &
529 : have_too_large_wind_mdot, have_too_small_wind_mdot, &
530 : ignored_first_step, was_in_implicit_wind_limit
531 11 : real(dp) :: w_div_w_crit, w_div_w_crit_prev, mstar_dot, mstar_dot_prev, abs_mstar_delta, &
532 11 : explicit_mdot, max_wind_mdot, wind_mdot, r_phot, kh_timescale, dmskhf, dmsfac, &
533 11 : too_large_wind_mdot, too_small_wind_mdot, boost, mstar_dot_nxt, &
534 11 : surf_omega_div_omega_crit_limit, dt
535 :
536 : integer :: ph_k, mdot_action
537 11 : real(dp) :: implicit_mdot, ph_L, iwind_tolerance, iwind_lambda
538 11 : real(dp) :: dummy1, dummy2, dummy3, dummy4, dummy5, dummy6, dummy7, dummy8
539 : integer :: iwind_redo_cnt, iwind_max_redo_cnt
540 : integer, parameter :: exit_loop = 1, cycle_loop = 0
541 :
542 : logical, parameter :: dbg = .false.
543 :
544 : include 'formats'
545 :
546 11 : ierr = 0
547 : call get_star_ptr(id, s, ierr)
548 11 : if (ierr /= 0) return
549 :
550 11 : if (s% dt <= 0d0) then
551 0 : do_step_part2 = terminate
552 0 : s% termination_code = t_dt_is_zero
553 0 : s% result_reason = dt_is_zero
554 0 : return
555 : end if
556 :
557 11 : time0 = s% starting_system_clock_time
558 11 : clock_rate = s% system_clock_rate
559 11 : trace = s% trace_evolve
560 11 : nz = s% nz
561 :
562 11 : call setup_for_implicit_mdot_loop
563 :
564 44 : implicit_mdot_loop: do
565 :
566 11 : dt = s% dt
567 11 : s% time = s% time_old + dt
568 11 : s% star_age = s% time/secyer
569 11 : s% okay_to_set_mixing_info = .true.
570 :
571 11 : if (s% v_center /= 0d0) then ! adjust R_center
572 0 : s% R_center = s% R_center_old + dt*s% v_center
573 0 : if (s% R_center < 0d0) then
574 0 : write(*,2) 's% R_center', s% model_number, s% R_center
575 0 : do_step_part2 = retry
576 0 : return
577 : end if
578 : end if
579 :
580 : call set_vars_if_needed(s, dt, 'start of implicit_mdot_loop', ierr)
581 11 : if (failed('set_vars_if_needed start of implicit_mdot_loop')) return
582 :
583 22 : if (s% RSP_flag) then
584 :
585 : call set_cgrav(s, ierr)
586 0 : if (failed('set_cgrav')) return
587 :
588 : else
589 :
590 : call do_adjust_mass(s, s% species, ierr)
591 11 : if (failed('do_adjust_mass')) return
592 11 : s% star_mdot = s% mstar_dot/(Msun/secyer) ! dm/dt in msolar per year
593 : call set_vars_if_needed(s, dt, 'after do_adjust_mass', ierr)
594 11 : if (failed('set_vars_if_needed after do_adjust_mass')) return
595 :
596 11 : call calculate_eps_mdot(s, dt, ierr)
597 11 : if (failed('calculate_eps_mdot')) return
598 :
599 11 : if (s% mstar_dot /= 0d0) then
600 : s% energy_change_from_do_adjust_mass_and_calculate_eps_mdot = &
601 0 : s% total_energy_after_adjust_mass - s% total_energy_before_adjust_mass
602 : else
603 11 : s% energy_change_from_do_adjust_mass_and_calculate_eps_mdot = 0d0
604 : end if
605 :
606 : call set_vars_if_needed(s, dt, 'after calculate_eps_mdot', ierr)
607 11 : if (failed('set_vars_if_needed after calculate_eps_mdot')) return
608 :
609 11 : if (s% do_conv_premix) then
610 0 : do k=1,s% nz
611 0 : s% eps_pre_mix(k) = s% energy(k)
612 : end do
613 : call do_conv_premix(s, ierr)
614 0 : if (failed('do_conv_premix')) return
615 : call set_vars_if_needed(s, dt, 'after do_conv_premix', ierr)
616 0 : if (failed('set_vars_if_needed after do_conv_premix')) return
617 0 : do k=1,s% nz ! for use by energy equation
618 0 : s% eps_pre_mix(k) = (s% eps_pre_mix(k) - s% energy(k)) / dt
619 : end do
620 : end if
621 :
622 11 : if(s% do_phase_separation) then
623 : call do_phase_separation(s, dt, ierr)
624 0 : if (failed('do_phase_separation')) return
625 :
626 : call set_vars_if_needed(s, dt, 'after phase separation', ierr)
627 0 : if (failed('set_vars_if_needed after phase separation')) return
628 : else
629 11 : s% crystal_core_boundary_mass = -1d0
630 : end if
631 :
632 11 : s% okay_to_set_mixing_info = .false. ! no mixing changes in set_vars after this point
633 :
634 11 : if (s% do_element_diffusion) then
635 : call do_element_diffusion(s, dt, ierr)
636 0 : if (failed('do_element_diffusion')) return
637 : call set_vars_if_needed(s, dt, 'after element diffusion', ierr)
638 0 : if (failed('set_vars_if_needed after element diffusion')) return
639 0 : call finish_element_diffusion(s,dt) ! calculates eps_diffusion from energy changes
640 : if (.false.) then
641 : write(*,1) 'dt*dm*eps_diffusion/total_energy_old', &
642 : dt*dot_product(s% dm(1:s% nz), s% eps_diffusion(1:s% nz))/s% total_energy_old
643 : end if
644 : end if
645 :
646 : end if
647 :
648 11 : if (s% rotation_flag .and. s% premix_omega .and. .not. s% j_rot_flag) then
649 0 : do_step_part2 = do_solve_omega_mix(s, 0.5d0*dt)
650 0 : if (do_step_part2 /= keep_going) return
651 0 : call set_rotation_info(s, .false., ierr)
652 0 : if (failed('set_rotation_info')) return
653 : call set_vars_if_needed(s, dt, 'after do_solve_omega_mix', ierr)
654 0 : if (failed('after do_solve_omega_mix')) return
655 : end if
656 :
657 11 : if (s% use_other_pressure) then
658 0 : call s% other_pressure(s% id, ierr)
659 0 : if (failed('other_pressure returned ierr')) return
660 : end if
661 : call set_vars_if_needed(s, dt, 'after other_pressure', ierr)
662 11 : if (failed('set_vars_if_needed after other_pressure')) return
663 :
664 : call check_for_extra_heat(s, ierr)
665 11 : if (failed('check_for_extra_heat')) return
666 : call set_vars_if_needed(s, dt, 'after check_for_extra_heat', ierr)
667 11 : if (failed('set_vars_if_needed after check_for_extra_heat')) return
668 :
669 11 : if (.not. s% have_new_cz_bdy_info) then
670 : call set_cz_bdy_mass(s, ierr)
671 0 : if (failed('set_cz_bdy_mass')) return
672 : end if
673 :
674 : skip_global_corr_coeff_limit = (first_try .or. &
675 11 : s% model_number_for_last_retry /= s% model_number) ! last alternative is for redo's
676 :
677 11 : s% doing_struct_burn_mix = .true.
678 11 : do_step_part2 = do_struct_burn_mix(s, skip_global_corr_coeff_limit)
679 11 : s% doing_struct_burn_mix = .false.
680 11 : if (do_step_part2 /= keep_going) return
681 : ! when reach here, have taken the step successfully
682 : ! but might not satisfy the implicit mdot requirements.
683 :
684 11 : mdot_action = select_mdot_action(ierr)
685 11 : if (failed('select_mdot_action')) return
686 11 : if (do_step_part2 /= keep_going) return
687 11 : if (mdot_action == exit_loop) exit implicit_mdot_loop
688 11 : if (s% trace_evolve) write(*,*) 'cycle implicit_mdot_loop'
689 :
690 : end do implicit_mdot_loop
691 :
692 11 : s% solver_iter = 0 ! to indicate that no longer doing solver iterations
693 :
694 11 : if (.not. s% RSP_flag) then
695 : call set_final_vars(s, dt, ierr)
696 11 : if (failed('set_final_vars')) return
697 11 : if (s% okay_to_set_mlt_vc .and. .not. s% have_mlt_vc) then
698 1 : s% have_mlt_vc = .true.
699 : end if
700 11 : s% okay_to_set_mlt_vc = .false.
701 : end if
702 :
703 11 : if (.not. okay_energy_conservation()) return
704 :
705 11 : if (s% max_timestep_hi_T_limit > 0 .and. &
706 : s% max_years_for_timestep /= s% hi_T_max_years_for_timestep) then
707 0 : if (maxval(s% T(1:nz)) >= s% max_timestep_hi_T_limit) then
708 0 : write(*,1) 'switch to high T max timesteps'
709 0 : s% max_years_for_timestep = s% hi_T_max_years_for_timestep
710 0 : s% max_timestep = secyer*s% max_years_for_timestep
711 : end if
712 : end if
713 :
714 11 : call set_luminosity_by_category(s) ! final values for use in selecting timestep
715 11 : call set_power_info(s)
716 :
717 11 : s% total_angular_momentum = total_angular_momentum(s)
718 :
719 : ! do_report needs to be called now because checks for redo and retry rely on
720 : ! data stored by do_report.
721 : call do_report(s, ierr)
722 11 : if (failed('do_report')) return
723 11 : call set_phase_of_evolution(s)
724 :
725 11 : call system_clock(time0,clock_rate)
726 11 : s% current_system_clock_time = time0
727 : s% total_elapsed_time = &
728 33 : dble(time0 - s% starting_system_clock_time)/dble(clock_rate)
729 :
730 :
731 : contains
732 :
733 :
734 132 : logical function failed(str)
735 : character (len=*), intent(in) :: str
736 121 : if (ierr == 0) then
737 121 : failed = .false.
738 : return
739 : end if
740 0 : failed = .true.
741 0 : do_step_part2 = retry
742 0 : if (s% report_ierr) write(*, *) 'do_step_part2: ' // trim(str)
743 0 : s% result_reason = nonzero_ierr
744 11 : end function failed
745 :
746 :
747 11 : subroutine setup_for_implicit_mdot_loop
748 :
749 11 : ignored_first_step = .false.
750 :
751 11 : mstar_dot = 0
752 11 : w_div_w_crit = -1
753 11 : surf_omega_div_omega_crit_limit = s% surf_omega_div_omega_crit_limit
754 11 : mdot_redo_cnt = 0
755 11 : max_mdot_redo_cnt = s% max_mdot_redo_cnt
756 :
757 11 : max_wind_mdot = 10*Msun/secyer
758 11 : have_too_large_wind_mdot = .false.
759 11 : have_too_small_wind_mdot = .false.
760 11 : too_large_wind_mdot = 0
761 11 : too_small_wind_mdot = 0
762 :
763 11 : explicit_mdot = s% mstar_dot
764 :
765 11 : was_in_implicit_wind_limit = s% was_in_implicit_wind_limit
766 11 : if(abs(s% mstar_dot_old) > 0) then
767 : if (was_in_implicit_wind_limit .and. &
768 0 : s% generations == 2 .and. &
769 : abs((s% mstar_dot-s% mstar_dot_old)/s% mstar_dot_old)+1 > &
770 : s% mdot_revise_factor) then
771 0 : write(*,*) "Skipping first step in implicit mdot"
772 0 : s% mstar_dot = s% mstar_dot_old
773 0 : mdot_redo_cnt = 1
774 0 : ignored_first_step = .true.
775 : end if
776 : end if
777 :
778 11 : abs_mstar_delta = 0
779 :
780 11 : iwind_redo_cnt = 0
781 11 : iwind_max_redo_cnt = s% max_tries_for_implicit_wind
782 11 : iwind_tolerance = s% iwind_tolerance
783 11 : iwind_lambda = s% iwind_lambda
784 :
785 11 : end subroutine setup_for_implicit_mdot_loop
786 :
787 :
788 11 : integer function select_mdot_action(ierr)
789 : use hydro_rotation, only: set_surf_avg_rotation_info
790 : integer, intent(out) :: ierr
791 : include 'formats'
792 :
793 11 : select_mdot_action = exit_loop
794 11 : if (s% mstar_dot == 0 .or. max_mdot_redo_cnt <= 0) return
795 : ! the test of max_mdot_redo_cnt <= 0 belongs here. it was erroneously placed
796 : ! after possible select_mdot_action = cycle_loop, return. BP: Apr 25, 2021.
797 :
798 0 : if (iwind_redo_cnt < iwind_max_redo_cnt .and. iwind_lambda > 0d0) then
799 : ! check mdot calculated at end of step
800 : call get_phot_info(s, dummy1, dummy2, dummy3, ph_L, dummy4, dummy5, dummy6, dummy7, dummy8, ph_k)
801 0 : call set_mdot(s, ph_L, s% mstar, s% Teff, ierr)
802 0 : if (ierr /= 0) then
803 0 : do_step_part2 = retry
804 0 : s% result_reason = nonzero_ierr
805 0 : if (s% report_ierr) write(*, *) 'do_step_part2: set_mdot ierr'
806 0 : return
807 : end if
808 0 : implicit_mdot = s% mstar_dot
809 0 : if (abs(explicit_mdot - implicit_mdot) > &
810 : abs(implicit_mdot)*iwind_tolerance) then
811 0 : call set_current_to_old(s) ! preparing for redo
812 : s% mstar_dot = explicit_mdot + &
813 0 : iwind_lambda*(implicit_mdot - explicit_mdot)
814 0 : explicit_mdot = s% mstar_dot
815 0 : do_step_part2 = prepare_for_new_try(s)
816 0 : if (do_step_part2 /= keep_going) return
817 0 : iwind_redo_cnt = iwind_redo_cnt + 1
818 0 : select_mdot_action = cycle_loop
819 0 : return
820 : end if
821 0 : iwind_max_redo_cnt = iwind_redo_cnt ! done with implicit wind
822 : end if
823 :
824 0 : if (.not. s% rotation_flag) then
825 11 : select_mdot_action = exit_loop
826 : return
827 : end if
828 :
829 : ! check for omega > omega_crit
830 :
831 0 : mstar_dot_prev = mstar_dot
832 0 : mstar_dot = s% mstar_dot
833 0 : wind_mdot = -s% mstar_dot
834 :
835 0 : if (mdot_redo_cnt == 1 .or. ignored_first_step) then
836 : ! this is the 1st correction to mdot
837 0 : r_phot = sqrt(s% L(1)/(pi*crad*clight*pow4(s% Teff)))
838 0 : kh_timescale = eval_kh_timescale(s% cgrav(1), s% mstar, r_phot, s% L(1))
839 0 : dmskhf = s% rotational_mdot_kh_fac
840 0 : dmsfac = s% rotational_mdot_boost_fac
841 0 : max_wind_mdot = dmskhf*s% mstar/kh_timescale
842 0 : if (wind_mdot > 0) max_wind_mdot = min(max_wind_mdot, wind_mdot*dmsfac)
843 : end if
844 :
845 0 : w_div_w_crit_prev = w_div_w_crit
846 : ! check the new w_div_w_crit to make sure not too large
847 0 : call set_surf_avg_rotation_info(s)
848 0 : w_div_w_crit = s% w_div_w_crit_avg_surf
849 :
850 0 : if (wind_mdot >= max_wind_mdot) then
851 0 : if (mdot_redo_cnt == 0) then
852 0 : write(*,*) 'cannot fix omega >= omega_crit -- mass loss already at max'
853 : else
854 0 : write(*,2) 'retry: at max wind mass loss', s% model_number, &
855 0 : log10(max_wind_mdot/(Msun/secyer))
856 0 : do_step_part2 = retry
857 0 : s% result_reason = nonzero_ierr
858 0 : return
859 : end if
860 0 : write(*,'(A)')
861 0 : if (w_div_w_crit > surf_omega_div_omega_crit_limit) then
862 0 : write(*,1) 'retry: w_div_w_crit > surf_omega_div_omega_crit_limit', &
863 0 : w_div_w_crit, surf_omega_div_omega_crit_limit
864 0 : do_step_part2 = retry
865 0 : s% result_reason = nonzero_ierr
866 0 : return
867 : end if
868 11 : select_mdot_action = exit_loop
869 : return
870 : end if
871 :
872 : ! NOTE: we assume that if surface omega/omega_crit (w_div_w_crit) is too large,
873 : ! then mass loss needs to be made larger to fix the problem.
874 : ! if that assumption is wrong,
875 : ! i.e. if bigger mass loss makes w_div_w_crit worse,
876 : ! then in an unstable situation and will remove mass until regain stability.
877 :
878 : if (w_div_w_crit <= surf_omega_div_omega_crit_limit &
879 0 : .and. mdot_redo_cnt == 0) then
880 0 : s% was_in_implicit_wind_limit = .false.
881 0 : select_mdot_action = exit_loop
882 0 : return
883 : end if
884 :
885 : if (w_div_w_crit <= surf_omega_div_omega_crit_limit &
886 0 : .and. s% mstar_dot == explicit_mdot) then
887 : !exit implicit_mdot_loop
888 11 : select_mdot_action = exit_loop
889 : return
890 : ! implicit scheme reached the limit set by the explicit_mdot;
891 : ! no problem; no redo required.
892 : end if
893 :
894 0 : s% was_in_implicit_wind_limit = .true.
895 :
896 0 : if (dt/secyer < s% min_years_dt_for_redo_mdot) then
897 : if (.true.) write(*,1) &
898 0 : 'dt too small for fix to fix w > w_crit; min_years_dt_for_redo_mdot', &
899 0 : dt/secyer, s% min_years_dt_for_redo_mdot
900 0 : select_mdot_action = exit_loop
901 0 : return
902 : end if
903 :
904 : ! if get here, need to revise mdot to fix w_div_w_crit
905 :
906 0 : mdot_redo_cnt = mdot_redo_cnt + 1
907 :
908 0 : if (mdot_redo_cnt == 1) then ! this is the 1st correction to mdot
909 :
910 0 : call set_current_to_old(s)
911 0 : do_step_part2 = prepare_for_new_try(s)
912 0 : if (do_step_part2 /= keep_going) return
913 :
914 0 : have_too_small_wind_mdot = .true.
915 0 : too_small_wind_mdot = wind_mdot
916 0 : if (s% mstar_dot < 0) then
917 0 : s% mstar_dot = mstar_dot*s% mdot_revise_factor
918 : else
919 0 : s% mstar_dot = mstar_dot/s% mdot_revise_factor
920 : end if
921 :
922 0 : if (-s% mstar_dot > max_wind_mdot) s% mstar_dot = -max_wind_mdot
923 :
924 0 : write(*,3) 'w > w_crit: revise mdot and redo', &
925 0 : s% model_number, mdot_redo_cnt, w_div_w_crit, &
926 0 : log10(abs(s% mstar_dot)/(Msun/secyer))
927 :
928 : !abs_mstar_delta = max(abs(s% mstar_dot), 1d-6*Msun/secyer)
929 0 : abs_mstar_delta = abs(s% mstar_dot)
930 :
931 0 : s% need_to_setvars = .true.
932 0 : select_mdot_action = cycle_loop
933 0 : return
934 :
935 0 : else if (mdot_redo_cnt == 2 .and. ignored_first_step) then
936 0 : abs_mstar_delta = abs(s% mstar_dot_old)
937 : end if
938 :
939 : ! have already done at least one correction -- check if okay now
940 : if (w_div_w_crit <= surf_omega_div_omega_crit_limit .and. &
941 0 : have_too_small_wind_mdot .and. &
942 : abs((wind_mdot-too_small_wind_mdot)/wind_mdot) < &
943 : s% surf_omega_div_omega_crit_tol) then
944 0 : write(*,3) 'OKAY', s% model_number, mdot_redo_cnt, w_div_w_crit, &
945 0 : log10(abs(s% mstar_dot)/(Msun/secyer))
946 0 : write(*,'(A)')
947 0 : select_mdot_action = exit_loop ! in bounds so accept it
948 0 : return
949 : end if
950 :
951 0 : if (mdot_redo_cnt >= max_mdot_redo_cnt) then
952 0 : if (max_mdot_redo_cnt > 0) then
953 0 : write(*,3) 'failed to fix w > w_crit: too many tries', &
954 0 : s% model_number, mdot_redo_cnt, w_div_w_crit, &
955 0 : log10(abs(s% mstar_dot)/(Msun/secyer))
956 0 : do_step_part2 = retry
957 0 : s% result_reason = nonzero_ierr
958 0 : return
959 : end if
960 11 : select_mdot_action = exit_loop
961 : return
962 : end if
963 :
964 : if (w_div_w_crit > surf_omega_div_omega_crit_limit &
965 : .and. w_div_w_crit_prev >= surf_omega_div_omega_crit_limit &
966 0 : .and. -mstar_dot >= max_wind_mdot) then
967 0 : write(*,3) 'failed to fix w > w_crit', &
968 0 : s% model_number, mdot_redo_cnt, w_div_w_crit, &
969 0 : log10(abs(s% mstar_dot)/(Msun/secyer))
970 0 : write(*,'(A)')
971 0 : do_step_part2 = retry
972 0 : s% result_reason = nonzero_ierr
973 0 : return
974 : end if
975 :
976 0 : if (w_div_w_crit >= surf_omega_div_omega_crit_limit) then ! wind too small
977 : !write(*,*) "entering too small wind mdot"
978 0 : if (.not. have_too_small_wind_mdot) then
979 : !write(*,*) "setting too small wind mdot"
980 0 : too_small_wind_mdot = wind_mdot
981 0 : have_too_small_wind_mdot = .true.
982 0 : else if (wind_mdot > too_small_wind_mdot) then
983 : !write(*,*) "changing too small wind mdot"
984 0 : too_small_wind_mdot = wind_mdot
985 : end if
986 0 : else if (w_div_w_crit < surf_omega_div_omega_crit_limit) then ! wind too large
987 : !write(*,*) "entering too large wind mdot"
988 0 : if (.not. have_too_large_wind_mdot) then
989 : !write(*,*) "setting too large wind mdot"
990 0 : too_large_wind_mdot = wind_mdot
991 0 : have_too_large_wind_mdot = .true.
992 0 : else if (wind_mdot < too_large_wind_mdot) then
993 : !write(*,*) "changing too large wind mdot"
994 0 : too_large_wind_mdot = wind_mdot
995 : end if
996 : end if
997 :
998 0 : call set_current_to_old(s)
999 0 : do_step_part2 = prepare_for_new_try(s)
1000 0 : if (do_step_part2 /= keep_going) return
1001 :
1002 0 : if (have_too_large_wind_mdot .and. have_too_small_wind_mdot) then
1003 0 : if (abs((too_large_wind_mdot-too_small_wind_mdot)/too_large_wind_mdot) &
1004 : < s% surf_omega_div_omega_crit_tol) then
1005 0 : write(*,*) "too_large_wind_mdot good enough, using it"
1006 0 : s% mstar_dot = -too_large_wind_mdot
1007 : else
1008 : ! have bracketing mdots; bisect for next one.
1009 0 : s% mstar_dot = -0.5d0*(too_large_wind_mdot + too_small_wind_mdot)
1010 0 : write(*,3) 'fix w > w_crit: bisect mdots and redo', &
1011 0 : s% model_number, mdot_redo_cnt, w_div_w_crit, &
1012 0 : log10(abs(s% mstar_dot)/(Msun/secyer)), &
1013 0 : log10(abs(too_large_wind_mdot)/(Msun/secyer)), &
1014 0 : log10(abs(too_small_wind_mdot)/(Msun/secyer))
1015 : end if
1016 :
1017 : else ! still have wind too small so boost it again
1018 0 : if (have_too_small_wind_mdot) then
1019 0 : if (mod(mdot_redo_cnt,2) == 1) then
1020 0 : boost = s% implicit_mdot_boost
1021 : ! increase mass loss
1022 0 : mstar_dot_nxt = mstar_dot - boost*abs_mstar_delta
1023 : else
1024 0 : if (mstar_dot < 0) then ! increase mass loss
1025 0 : mstar_dot_nxt = mstar_dot*s% mdot_revise_factor
1026 : else ! decrease mass gain
1027 0 : mstar_dot_nxt = mstar_dot/s% mdot_revise_factor
1028 : end if
1029 : end if
1030 : else
1031 0 : if (mod(mdot_redo_cnt,2) == 1) then
1032 0 : boost = s% implicit_mdot_boost
1033 : ! decrease mass loss
1034 0 : mstar_dot_nxt = mstar_dot + boost*abs_mstar_delta
1035 : else
1036 0 : if (mstar_dot < 0) then ! decrease mass loss
1037 0 : mstar_dot_nxt = mstar_dot/s% mdot_revise_factor
1038 : else ! increase mass gain
1039 0 : mstar_dot_nxt = mstar_dot*s% mdot_revise_factor
1040 : end if
1041 : end if
1042 : end if
1043 0 : if (mstar_dot_prev /= explicit_mdot) &
1044 0 : mstar_dot_nxt = min(mstar_dot_nxt, explicit_mdot)
1045 0 : if (mstar_dot_nxt == explicit_mdot) &
1046 0 : write(*,*) "implicit mdot: reached explicit_mdot"
1047 0 : s% mstar_dot = mstar_dot_nxt
1048 0 : if (-s% mstar_dot > max_wind_mdot) s% mstar_dot = -max_wind_mdot
1049 : !abs_mstar_delta = max(abs_mstar_delta, abs(s% mstar_dot))
1050 0 : write(*,3) 'fix w > w_crit: change mdot and redo', &
1051 0 : s% model_number, mdot_redo_cnt, w_div_w_crit, &
1052 0 : log10(abs(s% mstar_dot)/(Msun/secyer))
1053 : end if
1054 :
1055 11 : select_mdot_action = cycle_loop ! cycle
1056 :
1057 11 : end function select_mdot_action
1058 :
1059 :
1060 : subroutine show_debug
1061 : integer :: k
1062 : real(dp) :: alfa, beta, gamma1, Cv, chiRho, chiT, Cp, grada, &
1063 : Pgas, Prad, P, opacity
1064 : include 'formats'
1065 : k = 1205
1066 :
1067 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
1068 : beta = 1 - alfa
1069 :
1070 : gamma1 = alfa*s% gamma1(k) + beta*s% gamma1(k-1)
1071 : Cv = alfa*s% Cv(k) + beta*s% Cv(k-1)
1072 : chiRho = alfa*s% chiRho(k) + beta*s% chiRho(k-1)
1073 : chiT = alfa*s% chiT(k) + beta*s% chiT(k-1)
1074 : Cp = alfa*s% Cp(k) + beta*s% Cp(k-1)
1075 : grada = alfa*s% grada(k) + beta*s% grada(k-1)
1076 : Pgas = alfa*s% Pgas(k) + beta*s% Pgas(k-1)
1077 : Prad = alfa*s% Prad(k) + beta*s% Prad(k-1)
1078 : P = alfa*s% Peos(k) + beta*s% Peos(k-1)
1079 : opacity = alfa*s% opacity(k) + beta*s% opacity(k-1)
1080 :
1081 : write(*,2) 'at end of step', s% model_number
1082 : write(*,2) 'gamma1', k, gamma1
1083 : write(*,2) 'Cv', k, Cv
1084 : write(*,2) 'chiRho', k, chiRho
1085 : write(*,2) 'chiT', k, chiT
1086 : write(*,2) 'Cp', k, Cp
1087 : write(*,2) 'grada', k, grada
1088 : write(*,2) 'Pgas', k, Pgas
1089 : write(*,2) 'Prad', k, Prad
1090 : write(*,2) 'P', k, P
1091 : write(*,2) 'opacity', k, opacity
1092 : write(*,2) 'L', k, s% L(k)
1093 : write(*,2) 'gradr', k, s% gradr(k)
1094 : write(*,2) 'gradr/grada', k, s% gradr(k)/grada
1095 : write(*,3) 'mixing_type', k, s% mixing_type(k)
1096 : write(*,'(A)')
1097 :
1098 11 : end subroutine show_debug
1099 :
1100 :
1101 11 : logical function okay_energy_conservation()
1102 : use rsp, only: rsp_total_energy_integrals
1103 : integer :: nz, k
1104 11 : real(dp) :: phase1_sources_and_sinks, phase2_sources_and_sinks, phase2_work, &
1105 11 : phase1_total_energy_from_mdot, phase2_total_energy_from_mdot, &
1106 11 : expected_sum_cell_others, expected_sum_cell_sources, L_theta, &
1107 11 : diff_total_gravitational_energy, diff_total_internal_energy, diff_total_kinetic_energy, &
1108 11 : diff_total_turbulent_energy, &
1109 11 : virial, total_radiation, L_surf, sum_cell_de, sum_cell_detrb, &
1110 11 : sum_cell_dke, sum_cell_dpe, sum_cell_dL, sum_cell_ergs_error, sum_cell_others, &
1111 11 : sum_cell_sources, sum_cell_terms, sum_cell_work, total_energy_from_pre_mixing,&
1112 11 : total_energy_from_fixed_m_grav
1113 :
1114 : include 'formats'
1115 :
1116 : ! phase1 := from end of previous step until start of solver
1117 : ! phase2 := from start of solver to end of step
1118 : !
1119 : ! total_energy_old = at beginning of phase1
1120 : ! total_energy_start = at beginning of phase2
1121 : ! total_energy_end = at end of phase2
1122 : !
1123 : ! phase1_energy_error = total_energy_start - (total_energy_old + phase1_sources_and_sinks)
1124 : ! phase2_energy_error = total_energy_end - (total_energy_start + phase2_sources_and_sinks)
1125 : !
1126 : ! total_energy_sources_and_sinks = phase1_sources_and_sinks + phase2_sources_and_sinks
1127 : !
1128 : ! error_in_energy_conservation = total_energy_end - (total_energy_old + total_energy_sources_and_sinks)
1129 : !
1130 : ! equivalently, error_in_energy_conservation = phase1_energy_error + phase2_energy_error
1131 :
1132 :
1133 11 : okay_energy_conservation = .false.
1134 :
1135 11 : nz = s% nz
1136 11 : if (s% RSP_flag) then
1137 : call rsp_total_energy_integrals(s, &
1138 : s% total_internal_energy_end, &
1139 : s% total_gravitational_energy_end, &
1140 : s% total_radial_kinetic_energy_end, &
1141 : s% total_rotational_kinetic_energy_end, &
1142 : s% total_turbulent_energy_end, &
1143 0 : s% total_energy_end, total_radiation)
1144 0 : if (s% RSP_just_set_velocities) then ! reset everything when 1st set velocities
1145 0 : s% total_internal_energy_old = s% total_internal_energy_end
1146 0 : s% total_gravitational_energy_old = s% total_gravitational_energy_end
1147 0 : s% total_radial_kinetic_energy_old = s% total_radial_kinetic_energy_end
1148 0 : s% total_rotational_kinetic_energy_old = s% total_rotational_kinetic_energy_end
1149 0 : s% total_turbulent_energy_old = s% total_turbulent_energy_end
1150 0 : s% total_energy_old = s% total_energy_end
1151 0 : total_radiation = 0d0
1152 : end if
1153 : else
1154 : call eval_total_energy_integrals(s, &
1155 : s% total_internal_energy_end, &
1156 : s% total_gravitational_energy_end, &
1157 : s% total_radial_kinetic_energy_end, &
1158 : s% total_rotational_kinetic_energy_end, &
1159 : s% total_turbulent_energy_end, &
1160 11 : s% total_energy_end)
1161 : end if
1162 :
1163 11 : if (s% mstar_dot == 0d0) then
1164 11 : s% total_energy_change_from_mdot = 0d0
1165 11 : s% total_eps_mdot = 0d0
1166 : else
1167 : s% total_energy_change_from_mdot = &
1168 0 : s% mstar_dot*dt*s% surface_cell_specific_total_energy_old
1169 0 : s% total_eps_mdot = dt*dot_product(s% dm(1:nz), s% eps_mdot(1:nz))
1170 : end if
1171 :
1172 13098 : virial = 3*sum(s% dm(1:nz)*s% Peos(1:nz)/s% rho(1:nz))
1173 11 : s% virial_thm_P_avg = virial
1174 :
1175 13098 : s% total_eps_grav = dt*dot_product(s% dm(1:nz), s% eps_grav_ad(1:nz)% val)
1176 :
1177 11 : if (s% u_flag .and. s% total_eps_grav /= 0d0) then
1178 0 : write(*,2) 'u_flag energy accounting ignores total_eps_grav', s% model_number, s% total_eps_grav
1179 0 : s% total_eps_grav = 0
1180 : end if
1181 :
1182 : ! notes from Adam:
1183 : ! When there are mass changes the total energy of the model changes.
1184 : ! We can split this change into three parts:
1185 : ! 1. Mass flows into or out of the model with some specific energy.
1186 : ! 2. There is work done at the surface of the model by pushing material past the pressure of the surface face.
1187 : ! 3. The mass in the model changes state. Near the surface matter changes to maintain the same rho(q) and T(q).
1188 : ! Below the surface regions there is a transition region where the state interpolates between fixed-m and fixed-q.
1189 : ! Still deeper the state is approximately maintained at fixed rho(m), T(m).
1190 : !
1191 : ! Change (1) is accounted for entirely by the term s% total_energy_change_from_mdot.
1192 : ! Change (2) is accounted for entirely by the term s% mdot_acoustic_surface.
1193 : !
1194 : ! Change (3) is accounted for by the term eps_mdot in the energy equation and the term
1195 : ! mdot_adiabatic_surface in the energy accounting.
1196 : !
1197 : ! The eps_mdot term accounts for the energy required to change the state of the matter
1198 : ! from its state before adjust_mass to its state after adjust_mass. Matter not present in the
1199 : ! model before adjust_mass is assumed to be in the same state as the surface cell, or else in the
1200 : ! state specified by the other_accreting_surface hook (if used). Matter not present in the model
1201 : ! after adjust_mass is in a state calculated by comparing the thermal and mass-loss time-scales,
1202 : ! and differences between this and the surface state are accounted for by the term mdot_adiabatic_surface.
1203 : !
1204 : ! By adding eps_mdot, we cause a change in energy during the Newton iterations which
1205 : ! cancels the change in (3). Thus eps_mdot does not enter into the energy *accounting*, just into
1206 : ! the energy equation. A consequence of this is that the sum
1207 : !
1208 : ! total_energy_change_from_mdot + (total energy before adjust_mass) + mdot_acoustic_surface
1209 : !
1210 : ! does not equal the total energy *after* adjust_mass and *before* the Newton iterations.
1211 : ! However it should equal the total energy at the end of the step.
1212 :
1213 : if (s% rotation_flag .and. &
1214 : (s% use_other_torque .or. s% use_other_torque_implicit .or. &
1215 : associated(s% binary_other_torque))) then
1216 : ! keep track of rotational kinetic energy
1217 : end if
1218 :
1219 11 : if (s% eps_nuc_factor == 0d0 .or. s% nonlocal_NiCo_decay_heat) then
1220 0 : s% total_nuclear_heating = 0d0
1221 11 : else if (s% op_split_burn) then
1222 0 : s% total_nuclear_heating = 0d0
1223 0 : do k = 1, nz
1224 0 : if (s% T_start(k) >= s% op_split_burn_min_T) then
1225 : s% total_nuclear_heating = s% total_nuclear_heating + &
1226 0 : dt*s% dm(k)*s% burn_avg_epsnuc(k)
1227 : else
1228 : s% total_nuclear_heating = s% total_nuclear_heating + &
1229 0 : dt*s% dm(k)*s% eps_nuc(k)
1230 : end if
1231 : end do
1232 : else
1233 13098 : s% total_nuclear_heating = dt*dot_product(s% dm(1:nz), s% eps_nuc(1:nz))
1234 : end if
1235 :
1236 11 : if (s% RSP_flag) then
1237 0 : s% total_non_nuc_neu_cooling = 0d0
1238 0 : s% total_irradiation_heating = 0d0
1239 : else
1240 : s% total_non_nuc_neu_cooling = dt*0.5d0* &
1241 13098 : sum((s% non_nuc_neu(1:nz) + s% non_nuc_neu_start(1:nz))*s% dm(1:nz))
1242 : s% total_irradiation_heating = &
1243 13098 : dt*dot_product(s% dm(1:nz), s% irradiation_heat(1:nz))
1244 : end if
1245 :
1246 11 : s% total_WD_sedimentation_heating = 0d0
1247 11 : if (s% do_element_diffusion .and. s% do_WD_sedimentation_heating) then
1248 : s% total_WD_sedimentation_heating = &
1249 0 : dt*dot_product(s% dm(1:nz), s% eps_WD_sedimentation(1:nz))
1250 : end if
1251 :
1252 11 : s% total_energy_from_diffusion = 0d0
1253 11 : if (s% do_element_diffusion .and. s% do_diffusion_heating) then
1254 : s% total_energy_from_diffusion = &
1255 0 : dt*dot_product(s% dm(1:nz), s% eps_diffusion(1:nz))
1256 : end if
1257 :
1258 11 : total_energy_from_pre_mixing = 0d0
1259 11 : if (s% do_conv_premix) then
1260 : total_energy_from_pre_mixing = &
1261 0 : dt*dot_product(s% dm(1:nz), s% eps_pre_mix(1:nz))
1262 : end if
1263 :
1264 11 : s% total_energy_from_phase_separation = 0d0
1265 11 : if (s% do_phase_separation) then
1266 : s% total_energy_from_phase_separation = &
1267 0 : dt*dot_product(s% dm(1:nz), s% eps_phase_separation(1:nz))
1268 : end if
1269 :
1270 : phase2_total_energy_from_mdot = &
1271 13098 : dt*dot_product(s% dm(1:nz), s% eps_mdot(1:nz))
1272 :
1273 13098 : s% total_extra_heating = dt*dot_product(s% dm(1:nz), s% extra_heat(1:nz)%val)
1274 :
1275 11 : phase2_work = dt*(s% work_outward_at_surface - s% work_inward_at_center)
1276 :
1277 11 : if (.not. s% RSP_flag) then
1278 11 : if (s% using_velocity_time_centering .and. &
1279 : s% include_L_in_velocity_time_centering) then
1280 0 : L_theta = s% L_theta_for_velocity_time_centering
1281 : else
1282 : L_theta = 1d0
1283 : end if
1284 11 : L_surf = L_theta*s% L(1) + (1d0 - L_theta)*s% L_start(1)
1285 11 : total_radiation = dt*(L_surf - s% L_center)
1286 : end if
1287 :
1288 : !note: phase1_total_energy_from_mdot = &
1289 : ! s% mdot_acoustic_surface &
1290 : ! + s% total_energy_change_from_mdot &
1291 : ! + s% mdot_adiabatic_surface &
1292 : ! - phase2_total_energy_from_mdot
1293 :
1294 : ! during the newton iterations, m_grav remains fixed,
1295 : ! even though the mass_corrections respond to the change in composition.
1296 : ! This means there is a term
1297 : ! -sum dm G Delta(m_grav) / r
1298 : ! being neglected.
1299 : ! after the iterations, we update m_grav before doing the energy accounting.
1300 : ! this changes the potential energy so we need to correct by that amount.
1301 11 : if (s% use_mass_corrections) then
1302 : total_energy_from_fixed_m_grav = &
1303 0 : -dot_product(s% cgrav(1:nz)/s%r(1:nz)*(s%m_grav(1:nz)-s%m_grav_start(1:nz)), s% dm(1:nz))
1304 : else
1305 : total_energy_from_fixed_m_grav = 0
1306 : end if
1307 :
1308 : phase1_total_energy_from_mdot = &
1309 : s% energy_change_from_do_adjust_mass_and_calculate_eps_mdot &
1310 11 : + s% mdot_adiabatic_surface ! ??
1311 :
1312 : phase1_sources_and_sinks = &
1313 : phase1_total_energy_from_mdot &
1314 : + total_energy_from_pre_mixing &
1315 : + s% total_energy_from_phase_separation &
1316 : + s% total_WD_sedimentation_heating &
1317 : + s% total_energy_from_diffusion &
1318 11 : + s% non_epsnuc_energy_change_from_split_burn
1319 :
1320 : phase2_sources_and_sinks = &
1321 : - total_energy_from_pre_mixing &
1322 : - s% total_energy_from_phase_separation &
1323 : - s% total_WD_sedimentation_heating &
1324 : - s% total_energy_from_diffusion &
1325 : + phase2_total_energy_from_mdot &
1326 : + s% total_nuclear_heating &
1327 : - s% total_non_nuc_neu_cooling &
1328 : + s% total_irradiation_heating &
1329 : + s% total_extra_heating &
1330 : - total_radiation &
1331 11 : - phase2_work
1332 :
1333 : s% total_energy_sources_and_sinks = &
1334 11 : phase1_sources_and_sinks + phase2_sources_and_sinks
1335 11 : if (is_bad(s% total_energy_sources_and_sinks)) then
1336 0 : write(*,2) 's% total_energy_sources_and_sinks', s% model_number, s% total_energy_sources_and_sinks
1337 0 : write(*,2) 'phase1_sources_and_sinks', s% model_number, phase1_sources_and_sinks
1338 0 : write(*,2) 'phase2_sources_and_sinks', s% model_number, phase2_sources_and_sinks
1339 0 : write(*,2) 'total_energy_from_pre_mixing', s% model_number, total_energy_from_pre_mixing
1340 0 : write(*,2) 's% total_WD_sedimentation_heating', s% model_number, s% total_WD_sedimentation_heating
1341 0 : write(*,2) 'phase2_total_energy_from_mdot', s% model_number, phase2_total_energy_from_mdot
1342 0 : write(*,2) 's% total_nuclear_heating', s% model_number, s% total_nuclear_heating
1343 0 : write(*,2) 's% total_non_nuc_neu_cooling', s% model_number, s% total_non_nuc_neu_cooling
1344 0 : write(*,2) 's% total_irradiation_heating', s% model_number, s% total_irradiation_heating
1345 0 : write(*,2) 's% total_extra_heating', s% model_number, s% total_extra_heating
1346 0 : write(*,2) 'phase2_work', s% model_number, phase2_work
1347 0 : write(*,2) 's% work_outward_at_surface', s% model_number, s% work_outward_at_surface
1348 0 : write(*,2) 's% work_inward_at_center', s% model_number, s% work_inward_at_center
1349 : !write(*,2) '', s% model_number,
1350 : !write(*,2) '', s% model_number,
1351 0 : call mesa_error(__FILE__,__LINE__,'okay_energy_conservation')
1352 : end if
1353 :
1354 : s% error_in_energy_conservation = &
1355 11 : s% total_energy_end - (s% total_energy_old + s% total_energy_sources_and_sinks)
1356 :
1357 11 : if (s% absolute_cumulative_energy_err) then
1358 : s% cumulative_energy_error = s% cumulative_energy_error_old + &
1359 11 : abs(s% error_in_energy_conservation)
1360 : else
1361 : s% cumulative_energy_error = s% cumulative_energy_error_old + &
1362 0 : s% error_in_energy_conservation
1363 : end if
1364 :
1365 11 : s% total_internal_energy = s% total_internal_energy_end
1366 11 : s% total_gravitational_energy = s% total_gravitational_energy_end
1367 11 : s% total_radial_kinetic_energy = s% total_radial_kinetic_energy_end
1368 11 : s% total_rotational_kinetic_energy = s% total_rotational_kinetic_energy_end
1369 11 : s% total_turbulent_energy = s% total_turbulent_energy_end
1370 11 : s% total_energy = s% total_energy_end
1371 :
1372 : ! provide info about non-conservation due to mass corrections
1373 11 : if (s% use_mass_corrections) then
1374 0 : write(*,2) 'INFO: use_mass_corrections incurred rel_E_err', &
1375 0 : s% model_number, -total_energy_from_fixed_m_grav/s% total_energy
1376 : end if
1377 :
1378 : if (s% model_number == s% energy_conservation_dump_model_number &
1379 11 : .and. .not. s% doing_relax) then
1380 :
1381 0 : write(*,'(A)')
1382 0 : write(*,2) 's% error_in_energy_conservation', s% model_number, s% error_in_energy_conservation
1383 0 : write(*,2) 'total_energy', s% model_number, s% total_energy
1384 0 : write(*,2) 'rel_E_err = error/total_energy', s% model_number, s% error_in_energy_conservation/s% total_energy
1385 0 : write(*,2) 'rel err phase1', s% model_number, &
1386 0 : (s% total_energy_start - (s% total_energy_old + phase1_sources_and_sinks))/s% total_energy
1387 0 : write(*,2) 'rel err phase2', s% model_number, &
1388 0 : (s% total_energy_end - (s% total_energy_start + phase2_sources_and_sinks))/s% total_energy
1389 0 : write(*,'(A)')
1390 0 : write(*,2) 's% total_energy_old', s% model_number, s% total_energy_old
1391 0 : write(*,2) 's% total_energy_start', s% model_number, s% total_energy_start
1392 0 : write(*,2) 's% total_energy_end', s% model_number, s% total_energy_end
1393 0 : write(*,2) 's% total_energy_sources_and_sinks', s% model_number, s% total_energy_sources_and_sinks
1394 0 : write(*,'(A)')
1395 :
1396 0 : if (trim(s% energy_eqn_option) == 'dedt') then
1397 :
1398 0 : write(*,'(A)')
1399 0 : write(*,*) 'for debugging phase1_sources_and_sinks'
1400 0 : write(*,'(A)')
1401 0 : write(*,2) 'total_energy_from_pre_mixing', s% model_number, total_energy_from_pre_mixing
1402 0 : write(*,2) 's% total_WD_sedimentation_heating', s% model_number, s% total_WD_sedimentation_heating
1403 0 : write(*,2) 's% total_energy_from_diffusion', s% model_number, s% total_energy_from_diffusion
1404 0 : write(*,2) 's% non_epsnuc_energy_change_from_split_burn', &
1405 0 : s% model_number, s% non_epsnuc_energy_change_from_split_burn
1406 0 : write(*,2) 'phase2 sum cell dt*dm*eps_mdot', s% model_number, phase2_total_energy_from_mdot
1407 0 : write(*,2) 'phase1_total_energy_from_mdot', s% model_number, phase1_total_energy_from_mdot
1408 0 : write(*,2) 'from_do_adjust_mass_and_eps_mdot', s% model_number, &
1409 0 : s% energy_change_from_do_adjust_mass_and_calculate_eps_mdot
1410 0 : write(*,2) 's% mdot_acoustic_surface', s% model_number, s% mdot_acoustic_surface
1411 0 : write(*,2) 's% mdot_adiabatic_surface', s% model_number, s% mdot_adiabatic_surface
1412 0 : write(*,2) 'phase2_total_energy_from_mdot', s% model_number, phase2_total_energy_from_mdot
1413 :
1414 0 : write(*,'(A)')
1415 0 : write(*,2) 's% mdot_acoustic_surface', s% model_number, s% mdot_acoustic_surface
1416 0 : write(*,2) 's% mdot_adiabatic_surface', s% model_number, s% mdot_adiabatic_surface
1417 0 : write(*,2) 's% total_energy_change_from_mdot', s% model_number, s% total_energy_change_from_mdot
1418 0 : write(*,2) 'phase1_sources_and_sinks', s% model_number, phase1_sources_and_sinks
1419 0 : write(*,*)
1420 0 : write(*,2) 'energy_start - energy_old', s% model_number, s% total_energy_start - s% total_energy_old
1421 0 : write(*,2) 'err phase1_sources_and_sinks', s% model_number, &
1422 0 : s% total_energy_start - (s% total_energy_old + phase1_sources_and_sinks)
1423 0 : write(*,2) 'rel err phase1_sources_and_sinks', s% model_number, &
1424 0 : (s% total_energy_start - (s% total_energy_old + phase1_sources_and_sinks))/s% total_energy
1425 0 : write(*,'(A)')
1426 0 : write(*,'(A)')
1427 :
1428 :
1429 0 : write(*,*) 'for debugging phase2_sources_and_sinks'
1430 0 : write(*,'(A)')
1431 :
1432 0 : write(*,2) 's% total_nuclear_heating', s% model_number, s% total_nuclear_heating
1433 0 : write(*,2) 's% total_non_nuc_neu_cooling', s% model_number, s% total_non_nuc_neu_cooling
1434 0 : write(*,2) 's% total_irradiation_heating', s% model_number, s% total_irradiation_heating
1435 0 : write(*,2) 's% total_extra_heating', s% model_number, s% total_extra_heating
1436 0 : write(*,'(A)')
1437 0 : write(*,2) 'total_energy_from_pre_mixing', s% model_number, total_energy_from_pre_mixing
1438 0 : write(*,2) 's% total_WD_sedimentation_heating', s% model_number, s% total_WD_sedimentation_heating
1439 0 : write(*,2) 's% total_energy_from_diffusion', s% model_number, s% total_energy_from_diffusion
1440 0 : write(*,'(A)')
1441 0 : write(*,2) 's% total_energy_change_from_mdot', s% model_number, s% total_energy_change_from_mdot
1442 0 : write(*,2) 's% mdot_acoustic_surface', s% model_number, s% mdot_acoustic_surface
1443 0 : write(*,2) 's% mdot_adiabatic_surface', s% model_number, s% mdot_adiabatic_surface
1444 : ! write(*,2) 'phase2_total_energy_from_mdot', s% model_number, phase2_total_energy_from_mdot
1445 0 : write(*,'(A)')
1446 0 : write(*,2) 'phase2_work', s% model_number, phase2_work
1447 0 : write(*,2) 'total_radiation', s% model_number, total_radiation
1448 0 : write(*,2) 's% non_epsnuc_energy_change_from_split_burn', s% model_number, &
1449 0 : s% non_epsnuc_energy_change_from_split_burn
1450 0 : write(*,'(A)')
1451 :
1452 0 : write(*,2) 's% work_outward_at_surface', s% model_number, s% work_outward_at_surface
1453 0 : write(*,2) 's% work_inward_at_center', s% model_number, s% work_inward_at_center
1454 0 : write(*,2) 'L_surf', s% model_number, L_surf
1455 0 : write(*,2) 'L_center', s% model_number, s% L_center
1456 0 : write(*,'(A)')
1457 :
1458 0 : sum_cell_dL = dt*dot_product(s% dm(1:nz), s% dL_dm(1:nz))
1459 0 : sum_cell_sources = dt*dot_product(s% dm(1:nz), s% energy_sources(1:nz))
1460 0 : sum_cell_others = dt*dot_product(s% dm(1:nz), s% energy_others(1:nz))
1461 0 : sum_cell_work = dt*dot_product(s% dm(1:nz), s% dwork_dm(1:nz))
1462 0 : sum_cell_detrb = dt*dot_product(s% dm(1:nz), s% detrbdt(1:nz))
1463 0 : sum_cell_dke = dt*dot_product(s% dm(1:nz), s% dkedt(1:nz))
1464 0 : sum_cell_dpe = dt*dot_product(s% dm(1:nz), s% dpedt(1:nz))
1465 0 : sum_cell_de = dt*dot_product(s% dm(1:nz), s% dedt(1:nz))
1466 : sum_cell_terms = &
1467 : - sum_cell_dL + sum_cell_sources + sum_cell_others - sum_cell_work &
1468 0 : - sum_cell_detrb - sum_cell_dke - sum_cell_dpe - sum_cell_de
1469 0 : sum_cell_terms = -sum_cell_terms ! to make it the same sign as sum_cell_ergs_error
1470 0 : sum_cell_ergs_error = sum(s% ergs_error(1:nz))
1471 :
1472 : expected_sum_cell_others = &
1473 : - total_energy_from_pre_mixing &
1474 : - s% total_energy_from_phase_separation &
1475 : - s% total_WD_sedimentation_heating &
1476 0 : - s% total_energy_from_diffusion
1477 : expected_sum_cell_sources = &
1478 : phase2_total_energy_from_mdot &
1479 : + s% total_nuclear_heating &
1480 : - s% total_non_nuc_neu_cooling &
1481 : + s% total_irradiation_heating &
1482 0 : + s% total_extra_heating
1483 :
1484 : !write(*,2) 'rel err sum all cell terms', s% model_number, &
1485 : ! (phase2_sources_and_sinks - &
1486 : ! (sum_cell_others + sum_cell_sources + sum_cell_dL + sum_cell_work))/s% total_energy
1487 0 : write(*,2) 'rel err sum_cell_others', s% model_number, &
1488 0 : (sum_cell_others - expected_sum_cell_others)/s% total_energy, &
1489 0 : sum_cell_others, expected_sum_cell_others
1490 0 : write(*,2) 'rel err sum_cell_sources', s% model_number, &
1491 0 : (sum_cell_sources - expected_sum_cell_sources)/s% total_energy, &
1492 0 : sum_cell_sources, expected_sum_cell_sources
1493 0 : write(*,2) 'rel err sum_cell_dL', s% model_number, &
1494 0 : (sum_cell_dL - total_radiation)/s% total_energy, sum_cell_dL, total_radiation
1495 0 : write(*,2) 'rel err sum_cell_work', s% model_number, &
1496 0 : (sum_cell_work - phase2_work)/s% total_energy, sum_cell_work, phase2_work
1497 0 : write(*,'(A)')
1498 :
1499 : diff_total_internal_energy = &
1500 0 : s% total_internal_energy_end - s% total_internal_energy_start
1501 : diff_total_gravitational_energy = &
1502 0 : s% total_gravitational_energy_end - s% total_gravitational_energy_start
1503 : diff_total_kinetic_energy = &
1504 0 : s% total_radial_kinetic_energy_end - s% total_radial_kinetic_energy_start
1505 : !diff_total_rotational_kinetic_energy = &
1506 : ! s% total_rotational_kinetic_energy_end - s% total_rotational_kinetic_energy_start
1507 : diff_total_turbulent_energy = &
1508 0 : s% total_turbulent_energy_end - s% total_turbulent_energy_start
1509 :
1510 0 : write(*,2) 'phase2 rel err sum_cell_de', s% model_number, &
1511 0 : (sum_cell_de - diff_total_internal_energy)/s% total_energy, &
1512 0 : sum_cell_de, diff_total_internal_energy
1513 0 : write(*,2) 'phase2 rel err sum_cell_dpe', s% model_number, &
1514 0 : (sum_cell_dpe - diff_total_gravitational_energy)/s% total_energy, &
1515 0 : sum_cell_dpe, diff_total_gravitational_energy
1516 0 : write(*,2) 'phase2 rel err sum_cell_dke', s% model_number, &
1517 0 : (sum_cell_dke - diff_total_kinetic_energy)/s% total_energy, &
1518 0 : sum_cell_dke, diff_total_kinetic_energy
1519 : !write(*,2) 'rel err ', s% model_number, &
1520 : ! ( - diff_total_rotational_kinetic_energy)/s% total_energy, &
1521 : ! , diff_total_rotational_kinetic_energy
1522 0 : write(*,2) 'rel err sum_cell_detrb', s% model_number, &
1523 0 : (sum_cell_detrb - diff_total_turbulent_energy)/s% total_energy, &
1524 0 : sum_cell_detrb, diff_total_turbulent_energy
1525 0 : write(*,'(A)')
1526 :
1527 :
1528 0 : write(*,2) 'expected rel sum_cell_ergs_error', s% model_number, &
1529 0 : sum_cell_ergs_error/s% total_energy, &
1530 0 : sum_cell_ergs_error, s% total_energy
1531 0 : write(*,2) 'actual rel err phase2_sources_and_sinks', s% model_number, &
1532 0 : (s% total_energy_end - (s% total_energy_start + phase2_sources_and_sinks))/s% total_energy
1533 0 : write(*,2) 'actual/expected', s% model_number, &
1534 0 : (s% total_energy_end - (s% total_energy_start + phase2_sources_and_sinks))/sum_cell_ergs_error
1535 0 : write(*,2) 'total rel_E_err', s% model_number, &
1536 0 : s% error_in_energy_conservation/s% total_energy, &
1537 0 : s% error_in_energy_conservation, s% total_energy
1538 0 : write(*,'(A)')
1539 : end if
1540 :
1541 0 : call mesa_error(__FILE__,__LINE__,'okay_energy_conservation')
1542 :
1543 : end if
1544 :
1545 :
1546 11 : if (is_bad_num(s% error_in_energy_conservation)) then
1547 0 : write(*,2) 's% error_in_energy_conservation', &
1548 0 : s% model_number, s% error_in_energy_conservation
1549 0 : write(*,2) 's% total_energy_end', &
1550 0 : s% model_number, s% total_energy_end
1551 0 : write(*,2) 's% total_energy_change_from_mdot', &
1552 0 : s% model_number, s% total_energy_change_from_mdot
1553 0 : write(*,2) 's% total_energy_start', &
1554 0 : s% model_number, s% total_energy_start
1555 0 : write(*,2) 's% total_energy_sources_and_sinks', &
1556 0 : s% model_number, s% total_energy_sources_and_sinks
1557 0 : write(*,2) 's% total_nuclear_heating', s% model_number, s% total_nuclear_heating
1558 0 : write(*,2) 's% total_non_nuc_neu_cooling', s% model_number, s% total_non_nuc_neu_cooling
1559 0 : write(*,2) 's% total_irradiation_heating', s% model_number, s% total_irradiation_heating
1560 0 : write(*,2) 's% total_extra_heating', s% model_number, s% total_extra_heating
1561 0 : write(*,2) 'dt*L_surf', s% model_number, dt*L_surf
1562 0 : write(*,2) 'dt*L_center', s% model_number, dt*s% L_center
1563 0 : write(*,2) 'L_surf', s% model_number, L_surf
1564 0 : write(*,2) 's% Fr(1)', s% model_number, s% Fr(1)
1565 0 : write(*,2) 's% Lc(1)', s% model_number, s% Lc(1)
1566 0 : write(*,2) 's% Lt(1)', s% model_number, s% Lt(1)
1567 0 : write(*,2) 'sum L', s% model_number, s% Fr(1)*pi4*s% r(1)*s% r(1)+s% Lc(1)+s% Lt(1)
1568 0 : okay_energy_conservation = .false.
1569 0 : call mesa_error(__FILE__,__LINE__,'okay_energy_conservation')
1570 0 : return
1571 : end if
1572 :
1573 11 : if (is_bad_num(s% cumulative_energy_error)) then
1574 0 : write(*,2) 's% cumulative_energy_error', &
1575 0 : s% model_number, s% cumulative_energy_error
1576 0 : write(*,2) 's% cumulative_energy_error_old', &
1577 0 : s% model_number, s% cumulative_energy_error_old
1578 0 : write(*,2) 's% error_in_energy_conservation', &
1579 0 : s% model_number, s% error_in_energy_conservation
1580 0 : write(*,2) 's% total_energy_sources_and_sinks', &
1581 0 : s% model_number, s% total_energy_sources_and_sinks
1582 0 : write(*,2) 's% total_nuclear_heating', s% model_number, s% total_nuclear_heating
1583 0 : write(*,2) 's% total_non_nuc_neu_cooling', s% model_number, s% total_non_nuc_neu_cooling
1584 0 : write(*,2) 's% total_irradiation_heating', s% model_number, s% total_irradiation_heating
1585 0 : write(*,2) 's% total_extra_heating', s% model_number, s% total_extra_heating
1586 0 : write(*,2) 's% work_inward_at_center', s% model_number, s% work_inward_at_center
1587 0 : write(*,2) 's% work_outward_at_surface', s% model_number, s% work_outward_at_surface
1588 0 : write(*,2) 's% L_center', s% model_number, s% L_center
1589 0 : okay_energy_conservation = .false.
1590 0 : call mesa_error(__FILE__,__LINE__,'okay_energy_conservation')
1591 0 : return
1592 : end if
1593 :
1594 11 : okay_energy_conservation = .true.
1595 :
1596 11 : end function okay_energy_conservation
1597 :
1598 : end function do_step_part2
1599 :
1600 :
1601 : subroutine debug(str, s)
1602 : use chem_def
1603 : character (len=*), intent(in) :: str
1604 : type(star_info), pointer :: s
1605 : integer :: k, j
1606 : include 'formats'
1607 :
1608 : return
1609 :
1610 : if (.not. s% rotation_flag) return
1611 : k = 1
1612 : write(*,3) trim(str) // ' s% omega(k)', k, s% model_number, s% omega(k)
1613 : return
1614 : j = 2
1615 : !do j=1,1 !s% species
1616 : if (.true. .or. s% xa(j,k) > 1d-9) &
1617 : write(*,1) trim(str) // ' xin(net_iso(i' // &
1618 : trim(chem_isos% name(s% chem_id(j))) // '))= ', &
1619 : s% xa(j,k), s% abar(k)
1620 : !end do
1621 : end subroutine debug
1622 :
1623 :
1624 11 : subroutine check_for_extra_heat(s, ierr)
1625 : use hydro_vars, only: set_vars
1626 : use utils_lib, only: is_bad
1627 : type (star_info), pointer :: s
1628 : integer, intent(out) :: ierr
1629 :
1630 11 : real(dp) :: start_time, end_time, left_to_inject, &
1631 11 : q00, qp1, qmin, qmax, qtop, qbot, extra, dt, &
1632 11 : target_injection_time, target_injection_ergs, &
1633 11 : kap_gamma, tau_gamma_sum
1634 : integer :: k, nz, k1
1635 :
1636 : include 'formats'
1637 :
1638 11 : ierr = 0
1639 11 : if (s% use_other_energy_implicit) return
1640 :
1641 11 : nz = s% nz
1642 11 : dt = s% dt
1643 13098 : do k=1,nz
1644 13098 : s% extra_heat(k) = s% extra_power_source
1645 : end do
1646 :
1647 11 : if (s% use_other_energy) then
1648 0 : call s% other_energy(s% id, ierr)
1649 0 : if (ierr /= 0) then
1650 0 : if (s% report_ierr) &
1651 0 : write(*, *) 'check_for_extra_heat: other_energy returned ierr', ierr
1652 0 : return
1653 : end if
1654 11 : else if (s% inject_uniform_extra_heat /= 0) then
1655 0 : qp1 = 0d0
1656 0 : qmin = s% min_q_for_uniform_extra_heat
1657 0 : qmax = s% max_q_for_uniform_extra_heat
1658 0 : extra = s% inject_uniform_extra_heat
1659 0 : do k=nz,1,-1
1660 0 : q00 = s% q(k)
1661 0 : if (qp1 >= qmin .and. q00 <= qmax) then ! all inside of region
1662 0 : s% extra_heat(k) = s% extra_heat(k) + extra
1663 : else
1664 0 : qtop = min(q00, qmax)
1665 0 : qbot = max(qp1, qmin)
1666 0 : if (qtop > qbot) then ! overlaps region
1667 0 : s% extra_heat(k) = s% extra_heat(k) + extra*(qtop - qbot)/s% dq(k)
1668 : end if
1669 : end if
1670 0 : qp1 = q00
1671 : end do
1672 0 : s% need_to_setvars = .true.
1673 11 : else if (s% nonlocal_NiCo_decay_heat) then
1674 0 : kap_gamma = s% nonlocal_NiCo_kap_gamma
1675 0 : do k1=1,nz
1676 : tau_gamma_sum = 0
1677 0 : do k=k1,1,-1 ! move eps_nuc outward from k1 to extra_heat at k
1678 : tau_gamma_sum = tau_gamma_sum + &
1679 0 : kap_gamma*s% dm(k)/(pi4*s% rmid(k)*s% rmid(k))
1680 0 : if (tau_gamma_sum >= s% dtau_gamma_NiCo_decay_heat) then
1681 : s% extra_heat(k) = s% extra_heat(k) + &
1682 0 : s% eps_nuc(k1)*s% dm(k1)/s% dm(k)
1683 0 : exit
1684 : end if
1685 : end do
1686 : end do
1687 0 : s% need_to_setvars = .true.
1688 : end if
1689 :
1690 : if (s% inject_until_reach_model_with_total_energy <= s% total_energy_initial &
1691 11 : .or. dt <= 0d0 .or. s% total_mass_for_inject_extra_ergs_sec <= 0d0) return
1692 :
1693 0 : start_time = s% start_time_for_inject_extra_ergs_sec
1694 0 : if (s% time < start_time) return
1695 :
1696 0 : if (s% duration_for_inject_extra_ergs_sec > 0) then
1697 0 : end_time = start_time + s% duration_for_inject_extra_ergs_sec
1698 0 : else if (s% max_age_in_seconds > 0) then
1699 : end_time = s% max_age_in_seconds
1700 0 : else if (s% max_age_in_days > 0) then
1701 0 : end_time = s% max_age_in_days*secday
1702 : else
1703 0 : end_time = s% max_age*secyer
1704 : end if
1705 0 : if (s% time_old > end_time) return
1706 :
1707 : target_injection_ergs = &
1708 0 : s% inject_until_reach_model_with_total_energy - s% total_energy_initial
1709 0 : target_injection_time = end_time - start_time
1710 0 : s% inject_extra_ergs_sec = target_injection_ergs/target_injection_time
1711 : left_to_inject = &
1712 0 : s% inject_until_reach_model_with_total_energy - s% total_energy_old
1713 0 : qp1 = 0d0
1714 0 : qmin = max(0d0, Msun*s% base_of_inject_extra_ergs_sec - s% M_center)/s% xmstar
1715 0 : qmax = min(1d0, qmin + Msun*s% total_mass_for_inject_extra_ergs_sec/s% xmstar)
1716 0 : extra = s% inject_extra_ergs_sec/(s% xmstar*(qmax - qmin))
1717 0 : if (s% time > end_time .or. s% time_old < start_time) then
1718 0 : extra = extra*(min(s% time, end_time) - max(s% time_old, start_time))/dt
1719 : end if
1720 0 : if (left_to_inject < extra*dt*s% xmstar*(qmax - qmin)) then
1721 0 : extra = left_to_inject/(dt*s% xmstar*(qmax - qmin))
1722 : end if
1723 0 : if (is_bad(extra)) then
1724 0 : write(*,2) 'extra', s% model_number, extra
1725 0 : write(*,2) 'left_to_inject', s% model_number, left_to_inject
1726 0 : write(*,2) 's% total_energy_old', s% model_number, s% total_energy_old
1727 0 : call mesa_error(__FILE__,__LINE__,'check_for_extra_heat')
1728 : end if
1729 0 : do k=nz,1,-1
1730 0 : q00 = s% q(k)
1731 0 : if (qp1 >= qmin .and. q00 <= qmax) then ! all inside of region
1732 0 : s% extra_heat(k) = s% extra_heat(k) + extra
1733 : else
1734 0 : qtop = min(q00, qmax)
1735 0 : qbot = max(qp1, qmin)
1736 0 : if (qtop > qbot) then ! overlaps region
1737 0 : s% extra_heat(k) = s% extra_heat(k) + extra*(qtop - qbot)/s% dq(k)
1738 : end if
1739 : end if
1740 0 : qp1 = q00
1741 : end do
1742 :
1743 11 : end subroutine check_for_extra_heat
1744 :
1745 :
1746 21 : subroutine set_start_of_step_info(s, str, ierr)
1747 11 : use report, only: do_report
1748 : use hydro_vars, only: set_vars_if_needed
1749 : use turb_info, only: set_gradT_excess_alpha
1750 : use star_utils, only: min_dr_div_cs, get_remnant_mass, &
1751 : total_angular_momentum, eval_Ledd, set_luminosity_by_category
1752 :
1753 : type (star_info), pointer :: s
1754 : character (len=*), intent(in) :: str
1755 : integer, intent(out) :: ierr
1756 :
1757 : logical :: trace
1758 : integer :: nz
1759 :
1760 : include 'formats'
1761 :
1762 21 : ierr = 0
1763 21 : trace = s% trace_evolve
1764 21 : nz = s% nz
1765 :
1766 21 : if (.not. s% RSP_flag) then
1767 21 : call set_vars_if_needed(s, s% dt, str, ierr)
1768 21 : if (failed('set_vars_if_needed')) return
1769 226659 : s% edv(1:s% species, 1:s% nz) = 0
1770 21 : call set_luminosity_by_category(s)
1771 21 : s% total_angular_momentum = total_angular_momentum(s)
1772 21 : call do_report(s, ierr)
1773 21 : if (failed('do_report ierr')) return
1774 : end if
1775 :
1776 : ! save a few things from start of step that will need later
1777 21 : if (s% rotation_flag) then
1778 0 : s% surf_r_equatorial = s% r_equatorial(1)
1779 : else
1780 21 : s% surf_r_equatorial = s% r(1)
1781 : end if
1782 21 : s% starting_T_center = s% T(nz)
1783 21 : s% surf_csound = s% csound(1)
1784 21 : s% surf_rho = s% rho(1)
1785 21 : s% prev_Ledd = eval_Ledd(s,ierr)
1786 21 : if (failed('eval_Ledd ierr')) return
1787 :
1788 42 : if (.not. s% RSP_flag) then
1789 21 : call set_gradT_excess_alpha(s, ierr)
1790 21 : if (failed('set_gradT_excess_alpha ierr')) return
1791 : end if
1792 :
1793 : contains
1794 :
1795 84 : logical function failed(str)
1796 : character (len=*), intent(in) :: str
1797 84 : if (ierr == 0) then
1798 84 : failed = .false.
1799 : return
1800 : end if
1801 0 : failed = .true.
1802 0 : if (s% report_ierr) write(*, *) 'set_start_of_step_info: ' // trim(str)
1803 0 : s% result_reason = nonzero_ierr
1804 21 : end function failed
1805 :
1806 : end subroutine set_start_of_step_info
1807 :
1808 :
1809 11 : integer function prepare_for_new_step(s)
1810 : use evolve_support, only: new_generation
1811 : use chem_def
1812 : use star_utils, only: use_xh_to_set_rho_to_dm_div_dV, set_phot_info
1813 : use hydro_vars, only: set_vars_if_needed
1814 :
1815 : type (star_info), pointer :: s
1816 :
1817 : integer :: ierr, k
1818 11 : real(dp) :: force_timestep_min, force_timestep
1819 : logical :: trace
1820 :
1821 : include 'formats'
1822 :
1823 11 : ierr = 0
1824 11 : trace = s% trace_evolve
1825 :
1826 11 : prepare_for_new_step = keep_going
1827 :
1828 11 : if (s% dt_next <= 0) then
1829 0 : write(*, *) 's% dt_next', s% dt_next
1830 0 : prepare_for_new_step = terminate
1831 : if ((s% time >= s% max_age*secyer .and. s% max_age > 0) .or. &
1832 0 : (s% time >= s% max_age_in_days*secday .and. s% max_age_in_days > 0) .or. &
1833 : (s% time >= s% max_age_in_seconds .and. s% max_age_in_seconds > 0)) then
1834 0 : s% result_reason = result_reason_normal
1835 0 : s% termination_code = t_max_age
1836 : else
1837 0 : s% result_reason = dt_is_zero
1838 0 : s% termination_code = t_dt_is_zero
1839 : end if
1840 0 : return
1841 : end if
1842 :
1843 11 : if (s% dt_next < s% min_timestep_limit) then
1844 0 : write(*, *) 's% dt_next', s% dt_next
1845 0 : write(*, *) 's% min_timestep_limit', s% min_timestep_limit
1846 : prepare_for_new_step = terminate
1847 0 : s% termination_code = t_min_timestep_limit
1848 0 : s% result_reason = timestep_limits
1849 0 : return
1850 : end if
1851 :
1852 11 : if (s% set_rho_to_dm_div_dV .and. .not. (s% u_flag .or. s% RSP_flag)) then
1853 : call use_xh_to_set_rho_to_dm_div_dV(s,ierr)
1854 11 : if (failed('use_xh_to_set_rho_to_dm_div_dV ierr')) return
1855 : end if
1856 :
1857 11 : if (.not. s% RSP_flag) then ! store mesh info for following step eps_mdot
1858 14570 : do k=1, s% nz
1859 131031 : s% prev_mesh_xa(:,k) = s% xa(:,k)
1860 72795 : s% prev_mesh_xh(:,k) = s% xh(:,k)
1861 14559 : s% prev_mesh_j_rot(k) = s% j_rot(k)
1862 14559 : s% prev_mesh_omega(k) = s% omega(k)
1863 14559 : s% prev_mesh_dq(k) = s% dq(k)
1864 14559 : s% prev_mesh_mlt_vc(k) = s% mlt_vc(k)
1865 14570 : s% prev_mesh_species_or_nvar_hydro_changed = .false.
1866 : end do
1867 11 : s% prev_mesh_nz = s% nz
1868 : ! also store ST information for time smoothing
1869 11 : if (s% have_ST_start_info) then
1870 0 : do k=1, s% nz
1871 0 : s% prev_mesh_D_ST_start(k) = s% D_ST_start(k)
1872 0 : s% prev_mesh_nu_ST_start(k) = s% nu_ST_start(k)
1873 : end do
1874 0 : s% prev_mesh_have_ST_start_info = .true.
1875 : else
1876 11 : s% prev_mesh_have_ST_start_info = .false.
1877 : end if
1878 : end if
1879 :
1880 11 : if (s% okay_to_remesh) then
1881 11 : if (s% rsp_flag .or. .not. s% doing_first_model_of_run) then
1882 : call set_start_of_step_info(s, 'before do_mesh', ierr)
1883 10 : if (failed('set_start_of_step_info ierr')) return
1884 10 : prepare_for_new_step = do_mesh(s) ! sets s% need_to_setvars = .true. if changes anything
1885 10 : if (prepare_for_new_step /= keep_going) return
1886 : end if
1887 : end if
1888 :
1889 11 : call set_vars_if_needed(s, s% dt_next, 'prepare_for_new_step', ierr)
1890 11 : if (failed('set_vars_if_needed ierr')) return
1891 11 : call set_phot_info(s) ! this sets Teff and other stuff
1892 :
1893 : call new_generation(s, ierr)
1894 11 : if (failed('new_generation ierr')) return
1895 11 : s% generations = 2
1896 :
1897 11 : if ((s% time + s% dt_next) > s% max_age*secyer .and. s% max_age > 0) then
1898 0 : s% dt_next = max(0d0, s% max_age*secyer - s% time)
1899 0 : if ( s% dt_next == 0d0 ) then
1900 0 : write(*,*) 'WARNING: max_age reached'
1901 0 : write(*,1) 's% max_age*secyer', s% max_age*secyer
1902 0 : write(*,1) 's% time', s% time
1903 0 : write(*,1) 's% max_age', s% max_age
1904 : end if
1905 : else if ((s% time + s% dt_next) > s% max_age_in_seconds &
1906 11 : .and. s% max_age_in_seconds > 0) then
1907 0 : s% dt_next = max(0d0, s% max_age_in_seconds - s% time)
1908 : else if ((s% time + s% dt_next) > s% max_age_in_days*secday &
1909 11 : .and. s% max_age_in_days > 0) then
1910 0 : s% dt_next = max(0d0, s% max_age_in_days*secday - s% time)
1911 : end if
1912 :
1913 11 : s% dt = s% dt_next
1914 :
1915 11 : force_timestep_min = s% force_timestep_min
1916 11 : if (force_timestep_min == 0) &
1917 11 : force_timestep_min = secyer*s% force_timestep_min_years
1918 11 : if (s% dt < force_timestep_min) then
1919 0 : s% dt = min(s% dt*s% force_timestep_min_factor, force_timestep_min)
1920 0 : write(*,2) 'force increase in timestep', s% model_number, s% dt
1921 : end if
1922 :
1923 11 : force_timestep = s% force_timestep
1924 11 : if (force_timestep == 0) &
1925 11 : force_timestep = secyer*s% force_timestep_years
1926 11 : if (force_timestep > 0) s% dt = force_timestep
1927 :
1928 11 : s% dt_start = s% dt
1929 11 : s% time_step = s% dt/secyer
1930 :
1931 11 : if (is_bad(s% dt)) then
1932 0 : write(*,1) 's% dt', s% dt
1933 0 : call mesa_error(__FILE__,__LINE__,'prepare_for_new_step')
1934 : end if
1935 :
1936 11 : s% retry_cnt = 0
1937 11 : s% redo_cnt = 0
1938 :
1939 11 : s% need_to_save_profiles_now = .false.
1940 11 : s% need_to_update_history_now = s% doing_first_model_of_run
1941 :
1942 : contains
1943 :
1944 43 : logical function failed(str)
1945 : character (len=*), intent(in) :: str
1946 43 : if (ierr == 0) then
1947 43 : failed = .false.
1948 : return
1949 : end if
1950 0 : failed = .true.
1951 0 : prepare_for_new_step = terminate
1952 0 : if (s% report_ierr) write(*, *) 'prepare_for_new_step: ' // trim(str)
1953 0 : s% result_reason = nonzero_ierr
1954 11 : end function failed
1955 :
1956 : end function prepare_for_new_step
1957 :
1958 :
1959 10 : integer function do_mesh(s)
1960 : use adjust_mesh, only: remesh
1961 : use adjust_mesh_split_merge, only: remesh_split_merge
1962 : use star_utils, only: start_time, update_time
1963 : type (star_info), pointer :: s
1964 : integer(i8) :: time0
1965 10 : real(dp) :: total
1966 : include 'formats'
1967 10 : do_mesh = keep_going
1968 10 : if (.not. s% okay_to_remesh) return
1969 : if (s% restore_mesh_on_retry &
1970 10 : .and. s% model_number_for_last_retry > s% model_number - s% num_steps_to_hold_mesh_after_retry) return
1971 10 : s% need_to_setvars = .true.
1972 10 : if (s% doing_timing) call start_time(s, time0, total)
1973 10 : if (s% use_split_merge_amr) then
1974 0 : do_mesh = remesh_split_merge(s) ! sets s% need_to_setvars = .true. if changes anything
1975 0 : if (do_mesh /= keep_going .and. s% report_ierr) &
1976 0 : write(*, *) 'do_mesh: remesh_split_merge failed'
1977 10 : else if (.not. s% rsp_flag) then
1978 10 : do_mesh = remesh(s) ! sets s% need_to_setvars = .true. if changes anything
1979 10 : if (do_mesh /= keep_going .and. s% report_ierr) &
1980 0 : write(*, *) 'do_mesh: remesh failed'
1981 : end if
1982 10 : if (s% doing_timing) call update_time(s, time0, total, s% time_remesh)
1983 10 : if (do_mesh /= keep_going) then
1984 0 : s% result_reason = adjust_mesh_failed
1985 0 : return
1986 : end if
1987 10 : end function do_mesh
1988 :
1989 :
1990 11 : integer function prepare_for_new_try(s)
1991 : ! return keep_going, terminate, or retry
1992 : ! if don't return keep_going, then set result_reason to say why.
1993 10 : use net_lib, only: clean_up_fractions
1994 : use net, only: get_screening_mode
1995 : use hydro_rotation, only: use_xh_to_update_i_rot
1996 :
1997 : type (star_info), pointer :: s
1998 :
1999 : integer :: ierr
2000 11 : real(dp) :: screening
2001 :
2002 : include 'formats'
2003 :
2004 : ierr = 0
2005 11 : s% result_reason = result_reason_normal
2006 11 : s% termination_code = 0
2007 11 : s% solver_iter = 0
2008 11 : s% solver_adjust_iter = 0
2009 :
2010 11 : s% model_number = s% model_number_old + 1
2011 :
2012 11 : prepare_for_new_try = keep_going
2013 :
2014 : s% result_reason = result_reason_normal
2015 : s% model_number = s% model_number_old + 1
2016 : s% termination_code = 0
2017 : s% solver_iter = 0
2018 : s% solver_adjust_iter = 0
2019 :
2020 11 : if (.not. s% RSP_flag) then
2021 :
2022 11 : screening = get_screening_mode(s,ierr)
2023 11 : if (ierr /= 0) then
2024 0 : write(*,*) 'bad value for screening_mode ' // trim(s% screening_mode)
2025 0 : prepare_for_new_try = terminate
2026 0 : s% termination_code = t_failed_prepare_for_new_try
2027 0 : return
2028 : end if
2029 :
2030 : end if
2031 :
2032 11 : end function prepare_for_new_try
2033 :
2034 :
2035 11 : integer function pick_next_timestep(id)
2036 : ! determine what we want for the next timestep
2037 : ! if don't return keep_going, then set result_reason to say why.
2038 11 : use timestep, only: timestep_controller
2039 : integer, intent(in) :: id
2040 : integer :: ierr
2041 : type (star_info), pointer :: s
2042 : integer :: i, j, n
2043 11 : real(dp) :: max_timestep, remaining_years, min_max, prev_max_years
2044 : include 'formats'
2045 :
2046 11 : pick_next_timestep = terminate
2047 11 : call get_star_ptr(id, s, ierr)
2048 11 : if (ierr /= 0) return
2049 :
2050 11 : if (s% RSP_flag) then
2051 0 : pick_next_timestep = keep_going
2052 0 : s% dt_next = s% RSP_dt
2053 0 : s% dt_next_unclipped = s% dt_next
2054 0 : s% why_Tlim = Tlim_max_timestep
2055 0 : return
2056 : end if
2057 :
2058 11 : if (s% trace_evolve) write(*,'(/,a)') 'pick_next_timestep'
2059 :
2060 11 : if (s% max_years_for_timestep > 0) then
2061 0 : max_timestep = secyer*s% max_years_for_timestep
2062 0 : if (s% max_timestep > 0 .and. s% max_timestep < max_timestep) &
2063 0 : max_timestep = s% max_timestep
2064 : else
2065 11 : max_timestep = s% max_timestep
2066 : end if
2067 :
2068 11 : pick_next_timestep = timestep_controller(s, max_timestep)
2069 11 : if (pick_next_timestep /= keep_going) then
2070 0 : if (s% trace_evolve) &
2071 0 : write(*,*) 'pick_next_timestep: timestep_controller /= keep_going'
2072 0 : return
2073 : end if
2074 :
2075 11 : s% dt_next_unclipped = s% dt_next
2076 : ! write out the unclipped timestep in saved models
2077 11 : if (s% time < 0 .and. s% time + s% dt_next > 0) then
2078 0 : s% dt_next = -s% time
2079 11 : else if ((s% time + s% dt_next) > s% max_age*secyer .and. s% max_age > 0) then
2080 2 : s% dt_next = max(0d0, s% max_age*secyer - s% time)
2081 : else if ((s% time + s% dt_next) > s% max_age_in_days*secday &
2082 9 : .and. s% max_age_in_days > 0) then
2083 0 : s% dt_next = max(0d0, s% max_age_in_days*secday - s% time)
2084 : else if ((s% time + s% dt_next) > s% max_age_in_seconds &
2085 9 : .and. s% max_age_in_seconds > 0) then
2086 0 : s% dt_next = max(0d0, s% max_age_in_seconds - s% time)
2087 9 : else if (s% num_adjusted_dt_steps_before_max_age > 0 .and. &
2088 : s% max_years_for_timestep > 0) then
2089 0 : if (s% max_age > 0) then
2090 0 : remaining_years = s% max_age - s% star_age
2091 0 : else if (s% max_age_in_days > 0) then
2092 0 : remaining_years = (s% max_age_in_days*secday - s% time)/secyer
2093 0 : else if (s% max_age_in_seconds > 0) then
2094 0 : remaining_years = (s% max_age_in_seconds - s% time)/secyer
2095 : else
2096 0 : remaining_years = 1d99
2097 : end if
2098 0 : if (s% using_revised_max_yr_dt) &
2099 0 : s% max_years_for_timestep = s% revised_max_yr_dt
2100 0 : n = floor(remaining_years/s% max_years_for_timestep + 1d-6)
2101 0 : j = s% num_adjusted_dt_steps_before_max_age
2102 0 : if (remaining_years <= s% max_years_for_timestep) then
2103 0 : s% max_years_for_timestep = remaining_years
2104 0 : s% using_revised_max_yr_dt = .true.
2105 0 : s% revised_max_yr_dt = s% max_years_for_timestep
2106 0 : s% dt_next = s% max_years_for_timestep*secyer
2107 0 : write(*,3) 'remaining steps and years until max age', &
2108 0 : s% model_number, 1, remaining_years
2109 0 : else if (n <= j) then
2110 0 : prev_max_years = s% max_years_for_timestep
2111 0 : i = floor(remaining_years/s% dt_years_for_steps_before_max_age + 1d-6)
2112 0 : if ((i+1d-9)*s% dt_years_for_steps_before_max_age < remaining_years) then
2113 0 : s% max_years_for_timestep = remaining_years/(i+1)
2114 : else
2115 0 : s% max_years_for_timestep = remaining_years/i
2116 : end if
2117 0 : min_max = prev_max_years*s% reduction_factor_for_max_timestep
2118 0 : if (s% max_years_for_timestep < min_max) &
2119 0 : s% max_years_for_timestep = min_max
2120 0 : if (.not. s% using_revised_max_yr_dt) then
2121 0 : s% using_revised_max_yr_dt = .true.
2122 0 : write(*,2) 'begin reducing max timestep prior to max age', &
2123 0 : s% model_number, remaining_years
2124 0 : else if (s% revised_max_yr_dt > s% max_years_for_timestep) then
2125 0 : write(*,2) 'reducing max timestep prior to max age', &
2126 0 : s% model_number, remaining_years
2127 0 : else if (s% max_years_for_timestep <= &
2128 : s% dt_years_for_steps_before_max_age) then
2129 0 : i = floor(remaining_years/s% max_years_for_timestep + 1d-6)
2130 0 : write(*,3) 'remaining steps and years until max age', &
2131 0 : s% model_number, i, remaining_years
2132 : else
2133 0 : write(*,2) 'remaining_years until max age', &
2134 0 : s% model_number, remaining_years
2135 : end if
2136 0 : s% revised_max_yr_dt = s% max_years_for_timestep
2137 0 : if (s% dt_next/secyer > s% max_years_for_timestep) &
2138 0 : s% dt_next = s% max_years_for_timestep*secyer
2139 : end if
2140 :
2141 : end if
2142 :
2143 11 : end function pick_next_timestep
2144 :
2145 :
2146 0 : integer function prepare_to_redo(id)
2147 11 : use evolve_support, only: set_current_to_old
2148 : integer, intent(in) :: id
2149 : type (star_info), pointer :: s
2150 : integer :: ierr
2151 : include 'formats'
2152 : ierr = 0
2153 0 : call get_star_ptr(id, s, ierr)
2154 0 : if (ierr /= 0) then
2155 0 : prepare_to_redo = terminate
2156 : return
2157 : end if
2158 0 : s% redo_cnt = s% redo_cnt + 1
2159 0 : if (s% redo_limit > 0 .and. s% redo_cnt > s% redo_limit) then
2160 0 : write(*,2) 'redo_cnt', s% redo_cnt
2161 0 : write(*,2) 'redo_limit', s% redo_limit
2162 0 : call report_problems(s, '-- too many redos')
2163 0 : s% termination_code = t_redo_limit
2164 0 : prepare_to_redo = terminate
2165 0 : return
2166 : end if
2167 0 : prepare_to_redo = keep_going
2168 0 : if (s% trace_evolve) write(*,'(/,a)') 'prepare_to_redo'
2169 0 : call set_current_to_old(s)
2170 0 : s% need_to_setvars = .true.
2171 0 : end function prepare_to_redo
2172 :
2173 :
2174 0 : integer function prepare_to_retry(id)
2175 0 : use evolve_support, only: set_current_to_old
2176 : integer, intent(in) :: id
2177 0 : real(dp) :: retry_factor
2178 : type (star_info), pointer :: s
2179 : integer :: ierr, k
2180 : include 'formats'
2181 :
2182 0 : call get_star_ptr(id, s, ierr)
2183 0 : if (ierr /= 0) then
2184 0 : prepare_to_retry = terminate
2185 : return
2186 : end if
2187 :
2188 0 : s% need_to_setvars = .true.
2189 :
2190 0 : if (s% restore_mesh_on_retry .and. .not. s% RSP_flag) then
2191 0 : if (.not. s% prev_mesh_species_or_nvar_hydro_changed) then
2192 0 : do k=1, s% prev_mesh_nz
2193 0 : s% xh_old(:,k) = s% prev_mesh_xh(:,k)
2194 0 : s% xa_old(:,k) = s% prev_mesh_xa(:,k)
2195 0 : s% dq_old(k) = s% prev_mesh_dq(k)
2196 0 : s% omega_old(k) = s% prev_mesh_omega(k)
2197 0 : s% j_rot_old(k) = s% prev_mesh_j_rot(k)
2198 0 : s% mlt_vc_old(k) = s% prev_mesh_mlt_vc(k)
2199 : end do
2200 0 : call normalize_dqs(s, s% prev_mesh_nz, s% dq_old, ierr)
2201 0 : if (ierr /= 0) then
2202 0 : prepare_to_retry = terminate
2203 : return
2204 : end if
2205 0 : s% nz_old = s% prev_mesh_nz
2206 : end if
2207 : end if
2208 :
2209 0 : if (s% trace_evolve) write(*,'(/,a)') 'prepare_to_retry'
2210 :
2211 0 : s% retry_cnt = s% retry_cnt + 1
2212 0 : if (s% retry_limit > 0 .and. s% retry_cnt > s% retry_limit) then
2213 0 : s% dt_start = sqrt(s% dt*s% dt_start)
2214 0 : prepare_to_retry = terminate
2215 0 : return
2216 : end if
2217 :
2218 0 : prepare_to_retry = keep_going
2219 :
2220 0 : retry_factor = s% timestep_factor_for_retries
2221 0 : s% dt = s% dt*retry_factor
2222 0 : s% time_step = s% dt/secyer
2223 0 : if (len_trim(s% retry_message) > 0) then
2224 0 : if (s% retry_message_k > 0) then
2225 0 : write(*,'(a, 2i8)') ' retry: ' // trim(s% retry_message), s% retry_message_k, s% model_number
2226 : else
2227 0 : write(*,'(a, i8)') ' retry: ' // trim(s% retry_message), s% model_number
2228 : end if
2229 : else
2230 0 : write(*,'(a, i8)') ' retry', s% model_number
2231 : !if (.true.) call mesa_error(__FILE__,__LINE__,'failed to set retry_message')
2232 : end if
2233 0 : s% retry_message_k = 0
2234 0 : if (s% report_ierr) &
2235 0 : write(*,'(a50,2i6,3f16.6)') 'retry log10(dt/yr), log10(dt), retry_factor', &
2236 0 : s% retry_cnt, s% model_number, log10(s% dt*retry_factor/secyer), &
2237 0 : log10(s% dt*retry_factor), retry_factor
2238 0 : if (s% dt <= max(s% min_timestep_limit,0d0)) then
2239 0 : write(*,1) 'dt', s% dt
2240 0 : write(*,1) 'min_timestep_limit', s% min_timestep_limit
2241 0 : call report_problems(s, 'dt < min_timestep_limit')
2242 0 : prepare_to_retry = terminate
2243 0 : s% termination_code = t_min_timestep_limit
2244 0 : s% result_reason = timestep_limits
2245 0 : return
2246 : end if
2247 :
2248 0 : if (s% max_years_for_timestep > 0) &
2249 0 : s% max_timestep = secyer*s% max_years_for_timestep
2250 0 : if (s% max_timestep > 0 .and. s% dt > s% max_timestep) then
2251 0 : s% dt = s% max_timestep
2252 0 : s% time_step = s% dt/secyer
2253 : end if
2254 :
2255 0 : call set_current_to_old(s)
2256 :
2257 0 : s% num_retries = s% num_retries+1
2258 0 : if (s% num_retries > s% max_number_retries .and. s% max_number_retries >= 0) then
2259 0 : write(*,2) 'num_retries', s% num_retries
2260 0 : write(*,2) 'max_number_retries', s% max_number_retries
2261 0 : call report_problems(s, '-- too many retries')
2262 0 : s% termination_code = t_max_number_retries
2263 0 : prepare_to_retry = terminate; return
2264 : end if
2265 :
2266 0 : s% model_number_for_last_retry = s% model_number
2267 0 : if (s% why_Tlim == Tlim_neg_X) then
2268 : s% timestep_hold = s% model_number + &
2269 0 : max(s% retry_hold, s% neg_mass_fraction_hold)
2270 : else
2271 0 : s% timestep_hold = s% model_number + s% retry_hold
2272 : end if
2273 0 : s% why_Tlim = Tlim_retry
2274 :
2275 0 : end function prepare_to_retry
2276 :
2277 :
2278 0 : subroutine report_problems(s,str)
2279 : type (star_info), pointer :: s
2280 : character (len=*), intent(in) :: str
2281 0 : write(*,*) 'stopping because of problems ' // trim(str)
2282 0 : end subroutine report_problems
2283 :
2284 :
2285 11 : integer function finish_step(id, ierr)
2286 : ! returns keep_going or terminate
2287 : ! if don't return keep_going, then set result_reason to say why.
2288 : use evolve_support, only: output
2289 : use profile, only: do_save_profiles
2290 : use history, only: write_history_info
2291 : use utils_lib, only: free_iounit, number_iounits_allocated
2292 : use alloc, only: size_work_arrays
2293 :
2294 : integer, intent(in) :: id
2295 : integer, intent(out) :: ierr
2296 :
2297 : type (star_info), pointer :: s
2298 : integer, parameter :: nvals = 1, n_ivals = 0
2299 : integer :: nz, current_num_iounits_in_use, prev_num_iounits_in_use
2300 : logical :: trace, will_do_photo
2301 :
2302 : include 'formats'
2303 :
2304 11 : finish_step = terminate
2305 :
2306 11 : call get_star_ptr(id, s, ierr)
2307 11 : if (ierr /= 0) return
2308 :
2309 11 : nz = s% nz
2310 11 : trace = s% trace_evolve
2311 11 : prev_num_iounits_in_use = number_iounits_allocated()
2312 :
2313 11 : finish_step = keep_going
2314 11 : s% result_reason = result_reason_normal
2315 :
2316 11 : if (s% need_to_save_profiles_now .and. s% write_profiles_flag) then
2317 1 : call do_save_profiles(s, ierr)
2318 1 : s% need_to_save_profiles_now = .false.
2319 : end if
2320 :
2321 11 : call check(1)
2322 :
2323 11 : if (s% need_to_update_history_now .and. s% do_history_file) then
2324 :
2325 : call write_history_info( &
2326 4 : s, ierr)
2327 4 : if (ierr /= 0) then
2328 0 : finish_step = terminate
2329 0 : if (s% report_ierr) write(*, *) 'finish_step: write_history_info ierr', ierr
2330 0 : s% result_reason = nonzero_ierr
2331 0 : return
2332 : end if
2333 4 : s% need_to_update_history_now = .false.
2334 : end if
2335 :
2336 11 : call check(2)
2337 :
2338 11 : will_do_photo = .false.
2339 11 : if (.not. s% doing_relax) then
2340 11 : if(s% photo_interval > 0) then
2341 11 : if(mod(s% model_number, s% photo_interval) == 0) will_do_photo = .true.
2342 : end if
2343 11 : if(s% solver_save_photo_call_number > 0)then
2344 0 : if(s% solver_call_number == s% solver_save_photo_call_number - 1) will_do_photo = .true.
2345 : end if
2346 : end if
2347 :
2348 11 : if (will_do_photo) then
2349 :
2350 0 : call output(id, ierr)
2351 :
2352 0 : if (ierr /= 0) then
2353 0 : finish_step = terminate
2354 0 : if (s% report_ierr) write(*, *) 'finish_step: output ierr', ierr
2355 0 : s% result_reason = nonzero_ierr
2356 0 : return
2357 : end if
2358 :
2359 : end if
2360 :
2361 11 : call check(3)
2362 :
2363 11 : s% screening_mode_value = -1 ! force a new lookup for next step
2364 11 : s% doing_first_model_of_run = .false.
2365 :
2366 : contains
2367 :
2368 33 : subroutine check(i)
2369 : integer, intent(in) :: i
2370 : include 'formats'
2371 : !return
2372 :
2373 33 : current_num_iounits_in_use = number_iounits_allocated()
2374 33 : if (current_num_iounits_in_use > 3 .and. &
2375 : current_num_iounits_in_use > prev_num_iounits_in_use) then
2376 0 : write(*,2) 's% model_number', s% model_number
2377 0 : write(*,2) 'prev_num_iounits_in_use', prev_num_iounits_in_use
2378 0 : write(*,2) 'current_num_iounits_in_use', current_num_iounits_in_use
2379 0 : write(*,2) 'i', i
2380 0 : call mesa_error(__FILE__,__LINE__,'finish_step')
2381 : end if
2382 33 : prev_num_iounits_in_use = current_num_iounits_in_use
2383 11 : end subroutine check
2384 :
2385 : end function finish_step
2386 :
2387 :
2388 0 : subroutine set_age(id, age, ierr)
2389 : integer, intent(in) :: id
2390 : real(dp), intent(in) :: age
2391 : integer, intent(out) :: ierr
2392 : type (star_info), pointer :: s
2393 : include 'formats'
2394 0 : call get_star_ptr(id, s, ierr)
2395 0 : if (ierr /= 0) return
2396 0 : write(*,1) 'set_age', age
2397 0 : s% time = age*secyer
2398 0 : s% star_age = age
2399 : end subroutine set_age
2400 :
2401 : end module evolve
|
