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 timestep
21 :
22 : use star_private_def
23 : use utils_lib, only: is_bad
24 : use const_def, only: dp, ln10, secyer, msun, lsun, convective_mixing
25 : use chem_def
26 :
27 : implicit none
28 :
29 : private
30 : public :: timestep_controller
31 : public :: check_change
32 : public :: check_integer_limit
33 :
34 : logical, parameter :: dbg_timestep = .false.
35 : real(dp) :: max_dt
36 :
37 : contains
38 :
39 11 : integer function timestep_controller(s, max_timestep)
40 : ! if don't return keep_going, then set result_reason to say why.
41 : type (star_info), pointer :: s
42 : real(dp), intent(in) :: max_timestep
43 :
44 : include 'formats'
45 :
46 11 : timestep_controller = do_timestep_limits(s, s% dt)
47 11 : if (timestep_controller /= keep_going) s% result_reason = timestep_limits
48 :
49 11 : if (s% force_timestep > 0) then
50 0 : s% dt_next = s% force_timestep
51 0 : s% why_Tlim = Tlim_force_timestep
52 0 : return
53 : end if
54 :
55 : ! strictly enforce maximum timestep
56 11 : max_dt = max_timestep
57 11 : if (max_dt <= 0) max_dt = 1d99
58 11 : if (s% dt_next > max_dt) then
59 0 : s% dt_next = max_dt
60 0 : s% why_Tlim = Tlim_max_timestep
61 : end if
62 :
63 11 : if (s% timestep_hold > s% model_number .and. s% dt_next > s% dt) then
64 0 : s% dt_next = s% dt
65 0 : s% why_Tlim = Tlim_timestep_hold
66 0 : if (s% report_dt_hard_limit_retries .and. timestep_controller == keep_going) &
67 0 : write(*,3) 'timestep_hold > model_number, so no timestep increase', &
68 0 : s% timestep_hold, s% model_number
69 : end if
70 :
71 11 : if (is_bad_num(s% dt_next)) then
72 0 : write(*, *) 'timestep_controller: dt_next', s% dt_next
73 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'timestep_controller')
74 0 : timestep_controller = terminate
75 0 : s% termination_code = t_timestep_controller
76 0 : return
77 : end if
78 :
79 : if (dbg_timestep) then
80 : write(*,*) 'final result from timestep_controller for model_number', &
81 : timestep_controller, s% model_number
82 : write(*,1) 'lg dt/secyer', log10(s% dt/secyer)
83 : write(*,1) 'lg dt_next/secyer', log10(s% dt_next/secyer)
84 : write(*,1) 'dt_next/dt', s% dt_next/s% dt
85 : write(*,'(A)')
86 : write(*,'(A)')
87 : call mesa_error(__FILE__,__LINE__,'timestep_controller')
88 : end if
89 :
90 : end function timestep_controller
91 :
92 :
93 11 : integer function do_timestep_limits(s, dt)
94 : use rates_def, only: i_rate
95 : type (star_info), pointer :: s
96 : real(dp), intent(in) :: dt ! timestep just completed
97 :
98 : real(dp) :: dt_limit_ratio(numTlim), order, max_timestep_factor
99 : integer :: i_limit, nz, ierr
100 : logical :: skip_hard_limit
101 : integer :: num_mix_boundaries ! boundaries of regions with mixing_type /= no_mixing
102 11 : real(dp), pointer :: mix_bdy_q(:) ! (num_mix_boundaries)
103 : integer, pointer :: mix_bdy_loc(:) ! (num_mix_boundaries)
104 :
105 : include 'formats'
106 :
107 11 : if (s% never_skip_hard_limits) then
108 11 : skip_hard_limit = .false.
109 : else
110 : skip_hard_limit = (s% timestep_hold >= s% model_number) .or. &
111 : (s% relax_hard_limits_after_retry .and. &
112 0 : s% model_number_for_last_retry == s% model_number)
113 : end if
114 :
115 11 : nz = s% nz
116 :
117 : ! NOTE: when we get here, complete_model has called the report routine,
118 : ! so we can use information that it has calculated
119 :
120 11 : ierr = 0
121 :
122 11 : num_mix_boundaries = s% num_mix_boundaries
123 11 : mix_bdy_q => s% mix_bdy_q
124 11 : mix_bdy_loc => s% mix_bdy_loc
125 :
126 11 : dt_limit_ratio(:) = 0d0
127 :
128 : do_timestep_limits = check_varcontrol_limit( &
129 11 : s, dt_limit_ratio(Tlim_struc))
130 11 : if (return_now(Tlim_struc)) return
131 :
132 11 : if (.not. s% doing_first_model_of_run) then
133 :
134 10 : if (s% use_other_timestep_limit) then
135 : do_timestep_limits = s% other_timestep_limit( &
136 0 : s% id, skip_hard_limit, dt, dt_limit_ratio(Tlim_other_timestep_limit))
137 0 : if (return_now(Tlim_other_timestep_limit)) return
138 : end if
139 :
140 : do_timestep_limits = check_solver_iters_timestep_limit( &
141 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_solver))
142 10 : if (return_now(Tlim_solver)) return
143 :
144 : do_timestep_limits = check_burn_steps_limit( &
145 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_num_burn_steps))
146 10 : if (return_now(Tlim_num_burn_steps)) return
147 :
148 : do_timestep_limits = check_diffusion_steps_limit( &
149 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_num_diff_solver_steps))
150 10 : if (return_now(Tlim_num_diff_solver_steps)) return
151 :
152 : do_timestep_limits = check_diffusion_iters_limit( &
153 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_num_diff_solver_iters))
154 10 : if (return_now(Tlim_num_diff_solver_iters)) return
155 :
156 : do_timestep_limits = check_dX(s, skip_hard_limit, dt, &
157 : num_mix_boundaries, mix_bdy_loc, mix_bdy_q, &
158 10 : dt_limit_ratio(Tlim_dX), dt_limit_ratio(Tlim_dX_div_X))
159 10 : if (return_now(Tlim_dX_div_X)) return
160 :
161 : do_timestep_limits = check_dL_div_L( &
162 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dL_div_L))
163 10 : if (return_now(Tlim_dL_div_L)) return
164 :
165 : do_timestep_limits = check_dlgP_change( &
166 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_dlgP))
167 10 : if (return_now(Tlim_dlgP)) return
168 :
169 : do_timestep_limits = check_dlgRho_change( &
170 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_dlgRho))
171 10 : if (return_now(Tlim_dlgRho)) return
172 :
173 : do_timestep_limits = check_dlgT_change( &
174 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_dlgT))
175 10 : if (return_now(Tlim_dlgT)) return
176 :
177 : do_timestep_limits = check_dlgE_change( &
178 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_dlgE))
179 10 : if (return_now(Tlim_dlgE)) return
180 :
181 : do_timestep_limits = check_dlgR_change( &
182 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_dlgR))
183 10 : if (return_now(Tlim_dlgR)) return
184 :
185 : do_timestep_limits = check_lgL_nuc_cat_change( &
186 : s, num_mix_boundaries, mix_bdy_q, &
187 10 : skip_hard_limit, dt_limit_ratio(Tlim_dlgL_nuc_cat))
188 10 : if (return_now(Tlim_dlgL_nuc_cat)) return
189 :
190 : do_timestep_limits = check_lgL_H_change( &
191 20 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlgL_H))
192 10 : if (return_now(Tlim_dlgL_H)) return
193 :
194 : do_timestep_limits = check_lgL_He_change( &
195 20 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlgL_He))
196 10 : if (return_now(Tlim_dlgL_He)) return
197 :
198 : do_timestep_limits = check_lgL_z_change( &
199 20 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlgL_z))
200 10 : if (return_now(Tlim_dlgL_z)) return
201 :
202 : do_timestep_limits = check_lgL_power_photo_change( &
203 20 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlgL_power_photo))
204 10 : if (return_now(Tlim_dlgL_power_photo)) return
205 :
206 : do_timestep_limits = check_lgL_nuc_change( &
207 20 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlgL_nuc))
208 10 : if (return_now(Tlim_dlgL_nuc)) return
209 :
210 : do_timestep_limits = check_dlgTeff_change( &
211 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlgTeff))
212 10 : if (return_now(Tlim_dlgTeff)) return
213 :
214 : do_timestep_limits = check_dlgRho_cntr_change( &
215 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlgRho_cntr))
216 10 : if (return_now(Tlim_dlgRho_cntr)) return
217 :
218 : do_timestep_limits = check_dlgT_cntr_change( &
219 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlgT_cntr))
220 10 : if (return_now(Tlim_dlgT_cntr)) return
221 :
222 : do_timestep_limits = check_dlgT_max_change( &
223 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlgT_max))
224 10 : if (return_now(Tlim_dlgT_max)) return
225 :
226 : do_timestep_limits = check_dlgT_max_at_high_T_change( &
227 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlgT_max_at_high_T))
228 10 : if (return_now(Tlim_dlgT_max_at_high_T)) return
229 :
230 : do_timestep_limits = check_dlgP_cntr_change( &
231 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlgP_cntr))
232 10 : if (return_now(Tlim_dlgP_cntr)) return
233 :
234 : do_timestep_limits = check_dYe_highT_change( &
235 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_delta_Ye_highT))
236 10 : if (return_now(Tlim_delta_Ye_highT)) return
237 :
238 : do_timestep_limits = check_dlog_eps_nuc_change( &
239 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dlog_eps_nuc))
240 10 : if (return_now(Tlim_dlog_eps_nuc)) return
241 :
242 : do_timestep_limits = check_dX_div_X_cntr( &
243 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_dX_div_X_cntr))
244 10 : if (return_now(Tlim_dX_div_X_cntr)) return
245 :
246 : do_timestep_limits = check_lg_XH_cntr( &
247 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_lg_XH_cntr))
248 10 : if (return_now(Tlim_lg_XH_cntr)) return
249 :
250 : do_timestep_limits = check_lg_XHe_cntr( &
251 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_lg_XHe_cntr))
252 10 : if (return_now(Tlim_lg_XHe_cntr)) return
253 :
254 : do_timestep_limits = check_lg_XC_cntr( &
255 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_lg_XC_cntr))
256 10 : if (return_now(Tlim_lg_XC_cntr)) return
257 :
258 : do_timestep_limits = check_lg_XNe_cntr( &
259 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_lg_XNe_cntr))
260 10 : if (return_now(Tlim_lg_XNe_cntr)) return
261 :
262 : do_timestep_limits = check_lg_XO_cntr( &
263 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_lg_XO_cntr))
264 10 : if (return_now(Tlim_lg_XO_cntr)) return
265 :
266 : do_timestep_limits = check_lg_XSi_cntr( &
267 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_lg_XSi_cntr))
268 10 : if (return_now(Tlim_lg_XSi_cntr)) return
269 :
270 : do_timestep_limits = check_XH_cntr( &
271 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_XH_cntr))
272 10 : if (return_now(Tlim_XH_cntr)) return
273 :
274 : do_timestep_limits = check_XHe_cntr( &
275 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_XHe_cntr))
276 10 : if (return_now(Tlim_XHe_cntr)) return
277 :
278 : do_timestep_limits = check_XC_cntr( &
279 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_XC_cntr))
280 10 : if (return_now(Tlim_XC_cntr)) return
281 :
282 : do_timestep_limits = check_XNe_cntr( &
283 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_XNe_cntr))
284 10 : if (return_now(Tlim_XNe_cntr)) return
285 :
286 : do_timestep_limits = check_XO_cntr( &
287 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_XO_cntr))
288 10 : if (return_now(Tlim_XO_cntr)) return
289 :
290 : do_timestep_limits = check_XSi_cntr( &
291 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_XSi_cntr))
292 10 : if (return_now(Tlim_XSi_cntr)) return
293 :
294 : do_timestep_limits = check_delta_mstar( &
295 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dmstar))
296 10 : if (return_now(Tlim_dmstar)) return
297 :
298 : do_timestep_limits = check_delta_mdot( &
299 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_del_mdot))
300 10 : if (return_now(Tlim_del_mdot)) return
301 :
302 : do_timestep_limits = check_adjust_J_q( &
303 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_adjust_J_q))
304 10 : if (return_now(Tlim_adjust_J_q)) return
305 :
306 : do_timestep_limits = check_delta_lgL( &
307 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_lgL))
308 10 : if (return_now(Tlim_lgL)) return
309 :
310 : do_timestep_limits = check_dt_div_dt_cell_collapse( &
311 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dt_div_dt_cell_collapse))
312 10 : if (return_now(Tlim_dt_div_dt_cell_collapse)) return
313 :
314 : do_timestep_limits = check_dt_div_min_dr_div_cs( &
315 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dt_div_min_dr_div_cs))
316 10 : if (return_now(Tlim_dt_div_min_dr_div_cs)) return
317 :
318 : do_timestep_limits = check_delta_HR( &
319 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_delta_HR))
320 10 : if (return_now(Tlim_delta_HR)) return
321 :
322 : do_timestep_limits = check_rel_error_in_energy( &
323 10 : s, skip_hard_limit, dt_limit_ratio(Tlim_error_in_energy_conservation))
324 10 : if (return_now(Tlim_error_in_energy_conservation)) return
325 :
326 : do_timestep_limits = check_dX_nuc_drop( &
327 10 : s, skip_hard_limit, dt, dt_limit_ratio(Tlim_dX_nuc_drop))
328 10 : if (return_now(Tlim_dX_nuc_drop)) return
329 :
330 : end if
331 :
332 649 : i_limit = maxloc(dt_limit_ratio(1:numTlim), dim=1)
333 :
334 11 : order = 1
335 11 : call filter_dt_next(s, order, dt_limit_ratio(i_limit)) ! sets s% dt_next
336 :
337 11 : if (s% log_max_temperature > s% min_logT_for_max_timestep_factor_at_high_T) then
338 0 : max_timestep_factor = s% max_timestep_factor_at_high_T
339 : else
340 11 : max_timestep_factor = s% max_timestep_factor
341 : end if
342 :
343 11 : if (max_timestep_factor > 0 .and. s% dt_next > max_timestep_factor*s% dt) then
344 11 : s% dt_next = max_timestep_factor*s% dt
345 11 : if (s% report_solver_dt_info) then
346 0 : write(*,2) 's% dt', s% model_number, s% dt
347 0 : write(*,2) 'max_timestep_factor', s% model_number, max_timestep_factor
348 0 : write(*,2) 's% dt_next', s% model_number, s% dt_next
349 : end if
350 11 : if (s% dt_next == 0d0) call mesa_error(__FILE__,__LINE__,'filter_dt_next')
351 11 : if (i_limit == Tlim_struc) i_limit = Tlim_max_timestep_factor
352 : end if
353 :
354 11 : if (s% min_timestep_factor > 0 .and. s% dt_next < s% min_timestep_factor*s% dt) then
355 0 : s% dt_next = s% min_timestep_factor*s% dt
356 0 : if (s% report_solver_dt_info) then
357 0 : write(*,2) 's% dt', s% model_number, s% dt
358 0 : write(*,2) 'min_timestep_factor', s% model_number, s% min_timestep_factor
359 0 : write(*,2) 's% dt_next', s% model_number, s% dt_next
360 : end if
361 0 : if (s% dt_next == 0d0) call mesa_error(__FILE__,__LINE__,'filter_dt_next')
362 0 : if (i_limit == Tlim_struc) i_limit = Tlim_min_timestep_factor
363 : end if
364 :
365 11 : s% why_Tlim = i_limit
366 11 : if (i_limit > 0) s% dt_why_count(i_limit) = s% dt_why_count(i_limit) + 1
367 :
368 : contains
369 :
370 481 : logical function return_now(i_limit)
371 : integer, intent(in) :: i_limit
372 481 : if (do_timestep_limits == keep_going) then
373 481 : return_now = .false.
374 : return
375 : end if
376 0 : return_now = .true.
377 0 : s% why_Tlim = i_limit
378 0 : if (i_limit > 0) s% dt_why_retry_count(i_limit) = s% dt_why_retry_count(i_limit) + 1
379 : end function return_now
380 :
381 : end function do_timestep_limits
382 :
383 :
384 10 : integer function check_integer_limit( &
385 : s, limit, hard_limit, value, msg, skip_hard_limit, dt, dt_limit_ratio)
386 : type (star_info), pointer :: s
387 : integer, intent(in) :: limit, hard_limit, value
388 : character (len=*), intent(in) :: msg
389 : logical, intent(in) :: skip_hard_limit
390 : real(dp), intent(in) :: dt
391 : real(dp), intent(inout) :: dt_limit_ratio
392 : include 'formats'
393 10 : if (value > hard_limit .and. hard_limit > 0 .and. (.not. skip_hard_limit)) then
394 0 : check_integer_limit = retry
395 0 : s% retry_message = trim(msg) // ' hard limit'
396 0 : if (s% report_dt_hard_limit_retries) then
397 0 : write(*,*) trim(msg) // ' hard limit', hard_limit, value
398 0 : write(*,3) trim(msg), s% model_number, value
399 0 : write(*,3) trim(msg) // ' hard limit', s% model_number, hard_limit
400 : end if
401 0 : return
402 : end if
403 10 : check_integer_limit = keep_going
404 10 : if (value <= 0 .or. limit <= 0) return
405 10 : dt_limit_ratio = dble(value)/dble(limit) - 0.05d0
406 : ! subtract a bit so that allow dt to grow even if value == limit
407 10 : if (dt_limit_ratio <= 1d0) dt_limit_ratio = 0
408 : end function check_integer_limit
409 :
410 :
411 10 : integer function check_burn_steps_limit(s, skip_hard_limit, dt, dt_limit_ratio)
412 : type (star_info), pointer :: s
413 : logical, intent(in) :: skip_hard_limit
414 : real(dp), intent(in) :: dt
415 : real(dp), intent(inout) :: dt_limit_ratio
416 : integer :: max_steps
417 10 : check_burn_steps_limit = keep_going
418 10633 : if (.not. s% op_split_burn .or. maxval(s% T_start(1:s%nz)) < s% op_split_burn_min_T) return
419 :
420 0 : max_steps = maxval(s% burn_num_iters(1:s% nz),mask=s% T(1:s%nz)>s% op_split_burn_min_T)
421 : check_burn_steps_limit = check_integer_limit( &
422 : s, s% burn_steps_limit, s% burn_steps_hard_limit, max_steps, &
423 0 : 'num_burn_solver_steps', skip_hard_limit, dt, dt_limit_ratio)
424 0 : end function check_burn_steps_limit
425 :
426 :
427 10 : integer function check_solver_iters_timestep_limit(s, skip_hard_limit, dt, dt_limit_ratio)
428 : type (star_info), pointer :: s
429 : logical, intent(in) :: skip_hard_limit
430 : real(dp), intent(in) :: dt
431 : real(dp), intent(inout) :: dt_limit_ratio
432 : integer :: iters, limit, hard_limit
433 : include 'formats'
434 10 : check_solver_iters_timestep_limit = keep_going
435 10 : iters = s% num_solver_iterations
436 10 : limit = s% solver_iters_timestep_limit
437 10 : hard_limit = -1
438 10 : if (s% using_gold_tolerances) then
439 10 : limit = s% gold_solver_iters_timestep_limit
440 10 : if ((s% use_gold2_tolerances .and. s% steps_before_use_gold2_tolerances < 0) .or. &
441 : (s% steps_before_use_gold2_tolerances >= 0 .and. &
442 : s% model_number > s% steps_before_use_gold2_tolerances + max(0,s% init_model_number))) then
443 0 : limit = s% gold2_solver_iters_timestep_limit
444 : end if
445 : end if
446 10 : if (s% used_extra_iter_in_solver_for_accretion) iters = iters - 1
447 : check_solver_iters_timestep_limit = check_integer_limit( &
448 : s, limit, hard_limit, &
449 10 : iters, 'num_solver_iterations', skip_hard_limit, dt, dt_limit_ratio)
450 10 : end function check_solver_iters_timestep_limit
451 :
452 :
453 10 : integer function check_diffusion_steps_limit(s, skip_hard_limit, dt, dt_limit_ratio)
454 : type (star_info), pointer :: s
455 : logical, intent(in) :: skip_hard_limit
456 : real(dp), intent(in) :: dt
457 : real(dp), intent(inout) :: dt_limit_ratio
458 10 : check_diffusion_steps_limit = keep_going
459 10 : if (.not. s% do_element_diffusion) return
460 : check_diffusion_steps_limit = check_integer_limit( &
461 : s, s% diffusion_steps_limit, s% diffusion_steps_hard_limit, &
462 : s% num_diffusion_solver_steps, &
463 0 : 'num_diffusion_solver_steps', skip_hard_limit, dt, dt_limit_ratio)
464 0 : end function check_diffusion_steps_limit
465 :
466 :
467 10 : integer function check_diffusion_iters_limit(s, skip_hard_limit, dt, dt_limit_ratio)
468 : type (star_info), pointer :: s
469 : logical, intent(in) :: skip_hard_limit
470 : real(dp), intent(in) :: dt
471 : real(dp), intent(inout) :: dt_limit_ratio
472 10 : check_diffusion_iters_limit = keep_going
473 10 : if (.not. s% do_element_diffusion) return
474 : check_diffusion_iters_limit = check_integer_limit( &
475 : s, s% diffusion_iters_limit, s% diffusion_iters_hard_limit, &
476 : s% num_diffusion_solver_iters, &
477 0 : 'num_diffusion_solver_iters', skip_hard_limit, dt, dt_limit_ratio)
478 0 : end function check_diffusion_iters_limit
479 :
480 :
481 10 : integer function check_dX(s, skip_hard_limit, dt, &
482 : n_mix_bdy, mix_bdy_loc, mix_bdy_q, &
483 : dX_dt_limit_ratio, dX_div_X_dt_limit_ratio)
484 : use num_lib, only: binary_search
485 : logical, intent(in) :: skip_hard_limit
486 : integer, intent(in) :: n_mix_bdy, mix_bdy_loc(:)
487 : real(dp), intent(in) :: dt
488 : real(dp), intent(in), pointer :: mix_bdy_q(:)
489 : real(dp), intent(inout) :: dX_dt_limit_ratio, dX_div_X_dt_limit_ratio
490 :
491 : type (star_info), pointer :: s
492 : real(dp) :: X, X_old, delta_X, delta_X_div_X, max_dX, max_dX_div_X, &
493 : bdy_dist_dm, max_dX_bdy_dist_dm, max_dX_div_X_bdy_dist_dm, cz_dist_limit, &
494 : D_mix_cutoff, ratio_tmp_dX, ratio_tmp_dX_div_X
495 : integer :: i, j, k, cid, bdy, max_dX_j, max_dX_k, max_dX_div_X_j, max_dX_div_X_k
496 : real(dp), dimension(max_dX_limit_ctrls) :: dX_limit, dX_hard_limit, &
497 : dX_div_X_limit, dX_div_X_hard_limit
498 : character (len=strlen) :: sp
499 :
500 : include 'formats'
501 :
502 10 : check_dX = keep_going
503 10 : dX_dt_limit_ratio = 0d0
504 10 : dX_div_X_dt_limit_ratio = 0d0
505 :
506 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
507 :
508 1010 : do i=1, max_dx_limit_ctrls ! go over all potential species + XYZ
509 : if (s% dX_limit(i) >= 1 .and. & ! as soon as all of these are >= 1
510 : s% dX_hard_limit(i) >= 1 .and. & ! we'd have nothing to do
511 1000 : s% dX_div_X_limit(i) >= 1 .and. &
512 : s% dX_div_X_hard_limit(i) >= 1) then
513 : cycle ! go to next
514 : end if
515 :
516 2020 : dX_limit = s% dX_limit(i) * s% time_delta_coeff
517 2020 : dX_hard_limit = s% dX_hard_limit(i) * s% time_delta_coeff
518 :
519 20 : if (s% log_max_temperature > s% dX_div_X_at_high_T_limit_lgT_min(i)) then
520 0 : dX_div_X_limit = s% dX_div_X_at_high_T_limit(i)
521 0 : dX_div_X_hard_limit = s% dX_div_X_at_high_T_hard_limit(i)
522 : else
523 2020 : dX_div_X_limit = s% dX_div_X_limit(i)
524 2020 : dX_div_X_hard_limit = s% dX_div_X_hard_limit(i)
525 : end if
526 :
527 2020 : dX_div_X_limit = dX_div_X_limit * s% time_delta_coeff
528 2020 : dX_div_X_hard_limit = dX_div_X_hard_limit * s% time_delta_coeff
529 :
530 20 : max_dX = -1; max_dX_j = -1; max_dX_k = -1
531 20 : max_dX_div_X = -1; max_dX_div_X_j = -1; max_dX_div_X_k = -1
532 20 : bdy = 0
533 20 : max_dX_bdy_dist_dm = 0
534 20 : max_dX_div_X_bdy_dist_dm = 0
535 20 : cz_dist_limit = s% dX_mix_dist_limit*Msun
536 :
537 20 : if (s% set_min_D_mix .and. s% ye(s% nz) >= s% min_center_Ye_for_min_D_mix) then
538 0 : D_mix_cutoff = s% min_D_mix
539 : else
540 : D_mix_cutoff = 0
541 : end if
542 :
543 20 : sp = s% dX_limit_species(i)
544 :
545 21266 : do k = 1, s% nz
546 :
547 21246 : if (s% D_mix(k) > D_mix_cutoff) then
548 : cycle
549 : end if
550 18378 : if (k < s% nz) then
551 18378 : if (s% D_mix(k+1) > D_mix_cutoff) then
552 : cycle
553 : end if
554 : end if
555 :
556 : ! find the nearest mixing boundary
557 18338 : bdy = binary_search(n_mix_bdy, mix_bdy_q, bdy, s% q(k))
558 : ! don't check cells near a mixing boundary
559 18338 : if (bdy > 0 .and. bdy < n_mix_bdy) then
560 16398 : bdy_dist_dm = s% xmstar*abs(s% q(k) - mix_bdy_q(bdy))
561 16398 : if (bdy_dist_dm < cz_dist_limit) cycle
562 : else
563 : bdy_dist_dm = 0
564 : end if
565 :
566 165062 : do j = 1, s% species
567 :
568 146704 : cid = s% chem_id(j)
569 146704 : if (sp == 'X') then ! any hydrogen
570 0 : if (cid /= ih1 .or. cid /= ih2 .or. cid /= ih3) cycle
571 146704 : else if (sp == 'Y') then ! any helium
572 0 : if (cid /= ihe4 .or. cid /= ihe3) cycle
573 146704 : else if (sp == 'Z') then ! any metal
574 0 : if (chem_isos% Z(cid) <= 2) cycle
575 : else ! single isotope
576 146704 : if (trim(chem_isos% name(s% chem_id(j))) /= trim(sp)) cycle
577 : end if
578 :
579 18338 : X = s% xa(j,k)
580 18338 : X_old = s% xa_old(j,k)
581 18338 : delta_X = X_old - X ! decrease in abundance
582 :
583 18338 : if ((.not. s% dX_decreases_only(j)) .and. delta_X < 0) delta_X = -delta_X
584 :
585 18338 : if (X >= s% dX_limit_min_X(i)) then ! any check for dX_limit_* < 1 is useless since X <= 1 anyway
586 0 : if ((.not. skip_hard_limit) .and. delta_X > dX_hard_limit(i)) then
587 0 : check_dX = retry
588 0 : s% why_Tlim = Tlim_dX
589 0 : s% Tlim_dX_species = j
590 0 : s% Tlim_dX_cell = k
591 0 : s% retry_message = 'dX ' // trim(chem_isos% name(s% chem_id(j))) // ' hard limit'
592 0 : s% retry_message_k = k
593 0 : if (s% report_dt_hard_limit_retries) then
594 0 : write(*,2) 'old xa', s% model_number, X_old
595 0 : write(*,2) 'new xa', s% model_number, X
596 0 : write(*,2) 'delta xa', s% model_number, delta_X
597 0 : write(*,2) 'hard limit delta xa', s% model_number, dX_hard_limit
598 : end if
599 0 : return
600 : end if
601 0 : if (delta_X > max_dX) then
602 18338 : max_dX = delta_X
603 18338 : max_dX_j = j
604 18338 : max_dX_k = k
605 18338 : max_dX_bdy_dist_dm = bdy_dist_dm
606 : end if
607 : end if
608 39584 : if (X >= s% dX_div_X_limit_min_X(i)) then
609 18338 : delta_X_div_X = delta_X/X
610 18338 : if ((.not. skip_hard_limit) .and. delta_X_div_X > dX_div_X_hard_limit(i)) then
611 0 : check_dX = retry
612 0 : s% why_Tlim = Tlim_dX_div_X
613 0 : s% Tlim_dX_div_X_species = j
614 0 : s% Tlim_dX_div_X_cell = k
615 0 : s% retry_message = 'dX_div_X ' // trim(chem_isos% name(s% chem_id(j))) // ' hard limit'
616 0 : s% retry_message_k = k
617 0 : if (s% report_dt_hard_limit_retries) then
618 : write(*, '(a30, i5, 99(/,a30,e20.10))') &
619 0 : 'delta_X_div_X ' // trim(chem_isos% name(s% chem_id(j))), &
620 0 : k, 'delta_X_div_X', delta_X_div_X, &
621 0 : 'dX_div_X_hard_limit', dX_div_X_hard_limit
622 0 : write(*,2) 'old xa', s% model_number, X_old
623 0 : write(*,2) 'new xa', s% model_number, X
624 0 : write(*,2) 'delta_X_div_X', s% model_number, delta_X_div_X
625 0 : write(*,2) 'dX_div_X_hard_limit', s% model_number, dX_div_X_hard_limit
626 : end if
627 0 : return
628 : end if
629 18338 : if (delta_X_div_X > max_dX_div_X) then
630 3787 : max_dX_div_X = delta_X_div_X
631 3787 : max_dX_div_X_j = j
632 3787 : max_dX_div_X_k = k
633 3787 : max_dX_div_X_bdy_dist_dm = bdy_dist_dm
634 : end if
635 : end if
636 : end do
637 : end do
638 :
639 20 : if (dX_limit(i) > 0) then
640 20 : ratio_tmp_dX = max_dX/dX_limit(i)
641 20 : if (ratio_tmp_dX > dX_dt_limit_ratio) then
642 0 : dX_dt_limit_ratio = ratio_tmp_dX
643 0 : if (dX_dt_limit_ratio <= 1d0) then
644 0 : dX_dt_limit_ratio = 0
645 : else
646 0 : s% Tlim_dX_species = max_dX_j
647 0 : s% Tlim_dX_cell = max_dX_k
648 : ! write(*, '(a30, i5, 99e20.10)') &
649 : ! ' limit dt because of large dX '// &
650 : ! trim(chem_isos% name(s% chem_id(max_dX_j))) // &
651 : ! ' k, max, lim, m ', &
652 : ! max_dX_k, max_dX, dX_limit(i), &
653 : ! max_dX_bdy_dist_dm/Msun
654 : end if
655 : end if
656 : end if
657 :
658 30 : if (dX_div_X_limit(i) > 0) then
659 20 : ratio_tmp_dX_div_X = max_dX_div_X/dX_div_X_limit(i)
660 20 : if (ratio_tmp_dX_div_X > dX_div_X_dt_limit_ratio) then ! pick out largest culprit only!
661 10 : dX_div_X_dt_limit_ratio = ratio_tmp_dX_div_X
662 10 : if (dX_div_X_dt_limit_ratio <= 1d0) then
663 10 : dX_div_X_dt_limit_ratio = 0
664 : else
665 0 : s% Tlim_dX_div_X_species = max_dX_div_X_j
666 0 : s% Tlim_dX_div_X_cell = max_dX_div_X_k
667 : ! write(*, '(a35, i5, 99e20.10)') & ! shouldn't be written as is isn't guarantueed
668 : ! this control will trigger timestep reduction
669 : ! ' limit dt because of large dX_div_X ' // &
670 : ! trim(chem_isos% name(s% chem_id(max_dX_div_X_j))) // &
671 : ! ' k, max, lim, m ', &
672 : ! max_dX_div_X_k, max_dX_div_X, dX_div_X_limit(i), &
673 : ! max_dX_div_X_bdy_dist_dm/Msun
674 : end if
675 : end if
676 : end if
677 : end do
678 : end function check_dX
679 :
680 :
681 10 : integer function check_dL_div_L(s, skip_hard_limit, dt, dL_div_L_dt_ratio)
682 : type (star_info), pointer :: s
683 : logical, intent(in) :: skip_hard_limit
684 : real(dp), intent(in) :: dt
685 : real(dp), intent(inout) :: dL_div_L_dt_ratio
686 :
687 : real(dp) :: L, abs_dL, abs_dL_div_L, max_dL_div_L
688 : integer :: k, max_dL_div_L_k
689 : real(dp) :: dL_div_L_limit_min_L, dL_div_L_limit, dL_div_L_hard_limit
690 :
691 : include 'formats'
692 :
693 10 : check_dL_div_L = keep_going
694 :
695 10 : dL_div_L_limit_min_L = Lsun*s% dL_div_L_limit_min_L
696 10 : dL_div_L_limit = s% dL_div_L_limit*s% time_delta_coeff
697 10 : dL_div_L_hard_limit = s% dL_div_L_hard_limit*s% time_delta_coeff
698 :
699 10 : if (dL_div_L_limit_min_L <= 0) return
700 10 : if (dL_div_L_limit <= 0 .and. dL_div_L_hard_limit <= 0) return
701 :
702 0 : max_dL_div_L = -1
703 0 : max_dL_div_L_k = -1
704 0 : abs_dL_div_L = 0; L=0 ! to quiet gfortran
705 :
706 0 : do k = 1, s% nz
707 0 : L = s% L(k)
708 0 : abs_dL = abs(L - s% L_start(k))
709 0 : if (L >= dL_div_L_limit_min_L) then
710 0 : abs_dL_div_L = abs_dL/L
711 : if (dL_div_L_hard_limit > 0 .and. (.not. skip_hard_limit) &
712 0 : .and. abs_dL_div_L > dL_div_L_hard_limit) then
713 0 : check_dL_div_L= retry
714 0 : s% retry_message = 'dL_div_L hard limit'
715 0 : s% retry_message_k = k
716 0 : if (s% report_dt_hard_limit_retries) then
717 0 : write(*,2) 'L', L
718 0 : write(*,2) 'L_start', s% L_start(k)
719 0 : write(*,2) 'abs_dL_div_L', abs_dL_div_L
720 0 : write(*,2) 'dL_div_L_hard_limit', dL_div_L_hard_limit
721 : end if
722 0 : return
723 : end if
724 0 : if (abs_dL_div_L > max_dL_div_L) then
725 0 : max_dL_div_L = abs_dL_div_L
726 0 : max_dL_div_L_k = k
727 : end if
728 : end if
729 : end do
730 :
731 0 : if (dL_div_L_limit <= 0) return
732 0 : dL_div_L_dt_ratio = max_dL_div_L/dL_div_L_limit
733 :
734 0 : end function check_dL_div_L
735 :
736 :
737 160 : integer function check_change( &
738 : s, delta_value, lim_in, hard_lim_in, max_k, msg, &
739 : skip_hard_limit, dt_limit_ratio, relative_excess)
740 : type (star_info), pointer :: s
741 : real(dp), intent(in) :: delta_value, lim_in, hard_lim_in
742 : integer, intent(in) :: max_k
743 : character (len=*), intent(in) :: msg
744 : logical, intent(in) :: skip_hard_limit
745 : real(dp), intent(inout) :: dt_limit_ratio
746 : real(dp), intent(out) :: relative_excess
747 : real(dp) :: abs_change, lim, hard_lim
748 : include 'formats'
749 160 : if (is_bad(delta_value)) then
750 0 : write(*,1) trim(msg) // ' delta_value', delta_value
751 0 : call mesa_error(__FILE__,__LINE__,'check_change')
752 : end if
753 160 : check_change = keep_going
754 160 : abs_change = abs(delta_value)
755 160 : lim = lim_in*s% time_delta_coeff
756 160 : hard_lim = hard_lim_in*s% time_delta_coeff
757 160 : if (hard_lim > 0 .and. abs_change > hard_lim .and. (.not. skip_hard_limit)) then
758 0 : s% retry_message = trim(msg) // ' hard limit'
759 0 : s% retry_message_k = max_k
760 0 : check_change = retry
761 0 : return
762 : end if
763 160 : if (lim <= 0) return
764 120 : relative_excess = (abs_change - lim) / lim
765 120 : dt_limit_ratio = abs_change/lim ! 1d0/(s% timestep_dt_factor**relative_excess)
766 120 : if (is_bad(dt_limit_ratio)) then
767 0 : write(*,1) trim(msg) // ' dt_limit_ratio', dt_limit_ratio, abs_change, lim
768 0 : call mesa_error(__FILE__,__LINE__,'check_change')
769 : end if
770 120 : if (dt_limit_ratio <= 1d0) dt_limit_ratio = 0
771 : end function check_change
772 :
773 :
774 10 : subroutine get_dlgP_info(s, i, max_dlnP)
775 : type (star_info), pointer :: s
776 : integer, intent(out) :: i
777 : real(dp), intent(out) :: max_dlnP
778 : real(dp) :: lim, dlnP
779 : integer :: k
780 : include 'formats'
781 10 : lim = ln10*s% delta_lgP_limit_min_lgP
782 10 : i = 0
783 10 : max_dlnP = 0
784 10633 : do k=1,s% nz
785 10623 : if (s% lnPeos(k) < lim) cycle
786 0 : dlnP = abs(s% lnPeos(k) - s% lnPeos_start(k))
787 10 : if (dlnP > max_dlnP) then
788 0 : max_dlnP = dlnP
789 0 : i = k
790 : end if
791 : end do
792 10 : end subroutine get_dlgP_info
793 :
794 :
795 10 : integer function check_dlgP_change(s, skip_hard_limit, dt_limit_ratio)
796 : type (star_info), pointer :: s
797 : logical, intent(in) :: skip_hard_limit
798 : real(dp), intent(inout) :: dt_limit_ratio
799 : real(dp) :: relative_excess, max_dlnP
800 : integer :: i
801 : include 'formats'
802 10 : check_dlgP_change = keep_going
803 10 : call get_dlgP_info(s, i, max_dlnP)
804 10 : if (i == 0) return
805 : check_dlgP_change = check_change(s, max_dlnP/ln10, &
806 : s% delta_lgP_limit, s% delta_lgP_hard_limit, &
807 0 : i, 'delta_lgP', skip_hard_limit, dt_limit_ratio, relative_excess)
808 0 : if (check_dlgP_change /= keep_going .and. s% report_dt_hard_limit_retries) then
809 0 : write(*,3) 'lgP', i, s% lnPeos(i)/ln10
810 0 : write(*,3) 'lgP_old', i, s% lnPeos_start(i)/ln10
811 0 : write(*,3) 'dlgP', i, (s% lnPeos(i) - s% lnPeos_start(i))/ln10
812 0 : write(*,3) 'hard_limit', i, s% delta_lgP_hard_limit
813 : end if
814 : end function check_dlgP_change
815 :
816 :
817 10 : subroutine get_dlgRho_info(s, i, max_dlnRho)
818 : type (star_info), pointer :: s
819 : integer, intent(out) :: i
820 : real(dp), intent(out) :: max_dlnRho
821 : real(dp) :: lim, dlnRho, max_abs_dlnRho
822 : integer :: k
823 : include 'formats'
824 10 : lim = ln10*s% delta_lgRho_limit_min_lgRho
825 10 : i = 0
826 10 : max_abs_dlnRho = 0
827 10633 : do k=1,s% nz
828 10623 : if (s% lnd(k) < lim) cycle
829 0 : dlnRho = s% lnd(k) - s% lnd_start(k)
830 10 : if (abs(dlnRho) > max_abs_dlnRho) then
831 0 : max_dlnRho = dlnRho
832 0 : max_abs_dlnRho = abs(dlnRho)
833 0 : i = k
834 : end if
835 : end do
836 10 : end subroutine get_dlgRho_info
837 :
838 :
839 10 : integer function check_dlgRho_change(s, skip_hard_limit, dt_limit_ratio)
840 : ! check max change in log10(density)
841 : type (star_info), pointer :: s
842 : logical, intent(in) :: skip_hard_limit
843 : real(dp), intent(inout) :: dt_limit_ratio
844 : real(dp) :: relative_excess, max_dlnd
845 : integer :: i
846 : include 'formats'
847 10 : check_dlgRho_change = keep_going
848 10 : call get_dlgRho_info(s, i, max_dlnd)
849 10 : if (i == 0) return
850 : check_dlgRho_change = check_change(s, abs(max_dlnd)/ln10, &
851 : s% delta_lgRho_limit, s% delta_lgRho_hard_limit, &
852 0 : i, 'delta_lgRho', skip_hard_limit, dt_limit_ratio, relative_excess)
853 0 : if (check_dlgRho_change /= keep_going .and. s% report_dt_hard_limit_retries) then
854 0 : write(*,3) 'lgRho', i, s% lnd(i)/ln10
855 0 : write(*,3) 'lgRho_old', i, s% lnd_start(i)/ln10
856 0 : write(*,3) 'dlgRho', i, (s% lnd(i) - s% lnd_start(i))/ln10
857 0 : write(*,3) 'hard_limit', i, s% delta_lgRho_hard_limit
858 : end if
859 : end function check_dlgRho_change
860 :
861 :
862 10 : subroutine get_dlgT_info(s, i, max_dlnT)
863 : type (star_info), pointer :: s
864 : integer, intent(out) :: i
865 : real(dp), intent(out) :: max_dlnT
866 : real(dp) :: lim, dlnT, abs_max_dlnT
867 : integer :: k
868 : include 'formats'
869 10 : lim = ln10*s% delta_lgT_limit_min_lgT
870 10 : i = 0
871 10 : abs_max_dlnT = 0
872 10633 : do k=1,s% nz
873 10623 : if (s% lnT(k) < lim) cycle
874 0 : dlnT = s% lnT(k) - s% lnT_start(k)
875 10 : if (abs(dlnT) > abs_max_dlnT) then
876 0 : max_dlnT = dlnT
877 0 : abs_max_dlnT = abs(dlnT)
878 0 : i = k
879 : end if
880 : end do
881 10 : end subroutine get_dlgT_info
882 :
883 :
884 10 : integer function check_dlgT_change(s, skip_hard_limit, dt_limit_ratio)
885 : ! check max change in log10(temperature)
886 : type (star_info), pointer :: s
887 : logical, intent(in) :: skip_hard_limit
888 : real(dp), intent(inout) :: dt_limit_ratio
889 : real(dp) :: relative_excess, max_dlnT
890 : integer :: i
891 : include 'formats'
892 10 : check_dlgT_change = keep_going
893 10 : call get_dlgT_info(s, i, max_dlnT)
894 10 : if (i == 0) return
895 : check_dlgT_change = check_change(s, abs(max_dlnT)/ln10, &
896 : s% delta_lgT_limit, s% delta_lgT_hard_limit, &
897 0 : i, 'delta_lgT', skip_hard_limit, dt_limit_ratio, relative_excess)
898 0 : if (check_dlgT_change /= keep_going .and. s% report_dt_hard_limit_retries) then
899 0 : write(*,3) 'lgT', i, s% lnT(i)/ln10
900 0 : write(*,3) 'lgT_old', i, s% lnT_start(i)/ln10
901 0 : write(*,3) 'dlgT', i, (s% lnT(i) - s% lnT_start(i))/ln10
902 0 : write(*,3) 'hard_limit', i, s% delta_lgT_hard_limit
903 : end if
904 : end function check_dlgT_change
905 :
906 :
907 10 : subroutine get_dlgE_info(s, i, max_dlnE)
908 : type (star_info), pointer :: s
909 : integer, intent(out) :: i
910 : real(dp), intent(out) :: max_dlnE
911 : real(dp) :: lim, dlnE
912 : integer :: k
913 : include 'formats'
914 10 : lim = ln10*s% delta_lgE_limit_min_lgE
915 10 : i = 0
916 10 : max_dlnE = 0
917 10633 : do k=1,s% nz
918 10623 : if (s% lnE(k) < lim) cycle
919 0 : dlnE = abs(s% energy(k) - s% energy_start(k))/s% energy(k)
920 10 : if (dlnE > max_dlnE) then
921 0 : max_dlnE = dlnE
922 0 : i = k
923 : end if
924 : end do
925 10 : end subroutine get_dlgE_info
926 :
927 :
928 10 : integer function check_dlgE_change(s, skip_hard_limit, dt_limit_ratio)
929 : ! check max change in log10(internal energy)
930 : type (star_info), pointer :: s
931 : logical, intent(in) :: skip_hard_limit
932 : real(dp), intent(inout) :: dt_limit_ratio
933 : real(dp) :: relative_excess, max_dlnE
934 : integer :: i
935 : include 'formats'
936 10 : check_dlgE_change = keep_going
937 10 : call get_dlgE_info(s, i, max_dlnE)
938 10 : if (i == 0) return
939 : check_dlgE_change = check_change(s, max_dlnE/ln10, &
940 : s% delta_lgE_limit, s% delta_lgE_hard_limit, &
941 0 : i, 'delta_lgE', skip_hard_limit, dt_limit_ratio, relative_excess)
942 0 : if (check_dlgE_change /= keep_going .and. s% report_dt_hard_limit_retries) then
943 0 : write(*,3) 'lgE', i, safe_log10(s% energy(i))
944 0 : write(*,3) 'lgE_old', i, safe_log10(s% energy_start(i))
945 0 : write(*,3) 'dlgE', i, (s% energy(i) - s% energy_start(i))/s% energy(i)/ln10
946 0 : write(*,3) 'hard_limit', i, s% delta_lgE_hard_limit
947 : end if
948 : end function check_dlgE_change
949 :
950 :
951 10 : subroutine get_dlgR_info(s, i, max_dlnR)
952 : type (star_info), pointer :: s
953 : integer, intent(out) :: i
954 : real(dp), intent(out) :: max_dlnR
955 : real(dp) :: lim, dlnR
956 : integer :: k
957 : include 'formats'
958 10 : lim = ln10*s% delta_lgR_limit_min_lgR
959 10 : i = 0
960 10 : max_dlnR = 0
961 10633 : do k=1,s% nz
962 10623 : if (s% lnR(k) < lim) cycle
963 0 : dlnR = abs(s% lnR(k) - s% lnR_start(k))
964 10 : if (dlnR > max_dlnR) then
965 0 : max_dlnR = dlnR
966 0 : i = k
967 : end if
968 : end do
969 10 : end subroutine get_dlgR_info
970 :
971 :
972 10 : integer function check_dlgR_change(s, skip_hard_limit, dt_limit_ratio)
973 : type (star_info), pointer :: s
974 : logical, intent(in) :: skip_hard_limit
975 : real(dp), intent(inout) :: dt_limit_ratio
976 : real(dp) :: relative_excess, max_dlnR
977 : integer :: i
978 : include 'formats'
979 10 : check_dlgR_change = keep_going
980 10 : call get_dlgR_info(s, i, max_dlnR)
981 10 : if (i == 0) return
982 : check_dlgR_change = check_change(s, max_dlnR/ln10, &
983 : s% delta_lgR_limit, s% delta_lgR_hard_limit, &
984 0 : i, 'delta_lgR', skip_hard_limit, dt_limit_ratio, relative_excess)
985 0 : if (check_dlgR_change /= keep_going .and. s% report_dt_hard_limit_retries) then
986 0 : write(*,3) 'lgR', i, s% lnR(i)/ln10
987 0 : write(*,3) 'lgR_old', i, s% lnR_start(i)/ln10
988 0 : write(*,3) 'dlgR', i, (s% lnR(i) - s% lnR_start(i))/ln10
989 0 : write(*,3) 'hard_limit', i, s% delta_lgR_hard_limit
990 : end if
991 : end function check_dlgR_change
992 :
993 :
994 50 : integer function check_lgL( &
995 : s, iso_in, msg, skip_hard_limit, dt, dt_limit_ratio)
996 : type (star_info), pointer :: s
997 : integer, intent(in) :: iso_in
998 : character (len=*), intent(in) :: msg
999 : logical, intent(in) :: skip_hard_limit
1000 : real(dp), intent(in) :: dt
1001 : real(dp), intent(inout) :: dt_limit_ratio
1002 :
1003 : real(dp) :: &
1004 : new_L, max_other_L, old_L, lim, hard_lim, lgL_min, &
1005 : drop_factor, relative_limit, lgL, max_other_lgL, lgL_old, &
1006 : abs_change, relative_excess
1007 : logical, parameter :: dbg = .false.
1008 : integer :: iso
1009 : include 'formats'
1010 50 : check_lgL = keep_going
1011 :
1012 50 : iso = iso_in
1013 50 : if (iso == iprot) then ! check_lgL_power_photo_change
1014 10 : if (s% log_max_temperature < s% min_lgT_for_lgL_power_photo_limit) return
1015 0 : new_L = abs(s% power_photo)
1016 0 : max_other_L = 0d0
1017 0 : old_L = abs(s% power_photo_old)
1018 0 : lim = s% delta_lgL_power_photo_limit
1019 0 : hard_lim = s% delta_lgL_power_photo_hard_limit
1020 0 : lgL_min = s% lgL_power_photo_burn_min
1021 0 : drop_factor = s% lgL_power_photo_drop_factor
1022 0 : relative_limit = 0d0
1023 40 : else if (iso == ineut) then ! check_lgL_nuc_change
1024 10 : if (s% log_max_temperature > s% max_lgT_for_lgL_nuc_limit) return
1025 10 : new_L = s% power_nuc_burn
1026 10 : max_other_L = 0d0
1027 10 : old_L = s% power_nuc_burn_old
1028 10 : if (s% log_max_temperature > &
1029 : s% delta_lgL_nuc_at_high_T_limit_lgT_min) then
1030 0 : lim = s% delta_lgL_nuc_at_high_T_limit
1031 0 : hard_lim = s% delta_lgL_nuc_at_high_T_hard_limit
1032 : else
1033 10 : lim = s% delta_lgL_nuc_limit
1034 10 : hard_lim = s% delta_lgL_nuc_hard_limit
1035 : end if
1036 10 : lgL_min = s% lgL_nuc_burn_min
1037 10 : drop_factor = s% lgL_nuc_drop_factor
1038 10 : relative_limit = 0d0
1039 30 : else if (iso == ih1) then
1040 10 : new_L = s% power_h_burn
1041 10 : max_other_L = max(s% power_he_burn, s% power_z_burn)
1042 10 : old_L = s% power_h_burn_old
1043 10 : lim = s% delta_lgL_H_limit
1044 10 : hard_lim = s% delta_lgL_H_hard_limit
1045 10 : lgL_min = s% lgL_H_burn_min
1046 10 : drop_factor = s% lgL_H_drop_factor
1047 10 : relative_limit = s% lgL_H_burn_relative_limit
1048 20 : else if (iso == ihe4) then
1049 10 : new_L = s% power_he_burn
1050 10 : max_other_L = max(s% power_h_burn, s% power_z_burn)
1051 10 : old_L = s% power_he_burn_old
1052 10 : lim = s% delta_lgL_He_limit
1053 10 : hard_lim = s% delta_lgL_He_hard_limit
1054 10 : lgL_min = s% lgL_He_burn_min
1055 10 : drop_factor = s% lgL_He_drop_factor
1056 10 : relative_limit = s% lgL_He_burn_relative_limit
1057 10 : else if (iso == isi28) then
1058 10 : new_L = s% power_z_burn
1059 10 : max_other_L = max(s% power_h_burn, s% power_he_burn)
1060 10 : old_L = s% power_z_burn_old
1061 10 : lim = s% delta_lgL_z_limit
1062 10 : hard_lim = s% delta_lgL_z_hard_limit
1063 10 : lgL_min = s% lgL_z_burn_min
1064 10 : drop_factor = s% lgL_z_drop_factor
1065 10 : relative_limit = s% lgL_z_burn_relative_limit
1066 : else
1067 0 : call mesa_error(__FILE__,__LINE__,'bad iso arg for check_lgL')
1068 : end if
1069 :
1070 40 : if (old_L < 0d0) return
1071 :
1072 33 : lim = lim*s% time_delta_coeff
1073 33 : hard_lim = hard_lim*s% time_delta_coeff
1074 :
1075 33 : if (new_L < old_L) then
1076 9 : lim = lim*drop_factor
1077 9 : hard_lim = hard_lim*drop_factor
1078 : end if
1079 :
1080 : if (dbg) write(*,*)
1081 : if (dbg) write(*,1) trim(msg) // ' new_L', new_L
1082 : if (dbg) write(*,1) 'old_L', old_L
1083 33 : if (new_L <= 0 .or. old_L <= 0) return
1084 :
1085 30 : lgL = safe_log10(new_L)
1086 : if (dbg) write(*,1) 'lgL', lgL
1087 30 : if (lgL < lgL_min) return
1088 :
1089 20 : if (max_other_L > 0) then
1090 10 : max_other_lgL = safe_log10(max_other_L)
1091 : if (dbg) write(*,1) 'max_other_lgL', max_other_lgL
1092 10 : if (max_other_lgL - relative_limit > lgL) return
1093 : end if
1094 :
1095 20 : lgL_old = safe_log10(old_L)
1096 : if (dbg) write(*,1) 'lgL_old', lgL_old
1097 20 : abs_change = abs(lgL - lgL_old)
1098 : if (dbg) write(*,1) 'abs_change', abs_change
1099 : if (dbg) write(*,1) 'hard_lim', hard_lim
1100 20 : if (hard_lim > 0 .and. abs_change > hard_lim .and. (.not. skip_hard_limit)) then
1101 0 : if (s% report_dt_hard_limit_retries) then
1102 0 : write(*,1) trim(msg) // ' end', lgL
1103 0 : write(*,1) trim(msg) // ' start', lgL_old
1104 0 : write(*,1) trim(msg) // ' delta', lgL - lgL_old
1105 0 : write(*,1) trim(msg) // ' hard_lim', hard_lim
1106 : end if
1107 0 : check_lgL = retry
1108 0 : s% retry_message = trim(msg) // ' hard limit'
1109 0 : return
1110 : end if
1111 :
1112 : if (dbg) write(*,1) 'lim', lim
1113 20 : if (lim <= 0) return
1114 :
1115 0 : relative_excess = (abs_change - lim) / lim
1116 : if (dbg) write(*,1) 'relative_excess', relative_excess
1117 0 : dt_limit_ratio = 1d0/pow(s% timestep_dt_factor,relative_excess)
1118 0 : if (dt_limit_ratio <= 1d0) dt_limit_ratio = 0
1119 :
1120 : end function check_lgL
1121 :
1122 :
1123 10 : integer function check_lgL_H_change(s, skip_hard_limit, dt, dt_limit_ratio)
1124 : type (star_info), pointer :: s
1125 : logical, intent(in) :: skip_hard_limit
1126 : real(dp), intent(in) :: dt
1127 : real(dp), intent(inout) :: dt_limit_ratio
1128 : check_lgL_H_change = check_lgL( &
1129 10 : s, ih1, 'check_lgL_H_change', skip_hard_limit, dt, dt_limit_ratio)
1130 : end function check_lgL_H_change
1131 :
1132 :
1133 10 : integer function check_lgL_He_change(s, skip_hard_limit, dt, dt_limit_ratio)
1134 : type (star_info), pointer :: s
1135 : logical, intent(in) :: skip_hard_limit
1136 : real(dp), intent(in) :: dt
1137 : real(dp), intent(inout) :: dt_limit_ratio
1138 : check_lgL_He_change = check_lgL( &
1139 10 : s, ihe4, 'check_lgL_He_change', skip_hard_limit, dt, dt_limit_ratio)
1140 : end function check_lgL_He_change
1141 :
1142 :
1143 10 : integer function check_lgL_z_change(s, skip_hard_limit, dt, dt_limit_ratio)
1144 : type (star_info), pointer :: s
1145 : logical, intent(in) :: skip_hard_limit
1146 : real(dp), intent(in) :: dt
1147 : real(dp), intent(inout) :: dt_limit_ratio
1148 : check_lgL_z_change = check_lgL( &
1149 10 : s, isi28, 'check_lgL_z_change', skip_hard_limit, dt, dt_limit_ratio)
1150 : end function check_lgL_z_change
1151 :
1152 :
1153 10 : integer function check_lgL_power_photo_change(s, skip_hard_limit, dt, dt_limit_ratio)
1154 : type (star_info), pointer :: s
1155 : logical, intent(in) :: skip_hard_limit
1156 : real(dp), intent(in) :: dt
1157 : real(dp), intent(inout) :: dt_limit_ratio
1158 : check_lgL_power_photo_change = check_lgL( &
1159 10 : s, iprot, 'check_lgL_power_photo_change', skip_hard_limit, dt, dt_limit_ratio)
1160 : end function check_lgL_power_photo_change
1161 :
1162 :
1163 10 : integer function check_lgL_nuc_change(s, skip_hard_limit, dt, dt_limit_ratio)
1164 : type (star_info), pointer :: s
1165 : logical, intent(in) :: skip_hard_limit
1166 : real(dp), intent(in) :: dt
1167 : real(dp), intent(inout) :: dt_limit_ratio
1168 : include 'formats'
1169 : check_lgL_nuc_change = check_lgL( &
1170 10 : s, ineut, 'check_lgL_nuc_change', skip_hard_limit, dt, dt_limit_ratio)
1171 : end function check_lgL_nuc_change
1172 :
1173 :
1174 10 : integer function check_lgL_nuc_cat_change( &
1175 : s, n_mix_bdy, mix_bdy_q, skip_hard_limit, dt_limit_ratio)
1176 : use rates_def
1177 : use num_lib, only: binary_search
1178 : type (star_info), pointer :: s
1179 : logical, intent(in) :: skip_hard_limit
1180 : integer, intent(in) :: n_mix_bdy
1181 : real(dp), intent(in), pointer :: mix_bdy_q(:)
1182 : real(dp), intent(inout) :: dt_limit_ratio
1183 :
1184 : integer :: k, max_j, max_k, bdy
1185 : real(dp) :: max_lgL_diff, relative_excess, max_diff, cat_burn_min, &
1186 : max_luminosity, max_luminosity_start
1187 :
1188 : include 'formats'
1189 :
1190 10 : check_lgL_nuc_cat_change = keep_going
1191 :
1192 10 : if (s% delta_lgL_nuc_cat_limit <= 0 .and. s% delta_lgL_nuc_cat_hard_limit <= 0) return
1193 :
1194 0 : cat_burn_min = exp10(s% lgL_nuc_cat_burn_min)
1195 0 : max_diff = 0
1196 0 : max_j = -1
1197 0 : max_k = -1
1198 0 : bdy = -1
1199 :
1200 0 : do k = 1, s% nz
1201 :
1202 : ! find the nearest mixing boundary
1203 0 : bdy = binary_search(n_mix_bdy, mix_bdy_q, bdy, s% q(k))
1204 : ! don't check cells near a mixing boundary
1205 0 : if (bdy > 0) then
1206 0 : if (abs(s% q(k) - mix_bdy_q(bdy)) < s% lgL_nuc_mix_dist_limit) cycle
1207 : end if
1208 :
1209 0 : call do1_category(ipp,k)
1210 0 : call do1_category(icno,k)
1211 0 : call do1_category(i3alf,k)
1212 0 : call do1_category(i_burn_c,k)
1213 0 : call do1_category(i_burn_n,k)
1214 0 : call do1_category(i_burn_o,k)
1215 0 : call do1_category(i_burn_ne,k)
1216 0 : call do1_category(i_burn_na,k)
1217 0 : call do1_category(i_burn_mg,k)
1218 0 : call do1_category(i_burn_si,k)
1219 0 : call do1_category(i_burn_s,k)
1220 0 : call do1_category(i_burn_ar,k)
1221 0 : call do1_category(i_burn_ca,k)
1222 0 : call do1_category(i_burn_ti,k)
1223 0 : call do1_category(i_burn_cr,k)
1224 0 : call do1_category(i_burn_fe,k)
1225 0 : call do1_category(icc,k)
1226 0 : call do1_category(ico,k)
1227 0 : call do1_category(ioo,k)
1228 :
1229 : end do
1230 :
1231 :
1232 0 : if (max_diff <= 0) return
1233 :
1234 0 : max_lgL_diff = log10(max_diff/Lsun)
1235 0 : s% Tlim_dlgL_nuc_category = max_j
1236 0 : s% Tlim_dlgL_nuc_cell = max_k
1237 :
1238 : check_lgL_nuc_cat_change = check_change(s, max_lgL_diff, &
1239 : s% delta_lgL_nuc_cat_limit, s% delta_lgL_nuc_cat_hard_limit, &
1240 0 : max_k, 'delta_lgL_nuc_cat', skip_hard_limit, dt_limit_ratio, relative_excess)
1241 0 : if (check_lgL_nuc_cat_change /= keep_going .and. s% report_dt_hard_limit_retries) then
1242 0 : write(*,1) 'max_luminosity ' // trim(category_name(max_j)), max_luminosity
1243 0 : write(*,1) 'max_luminosity_start ' // trim(category_name(max_j)), max_luminosity_start
1244 : end if
1245 :
1246 : contains
1247 :
1248 0 : subroutine do1_category(j, k)
1249 : integer, intent(in) :: j, k
1250 : real(dp) :: diff, abs_diff
1251 0 : if (s% luminosity_by_category(j,k) < cat_burn_min) return
1252 0 : if (s% luminosity_by_category_start(j,k) < cat_burn_min) return
1253 0 : diff = s% luminosity_by_category(j,k) - s% luminosity_by_category_start(j,k)
1254 0 : abs_diff = abs(diff)
1255 0 : if (abs_diff <= max_diff) return
1256 0 : max_luminosity = s% luminosity_by_category(j,k)
1257 0 : max_luminosity_start = s% luminosity_by_category_start(j,k)
1258 0 : max_diff = abs_diff
1259 0 : max_j = j
1260 0 : max_k = k
1261 : end subroutine do1_category
1262 :
1263 : end function check_lgL_nuc_cat_change
1264 :
1265 :
1266 10 : integer function check_dlgTeff_change(s, skip_hard_limit, dt, dt_limit_ratio)
1267 : type (star_info), pointer :: s
1268 : logical, intent(in) :: skip_hard_limit
1269 : real(dp), intent(in) :: dt
1270 : real(dp), intent(inout) :: dt_limit_ratio
1271 : real(dp) :: relative_excess
1272 : include 'formats'
1273 10 : check_dlgTeff_change = keep_going
1274 10 : dt_limit_ratio = 0d0
1275 10 : if (s% doing_relax .or. s% Teff_old <= 0 .or. s% Teff <= 0) return
1276 : check_dlgTeff_change = check_change(s, safe_log10(s% Teff/s% Teff_old), &
1277 : s% delta_lgTeff_limit, s% delta_lgTeff_hard_limit, &
1278 10 : 0, 'delta_lgTeff', skip_hard_limit, dt_limit_ratio, relative_excess)
1279 10 : if (check_dlgTeff_change /= keep_going .and. s% report_dt_hard_limit_retries) then
1280 0 : write(*,1) 'lgTeff', safe_log10(s% Teff)
1281 0 : write(*,1) 'lgTeff_old', safe_log10(s% Teff_old)
1282 : end if
1283 : end function check_dlgTeff_change
1284 :
1285 :
1286 10 : integer function check_dYe_highT_change(s, skip_hard_limit, dt_limit_ratio) ! check max change in Ye
1287 : type (star_info), pointer :: s
1288 : logical, intent(in) :: skip_hard_limit
1289 : real(dp), intent(inout) :: dt_limit_ratio
1290 : real(dp) :: relative_excess, max_diff, ye_diff, T_limit
1291 : integer :: i, k
1292 : include 'formats'
1293 10 : check_dYe_highT_change = keep_going
1294 10 : dt_limit_ratio = 0d0
1295 10 : if (s% doing_relax) return
1296 10 : i = 0
1297 10 : max_diff = 0
1298 10 : T_limit = s% minT_for_highT_Ye_limit
1299 10633 : do k=1, s% nz
1300 10623 : if (s% T(k) < T_limit) cycle
1301 0 : ye_diff = abs(s% ye(k) - s% ye_start(k))
1302 0 : if (ye_diff <= max_diff) cycle
1303 0 : max_diff = ye_diff
1304 10633 : i = k
1305 : end do
1306 : check_dYe_highT_change = check_change(s, max_diff, &
1307 : s% delta_Ye_highT_limit, s% delta_Ye_highT_hard_limit, &
1308 10 : i, 'delta_Ye_highT', .false., dt_limit_ratio, relative_excess)
1309 10 : if (check_dYe_highT_change /= keep_going .and. s% report_dt_hard_limit_retries) then
1310 0 : write(*,2) 'ye', i, s% ye(i)
1311 0 : write(*,2) 'ye_start', i, s% ye_start(i)
1312 : end if
1313 : end function check_dYe_highT_change
1314 :
1315 :
1316 10 : integer function check_dlgT_max_change(s, skip_hard_limit, dt, dt_limit_ratio)
1317 : type (star_info), pointer :: s
1318 : logical, intent(in) :: skip_hard_limit
1319 : real(dp), intent(in) :: dt
1320 : real(dp), intent(inout) :: dt_limit_ratio
1321 : real(dp) :: relative_excess, change, lnTmax, lnTmax_start
1322 : integer :: lnTmax_k
1323 : include 'formats'
1324 10 : check_dlgT_max_change = keep_going
1325 10 : dt_limit_ratio = 0d0
1326 10 : if (s% doing_relax) return
1327 10 : if (s% delta_lgT_max_limit_lgT_min < 0d0) return
1328 10 : if (s% delta_lgT_max_limit_only_after_near_zams) then
1329 0 : if (s% X(s% nz) > 0.1d0 .and. &
1330 : s% L_nuc_burn_total/s% L_phot < s% Lnuc_div_L_zams_limit ) return
1331 : end if
1332 10633 : lnTmax_k = maxloc(s% lnT(1:s% nz),dim=1)
1333 10 : lnTmax = s% lnT(lnTmax_k)
1334 10 : if (lnTmax < s% delta_lgT_max_limit_lgT_min*ln10) return
1335 0 : lnTmax_start = maxval(s% lnT_start(1:s% nz))
1336 0 : change = (lnTmax - lnTmax_start)/ln10
1337 : check_dlgT_max_change = check_change(s, change, &
1338 : s% delta_lgT_max_limit, s% delta_lgT_max_hard_limit, &
1339 0 : lnTmax_k, 'delta_lgT_max', skip_hard_limit, dt_limit_ratio, relative_excess)
1340 0 : if (check_dlgT_max_change /= keep_going .and. s% report_dt_hard_limit_retries) then
1341 0 : write(*,2) 'lgT_max', lnTmax_k, lnTmax/ln10
1342 0 : write(*,2) 'lgT_max_old', lnTmax_k, lnTmax_start/ln10
1343 : end if
1344 : end function check_dlgT_max_change
1345 :
1346 :
1347 10 : integer function check_dlgT_max_at_high_T_change(s, skip_hard_limit, dt, dt_limit_ratio)
1348 : type (star_info), pointer :: s
1349 : logical, intent(in) :: skip_hard_limit
1350 : real(dp), intent(in) :: dt
1351 : real(dp), intent(inout) :: dt_limit_ratio
1352 : real(dp) :: relative_excess, change, lnTmax_at_high_T, lnTmax_at_high_T_start
1353 : integer :: lnTmax_k
1354 : include 'formats'
1355 10 : check_dlgT_max_at_high_T_change = keep_going
1356 10 : dt_limit_ratio = 0d0
1357 10 : if (s% doing_relax) return
1358 10 : if (s% delta_lgT_max_at_high_T_limit_lgT_min < 0d0) return
1359 0 : lnTmax_k = maxloc(s% lnT(1:s% nz),dim=1)
1360 0 : lnTmax_at_high_T = s% lnT(lnTmax_k)
1361 0 : if (lnTmax_at_high_T < s% delta_lgT_max_at_high_T_limit_lgT_min*ln10) return
1362 0 : lnTmax_at_high_T_start = maxval(s% lnT_start(1:s% nz))
1363 0 : change = (lnTmax_at_high_T - lnTmax_at_high_T_start)/ln10
1364 : check_dlgT_max_at_high_T_change = check_change(s, change, &
1365 : s% delta_lgT_max_at_high_T_limit, s% delta_lgT_max_at_high_T_hard_limit, &
1366 0 : lnTmax_k, 'delta_lgT_max_at_high_T', skip_hard_limit, dt_limit_ratio, relative_excess)
1367 0 : if (check_dlgT_max_at_high_T_change /= keep_going .and. s% report_dt_hard_limit_retries) then
1368 0 : write(*,2) 'lgT_max', lnTmax_k, lnTmax_at_high_T/ln10
1369 0 : write(*,2) 'lgT_max_old', lnTmax_k, lnTmax_at_high_T_start/ln10
1370 : end if
1371 : end function check_dlgT_max_at_high_T_change
1372 :
1373 :
1374 10 : integer function check_dlgT_cntr_change(s, skip_hard_limit, dt, dt_limit_ratio)
1375 : type (star_info), pointer :: s
1376 : logical, intent(in) :: skip_hard_limit
1377 : real(dp), intent(in) :: dt
1378 : real(dp), intent(inout) :: dt_limit_ratio
1379 : real(dp) :: relative_excess, change
1380 : include 'formats'
1381 10 : check_dlgT_cntr_change = keep_going
1382 10 : dt_limit_ratio = 0d0
1383 10 : if (s% doing_relax) return
1384 10 : if (s% delta_lgT_cntr_limit_only_after_near_zams) then
1385 10 : if (s% X(s% nz) > 0.1d0 .and. &
1386 : s% L_nuc_burn_total/s% L_phot < s% Lnuc_div_L_zams_limit) return
1387 : end if
1388 10 : change = (s% lnT(s% nz) - s% lnT_start(s% nz))/ln10
1389 : check_dlgT_cntr_change = check_change(s, change, &
1390 : s% delta_lgT_cntr_limit, s% delta_lgT_cntr_hard_limit, &
1391 10 : s% nz, 'delta_lgT_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1392 10 : if (check_dlgT_cntr_change /= keep_going .and. s% report_dt_hard_limit_retries) then
1393 0 : write(*,1) 'lgT_cntr', s% lnT(s% nz)/ln10
1394 0 : write(*,1) 'lgT_cntr_old', s% lnT_start(s% nz)/ln10
1395 : end if
1396 : end function check_dlgT_cntr_change
1397 :
1398 :
1399 10 : integer function check_dlgP_cntr_change(s, skip_hard_limit, dt, dt_limit_ratio)
1400 : type (star_info), pointer :: s
1401 : logical, intent(in) :: skip_hard_limit
1402 : real(dp), intent(in) :: dt
1403 : real(dp), intent(inout) :: dt_limit_ratio
1404 : real(dp) :: relative_excess, change
1405 : include 'formats'
1406 10 : check_dlgP_cntr_change = keep_going
1407 10 : dt_limit_ratio = 0d0
1408 10 : if (s% doing_relax) return
1409 10 : change = (s% lnPeos(s% nz) - s% lnPeos_start(s% nz))/ln10
1410 : check_dlgP_cntr_change = check_change(s, change, &
1411 : s% delta_lgP_cntr_limit, s% delta_lgP_cntr_hard_limit, &
1412 10 : s% nz, 'delta_lgP_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1413 10 : if (check_dlgP_cntr_change /= keep_going .and. s% report_dt_hard_limit_retries) then
1414 0 : write(*,1) 'lgP_cntr', s% lnPeos(s% nz)/ln10
1415 0 : write(*,1) 'lgP_cntr_old', s% lnPeos_start(s% nz)/ln10
1416 : end if
1417 : end function check_dlgP_cntr_change
1418 :
1419 :
1420 10 : integer function check_dlgRho_cntr_change(s, skip_hard_limit, dt, dt_limit_ratio)
1421 : type (star_info), pointer :: s
1422 : logical, intent(in) :: skip_hard_limit
1423 : real(dp), intent(in) :: dt
1424 : real(dp), intent(inout) :: dt_limit_ratio
1425 : real(dp) :: relative_excess, dlgRho_cntr
1426 : integer :: nz
1427 : include 'formats'
1428 10 : check_dlgRho_cntr_change = keep_going
1429 10 : dt_limit_ratio = 0d0
1430 10 : if (s% doing_relax) return
1431 10 : nz = s% nz
1432 10 : dlgRho_cntr = (s% lnd(nz) - s% lnd_start(nz))/ln10
1433 : check_dlgRho_cntr_change = check_change(s, dlgRho_cntr, &
1434 : s% delta_lgRho_cntr_limit, s% delta_lgRho_cntr_hard_limit, &
1435 10 : nz, 'delta_lgRho_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1436 10 : if (check_dlgRho_cntr_change /= keep_going .and. s% report_dt_hard_limit_retries) then
1437 0 : write(*,1) 'lgRho_cntr', s% lnd(s% nz)/ln10
1438 0 : write(*,1) 'lgRho_cntr_old', s% lnd_start(s% nz)/ln10
1439 : end if
1440 : end function check_dlgRho_cntr_change
1441 :
1442 :
1443 10 : integer function check_dlog_eps_nuc_change(s, skip_hard_limit, dt, dt_limit_ratio)
1444 : type (star_info), pointer :: s
1445 : logical, intent(in) :: skip_hard_limit
1446 : real(dp), intent(in) :: dt
1447 : real(dp), intent(inout) :: dt_limit_ratio
1448 : real(dp) :: relative_excess, max_ratio, ratio, delta, &
1449 : limit_ratio, delta_r
1450 : integer :: k, j, nz, k_max
1451 : include 'formats'
1452 10 : check_dlog_eps_nuc_change = keep_going
1453 10 : nz = s% nz
1454 10 : limit_ratio = exp10(s% delta_log_eps_nuc_limit)
1455 10 : max_ratio = limit_ratio
1456 10 : k_max = 0
1457 10633 : zoneloop: do k=1,nz
1458 10623 : if (s% eps_nuc_start(k) < 1) cycle zoneloop
1459 2295 : ratio = s% eps_nuc(k)/s% eps_nuc_start(k)
1460 : if (s% mixing_type(k) /= convective_mixing .and. &
1461 2295 : s% mixing_type(min(nz,k+1)) /= convective_mixing .and. &
1462 10 : ratio > max_ratio) then
1463 674 : do j = 1, s% num_conv_boundaries
1464 541 : delta_r = abs(s% r(s% conv_bdy_loc(j)) - s% r(k))
1465 674 : if (delta_r <= s% scale_height(k)) then
1466 : cycle zoneloop ! skip ones that are too close to convection zone
1467 : end if
1468 : end do
1469 133 : max_ratio = ratio
1470 133 : k_max = k
1471 : end if
1472 : end do zoneloop
1473 10 : if (k_max > 0) then
1474 10 : delta = log10(max_ratio)
1475 : else
1476 0 : delta = 0
1477 : end if
1478 : check_dlog_eps_nuc_change = check_change(s, delta, &
1479 : s% delta_log_eps_nuc_limit, s% delta_log_eps_nuc_hard_limit, &
1480 : k_max, 'delta_log_eps_nuc', &
1481 10 : skip_hard_limit, dt_limit_ratio, relative_excess)
1482 10 : if (check_dlog_eps_nuc_change /= keep_going .and. s% report_dt_hard_limit_retries) then
1483 0 : write(*,2) 'log_eps_nuc', k_max, safe_log10(abs(s% eps_nuc(k_max)))
1484 0 : write(*,2) 'log_eps_nuc_old', k_max, safe_log10(abs(s% eps_nuc_start(k_max)))
1485 : end if
1486 10 : end function check_dlog_eps_nuc_change
1487 :
1488 :
1489 10 : integer function check_dX_div_X_cntr(s, skip_hard_limit, dt_limit_ratio)
1490 : type (star_info), pointer :: s
1491 : logical, intent(in) :: skip_hard_limit
1492 : real(dp), intent(inout) :: dt_limit_ratio
1493 : real(dp) :: relative_excess, max_abs_dX_div_X, X, dX, abs_dX_div_X
1494 : integer :: j, nz, j_max
1495 : include 'formats'
1496 10 : check_dX_div_X_cntr = keep_going
1497 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1498 10 : nz = s% nz
1499 10 : max_abs_dX_div_X = -1
1500 90 : do j=1,s% species
1501 80 : X = s% xa(j,nz)
1502 80 : if (X > s% delta_dX_div_X_cntr_max) cycle
1503 0 : if (X < s% delta_dX_div_X_cntr_min) cycle
1504 0 : if (X <= 0d0) cycle
1505 0 : dX = X - s% xa_old(j,nz)
1506 0 : if (s% delta_dX_div_X_drop_only .and. dX > 0) cycle
1507 0 : abs_dX_div_X = abs(dX/X)
1508 10 : if (abs_dX_div_X > max_abs_dX_div_X) then
1509 0 : max_abs_dX_div_X = abs_dX_div_X
1510 0 : j_max = j
1511 : end if
1512 : end do
1513 10 : if (max_abs_dX_div_X <= 0d0) return
1514 : check_dX_div_X_cntr = check_change(s, max_abs_dX_div_X, &
1515 : s% delta_dX_div_X_cntr_limit, s% delta_dX_div_X_cntr_hard_limit, &
1516 0 : 0, 'delta_dX_div_X_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1517 0 : if (check_dX_div_X_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1518 0 : write(*,1) chem_isos% name(s% chem_id(j_max)) // ' X', s% xa(j_max,nz)
1519 0 : write(*,1) chem_isos% name(s% chem_id(j_max)) // ' X old', s% xa_old(j_max,nz)
1520 : end if
1521 : end function check_dX_div_X_cntr
1522 :
1523 :
1524 10 : integer function check_lg_XH_cntr(s, skip_hard_limit, dt_limit_ratio)
1525 : type (star_info), pointer :: s
1526 : logical, intent(in) :: skip_hard_limit
1527 : real(dp), intent(inout) :: dt_limit_ratio
1528 : real(dp) :: relative_excess, lg_XH_cntr, lg_XH_cntr_old
1529 : integer :: h1, nz
1530 : include 'formats'
1531 10 : check_lg_XH_cntr = keep_going
1532 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1533 10 : h1 = s% net_iso(ih1)
1534 10 : if (h1 == 0) return
1535 10 : nz = s% nz
1536 10 : if (s% xa(h1,nz) < 1d-10) return
1537 10 : lg_XH_cntr = log10(s% xa(h1,nz))
1538 10 : if (lg_XH_cntr > s% delta_lg_XH_cntr_max) return
1539 10 : if (lg_XH_cntr < s% delta_lg_XH_cntr_min) return
1540 10 : lg_XH_cntr_old = safe_log10(s% xa_old(h1,nz))
1541 10 : if (s% delta_lg_XH_drop_only .and. lg_XH_cntr >= lg_XH_cntr_old) return
1542 : check_lg_XH_cntr = check_change(s, lg_XH_cntr - lg_XH_cntr_old, &
1543 : s% delta_lg_XH_cntr_limit, s% delta_lg_XH_cntr_hard_limit, &
1544 10 : nz, 'delta_lg_XH_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1545 10 : if (check_lg_XH_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1546 0 : write(*,1) 'lg_XH_cntr', lg_XH_cntr
1547 0 : write(*,1) 'lg_XH_cntr_old', lg_XH_cntr_old
1548 0 : write(*,1) 'delta', lg_XH_cntr - lg_XH_cntr_old
1549 : end if
1550 : end function check_lg_XH_cntr
1551 :
1552 :
1553 10 : integer function check_lg_XHe_cntr(s, skip_hard_limit, dt_limit_ratio)
1554 : type (star_info), pointer :: s
1555 : logical, intent(in) :: skip_hard_limit
1556 : real(dp), intent(inout) :: dt_limit_ratio
1557 : real(dp) :: relative_excess, lg_XHe_cntr, lg_XHe_cntr_old
1558 : integer :: h1, he4, nz
1559 : real(dp) :: xh1, xhe4
1560 : include 'formats'
1561 10 : check_lg_XHe_cntr = keep_going
1562 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1563 10 : h1 = s% net_iso(ih1)
1564 10 : he4 = s% net_iso(ihe4)
1565 10 : if (h1 == 0 .or. he4 == 0) return
1566 10 : nz = s% nz
1567 10 : xh1 = s% xa(h1,nz)
1568 10 : xhe4 = s% xa(he4,nz)
1569 10 : if (xhe4 < max(xh1, 1d-10)) return
1570 0 : lg_XHe_cntr = log10(xhe4)
1571 0 : if (lg_XHe_cntr > s% delta_lg_XHe_cntr_max) return
1572 0 : if (lg_XHe_cntr < s% delta_lg_XHe_cntr_min) return
1573 0 : lg_XHe_cntr_old = safe_log10(s% xa_old(he4,nz))
1574 0 : if (s% delta_lg_XHe_drop_only .and. lg_XHe_cntr >= lg_XHe_cntr_old) return
1575 : check_lg_XHe_cntr = check_change(s, lg_XHe_cntr - lg_XHe_cntr_old, &
1576 : s% delta_lg_XHe_cntr_limit, s% delta_lg_XHe_cntr_hard_limit, &
1577 0 : nz, 'delta_lg_XHe_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1578 0 : if (check_lg_XHe_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1579 0 : write(*,1) 'lg_XHe_cntr', lg_XHe_cntr
1580 0 : write(*,1) 'lg_XHe_cntr_old', lg_XHe_cntr_old
1581 0 : write(*,1) 'delta', lg_XHe_cntr - lg_XHe_cntr_old
1582 : end if
1583 : end function check_lg_XHe_cntr
1584 :
1585 :
1586 10 : integer function check_lg_XC_cntr(s, skip_hard_limit, dt_limit_ratio)
1587 : type (star_info), pointer :: s
1588 : logical, intent(in) :: skip_hard_limit
1589 : real(dp), intent(inout) :: dt_limit_ratio
1590 : real(dp) :: relative_excess, lg_XC_cntr, lg_XC_cntr_old
1591 : integer :: h1, he4, c12, nz
1592 : real(dp) :: xh1, xhe4, xc12
1593 : include 'formats'
1594 10 : check_lg_XC_cntr = keep_going
1595 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1596 10 : h1 = s% net_iso(ih1)
1597 10 : he4 = s% net_iso(ihe4)
1598 10 : c12 = s% net_iso(ic12)
1599 10 : if (h1 == 0 .or. he4 == 0 .or. c12 == 0) return
1600 10 : nz = s% nz
1601 10 : xh1 = s% xa(h1,nz)
1602 10 : xhe4 = s% xa(he4,nz)
1603 10 : xc12 = s% xa(c12,nz)
1604 10 : if (xc12 < max(xh1, xhe4, 1d-10)) return
1605 0 : if (s% xa(c12,nz) < 1d-10) return
1606 0 : lg_XC_cntr = log10(xc12)
1607 0 : if (lg_XC_cntr > s% delta_lg_XC_cntr_max) return
1608 0 : if (lg_XC_cntr < s% delta_lg_XC_cntr_min) return
1609 0 : lg_XC_cntr_old = safe_log10(s% xa_old(c12,nz))
1610 0 : if (s% delta_lg_XC_drop_only .and. lg_XC_cntr >= lg_XC_cntr_old) return
1611 : check_lg_XC_cntr = check_change(s, lg_XC_cntr - lg_XC_cntr_old, &
1612 : s% delta_lg_XC_cntr_limit, s% delta_lg_XC_cntr_hard_limit, &
1613 0 : nz, 'delta_lg_XC_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1614 0 : if (check_lg_XC_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1615 0 : write(*,1) 'lg_XC_cntr', lg_XC_cntr
1616 0 : write(*,1) 'lg_XC_cntr_old', lg_XC_cntr_old
1617 0 : write(*,1) 'delta', lg_XC_cntr - lg_XC_cntr_old
1618 : end if
1619 : end function check_lg_XC_cntr
1620 :
1621 :
1622 10 : integer function check_lg_XNe_cntr(s, skip_hard_limit, dt_limit_ratio)
1623 : type (star_info), pointer :: s
1624 : logical, intent(in) :: skip_hard_limit
1625 : real(dp), intent(inout) :: dt_limit_ratio
1626 : real(dp) :: relative_excess, lg_XNe_cntr, lg_XNe_cntr_old
1627 : integer :: h1, he4, c12, ne20, nz
1628 : real(dp) :: xh1, xhe4, xc12, xne20
1629 : include 'formats'
1630 10 : check_lg_XNe_cntr = keep_going
1631 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1632 10 : h1 = s% net_iso(ih1)
1633 10 : he4 = s% net_iso(ihe4)
1634 10 : c12 = s% net_iso(ic12)
1635 10 : ne20 = s% net_iso(ine20)
1636 10 : if (h1 == 0 .or. he4 == 0 .or. c12 == 0 .or. ne20 == 0) return
1637 10 : nz = s% nz
1638 10 : xh1 = s% xa(h1,nz)
1639 10 : xhe4 = s% xa(he4,nz)
1640 10 : xc12 = s% xa(c12,nz)
1641 10 : xne20 = s% xa(ne20,nz)
1642 10 : if (xne20 < max(xh1, xhe4, xc12, 1d-10)) return
1643 0 : lg_XNe_cntr = log10(xne20)
1644 0 : if (lg_XNe_cntr > s% delta_lg_XNe_cntr_max) return
1645 0 : if (lg_XNe_cntr < s% delta_lg_XNe_cntr_min) return
1646 0 : lg_XNe_cntr_old = safe_log10(s% xa_old(ne20,nz))
1647 0 : if (s% delta_lg_XNe_drop_only .and. lg_XNe_cntr >= lg_XNe_cntr_old) return
1648 : check_lg_XNe_cntr = check_change(s, lg_XNe_cntr - lg_XNe_cntr_old, &
1649 : s% delta_lg_XNe_cntr_limit, s% delta_lg_XNe_cntr_hard_limit, &
1650 0 : nz, 'delta_lg_XNe_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1651 0 : if (check_lg_XNe_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1652 0 : write(*,1) 'lg_XNe_cntr', lg_XNe_cntr
1653 0 : write(*,1) 'lg_XNe_cntr_old', lg_XNe_cntr_old
1654 0 : write(*,1) 'delta', lg_XNe_cntr - lg_XNe_cntr_old
1655 : end if
1656 : end function check_lg_XNe_cntr
1657 :
1658 :
1659 10 : integer function check_lg_XO_cntr(s, skip_hard_limit, dt_limit_ratio)
1660 : type (star_info), pointer :: s
1661 : logical, intent(in) :: skip_hard_limit
1662 : real(dp), intent(inout) :: dt_limit_ratio
1663 : real(dp) :: relative_excess, lg_XO_cntr, lg_XO_cntr_old
1664 : integer :: h1, he4, c12, o16, nz
1665 : real(dp) :: xh1, xhe4, xc12, xo16
1666 : include 'formats'
1667 10 : check_lg_XO_cntr = keep_going
1668 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1669 10 : h1 = s% net_iso(ih1)
1670 10 : he4 = s% net_iso(ihe4)
1671 10 : c12 = s% net_iso(ic12)
1672 10 : o16 = s% net_iso(io16)
1673 10 : if (h1 == 0 .or. he4 == 0 .or. c12 == 0 .or. o16 == 0) return
1674 10 : nz = s% nz
1675 10 : xh1 = s% xa(h1,nz)
1676 10 : xhe4 = s% xa(he4,nz)
1677 10 : xc12 = s% xa(c12,nz)
1678 10 : xo16 = s% xa(o16,nz)
1679 10 : if (xo16 < max(xh1, xhe4, xc12, 1d-10)) return
1680 0 : lg_XO_cntr = log10(xo16)
1681 0 : if (lg_XO_cntr > s% delta_lg_XO_cntr_max) return
1682 0 : if (lg_XO_cntr < s% delta_lg_XO_cntr_min) return
1683 0 : lg_XO_cntr_old = safe_log10(s% xa_old(o16,nz))
1684 0 : if (s% delta_lg_XO_drop_only .and. lg_XO_cntr >= lg_XO_cntr_old) return
1685 : check_lg_XO_cntr = check_change(s, lg_XO_cntr - lg_XO_cntr_old, &
1686 : s% delta_lg_XO_cntr_limit, s% delta_lg_XO_cntr_hard_limit, &
1687 0 : nz, 'delta_lg_XO_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1688 0 : if (check_lg_XO_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1689 0 : write(*,1) 'lg_XO_cntr', lg_XO_cntr
1690 0 : write(*,1) 'lg_XO_cntr_old', lg_XO_cntr_old
1691 0 : write(*,1) 'delta', lg_XO_cntr - lg_XO_cntr_old
1692 : end if
1693 : end function check_lg_XO_cntr
1694 :
1695 :
1696 10 : integer function check_lg_XSi_cntr(s, skip_hard_limit, dt_limit_ratio)
1697 : type (star_info), pointer :: s
1698 : logical, intent(in) :: skip_hard_limit
1699 : real(dp), intent(inout) :: dt_limit_ratio
1700 : real(dp) :: relative_excess, lg_XSi_cntr, lg_XSi_cntr_old
1701 : integer :: h1, he4, c12, si28, nz
1702 : real(dp) :: xh1, xhe4, xc12, xsi28
1703 : include 'formats'
1704 10 : check_lg_XSi_cntr = keep_going
1705 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1706 10 : h1 = s% net_iso(ih1)
1707 10 : he4 = s% net_iso(ihe4)
1708 10 : c12 = s% net_iso(ic12)
1709 10 : si28 = s% net_iso(isi28)
1710 10 : if (h1 == 0 .or. he4 == 0 .or. c12 == 0 .or. si28 == 0) return
1711 0 : nz = s% nz
1712 0 : xh1 = s% xa(h1,nz)
1713 0 : xhe4 = s% xa(he4,nz)
1714 0 : xc12 = s% xa(c12,nz)
1715 0 : xsi28 = s% xa(si28,nz)
1716 0 : if (xsi28 < max(xh1, xhe4, xc12, 1d-10)) return
1717 0 : lg_XSi_cntr = log10(xsi28)
1718 0 : if (lg_XSi_cntr > s% delta_lg_XSi_cntr_max) return
1719 0 : if (lg_XSi_cntr < s% delta_lg_XSi_cntr_min) return
1720 0 : lg_XSi_cntr_old = safe_log10(s% xa_old(si28,nz))
1721 0 : if (s% delta_lg_XSi_drop_only .and. lg_XSi_cntr >= lg_XSi_cntr_old) return
1722 : check_lg_XSi_cntr = check_change(s, lg_XSi_cntr - lg_XSi_cntr_old, &
1723 : s% delta_lg_XSi_cntr_limit, s% delta_lg_XSi_cntr_hard_limit, &
1724 0 : nz, 'delta_lg_XSi_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1725 0 : if (check_lg_XSi_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1726 0 : write(*,1) 'lg_XSi_cntr', lg_XSi_cntr
1727 0 : write(*,1) 'lg_XSi_cntr_old', lg_XSi_cntr_old
1728 0 : write(*,1) 'delta', lg_XSi_cntr - lg_XSi_cntr_old
1729 : end if
1730 : end function check_lg_XSi_cntr
1731 :
1732 :
1733 10 : integer function check_XH_cntr(s, skip_hard_limit, dt_limit_ratio)
1734 : type (star_info), pointer :: s
1735 : logical, intent(in) :: skip_hard_limit
1736 : real(dp), intent(inout) :: dt_limit_ratio
1737 : real(dp) :: relative_excess, XH_cntr, XH_cntr_old
1738 : integer :: h1, nz
1739 : include 'formats'
1740 10 : check_XH_cntr = keep_going
1741 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1742 10 : h1 = s% net_iso(ih1)
1743 10 : if (h1 == 0) return
1744 10 : nz = s% nz
1745 10 : XH_cntr = s% xa(h1,nz)
1746 10 : XH_cntr_old = s% xa_old(h1,nz)
1747 10 : if (s% delta_XH_drop_only .and. XH_cntr >= XH_cntr_old) return
1748 : check_XH_cntr = check_change(s, XH_cntr - XH_cntr_old, &
1749 : s% delta_XH_cntr_limit, s% delta_XH_cntr_hard_limit, &
1750 10 : nz, 'delta_XH_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1751 10 : if (check_XH_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1752 0 : write(*,1) 'XH_cntr', XH_cntr
1753 0 : write(*,1) 'XH_cntr_old', XH_cntr_old
1754 0 : write(*,1) 'delta', XH_cntr - XH_cntr_old
1755 : end if
1756 : end function check_XH_cntr
1757 :
1758 :
1759 10 : integer function check_XHe_cntr(s, skip_hard_limit, dt_limit_ratio)
1760 : type (star_info), pointer :: s
1761 : logical, intent(in) :: skip_hard_limit
1762 : real(dp), intent(inout) :: dt_limit_ratio
1763 : real(dp) :: relative_excess, XHe_cntr, XHe_cntr_old
1764 : integer :: he4, nz
1765 : include 'formats'
1766 10 : check_XHe_cntr = keep_going
1767 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1768 10 : he4 = s% net_iso(ihe4)
1769 10 : if (he4 == 0) return
1770 10 : nz = s% nz
1771 10 : XHe_cntr = s% xa(he4,nz)
1772 10 : XHe_cntr_old = s% xa_old(he4,nz)
1773 10 : if (s% delta_XHe_drop_only .and. XHe_cntr >= XHe_cntr_old) return
1774 : check_XHe_cntr = check_change(s, XHe_cntr - XHe_cntr_old, &
1775 : s% delta_XHe_cntr_limit, s% delta_XHe_cntr_hard_limit, &
1776 10 : nz, 'delta_XHe_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1777 10 : if (check_XHe_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1778 0 : write(*,1) 'XHe_cntr', XHe_cntr
1779 0 : write(*,1) 'XHe_cntr_old', XHe_cntr_old
1780 0 : write(*,1) 'delta', XHe_cntr - XHe_cntr_old
1781 : end if
1782 : end function check_XHe_cntr
1783 :
1784 :
1785 10 : integer function check_XC_cntr(s, skip_hard_limit, dt_limit_ratio)
1786 : type (star_info), pointer :: s
1787 : logical, intent(in) :: skip_hard_limit
1788 : real(dp), intent(inout) :: dt_limit_ratio
1789 : real(dp) :: relative_excess, XC_cntr, XC_cntr_old
1790 : integer :: c12, nz
1791 : include 'formats'
1792 10 : check_XC_cntr = keep_going
1793 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1794 10 : c12 = s% net_iso(ic12)
1795 10 : if (c12 == 0) return
1796 10 : nz = s% nz
1797 10 : XC_cntr = s% xa(c12,nz)
1798 10 : XC_cntr_old = s% xa_old(c12,nz)
1799 10 : if (s% delta_XC_drop_only .and. XC_cntr >= XC_cntr_old) return
1800 : check_XC_cntr = check_change(s, XC_cntr - XC_cntr_old, &
1801 : s% delta_XC_cntr_limit, s% delta_XC_cntr_hard_limit, &
1802 10 : nz, 'delta_XC_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1803 10 : if (check_XC_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1804 0 : write(*,1) 'XC_cntr', XC_cntr
1805 0 : write(*,1) 'XC_cntr_old', XC_cntr_old
1806 0 : write(*,1) 'delta', XC_cntr - XC_cntr_old
1807 : end if
1808 : end function check_XC_cntr
1809 :
1810 :
1811 10 : integer function check_XNe_cntr(s, skip_hard_limit, dt_limit_ratio)
1812 : type (star_info), pointer :: s
1813 : logical, intent(in) :: skip_hard_limit
1814 : real(dp), intent(inout) :: dt_limit_ratio
1815 : real(dp) :: relative_excess, XNe_cntr, XNe_cntr_old
1816 : integer :: ne20, nz
1817 : include 'formats'
1818 10 : check_XNe_cntr = keep_going
1819 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1820 10 : ne20 = s% net_iso(ine20)
1821 10 : if (ne20 == 0) return
1822 10 : nz = s% nz
1823 10 : XNe_cntr = s% xa(ne20,nz)
1824 10 : XNe_cntr_old = s% xa_old(ne20,nz)
1825 10 : if (s% delta_XNe_drop_only .and. XNe_cntr >= XNe_cntr_old) return
1826 : check_XNe_cntr = check_change(s, XNe_cntr - XNe_cntr_old, &
1827 : s% delta_XNe_cntr_limit, s% delta_XNe_cntr_hard_limit, &
1828 10 : nz, 'delta_XNe_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1829 10 : if (check_XNe_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1830 0 : write(*,1) 'XNe_cntr', XNe_cntr
1831 0 : write(*,1) 'XNe_cntr_old', XNe_cntr_old
1832 0 : write(*,1) 'delta', XNe_cntr - XNe_cntr_old
1833 : end if
1834 : end function check_XNe_cntr
1835 :
1836 :
1837 10 : integer function check_XO_cntr(s, skip_hard_limit, dt_limit_ratio)
1838 : type (star_info), pointer :: s
1839 : logical, intent(in) :: skip_hard_limit
1840 : real(dp), intent(inout) :: dt_limit_ratio
1841 : real(dp) :: relative_excess, XO_cntr, XO_cntr_old
1842 : integer :: o16, nz
1843 : include 'formats'
1844 10 : check_XO_cntr = keep_going
1845 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1846 10 : o16 = s% net_iso(io16)
1847 10 : if (o16 == 0) return
1848 10 : nz = s% nz
1849 10 : XO_cntr = s% xa(o16,nz)
1850 10 : XO_cntr_old = s% xa_old(o16,nz)
1851 10 : if (s% delta_XO_drop_only .and. XO_cntr >= XO_cntr_old) return
1852 : check_XO_cntr = check_change(s, XO_cntr - XO_cntr_old, &
1853 : s% delta_XO_cntr_limit, s% delta_XO_cntr_hard_limit, &
1854 10 : nz, 'delta_XO_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1855 10 : if (check_XO_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1856 0 : write(*,1) 'XO_cntr', XO_cntr
1857 0 : write(*,1) 'XO_cntr_old', XO_cntr_old
1858 0 : write(*,1) 'delta', XO_cntr - XO_cntr_old
1859 : end if
1860 : end function check_XO_cntr
1861 :
1862 :
1863 10 : integer function check_XSi_cntr(s, skip_hard_limit, dt_limit_ratio)
1864 : type (star_info), pointer :: s
1865 : logical, intent(in) :: skip_hard_limit
1866 : real(dp), intent(inout) :: dt_limit_ratio
1867 : real(dp) :: relative_excess, XSi_cntr, XSi_cntr_old
1868 : integer :: si28, nz
1869 : include 'formats'
1870 10 : check_XSi_cntr = keep_going
1871 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
1872 10 : si28 = s% net_iso(isi28)
1873 10 : if (si28 == 0) return
1874 0 : nz = s% nz
1875 0 : XSi_cntr = s% xa(si28,nz)
1876 0 : XSi_cntr_old = s% xa_old(si28,nz)
1877 0 : if (s% delta_XSi_drop_only .and. XSi_cntr >= XSi_cntr_old) return
1878 : check_XSi_cntr = check_change(s, XSi_cntr - XSi_cntr_old, &
1879 : s% delta_XSi_cntr_limit, s% delta_XSi_cntr_hard_limit, &
1880 0 : nz, 'delta_XSi_cntr', skip_hard_limit, dt_limit_ratio, relative_excess)
1881 0 : if (check_XSi_cntr /= keep_going .and. s% report_dt_hard_limit_retries) then
1882 0 : write(*,1) 'XSi_cntr', XSi_cntr
1883 0 : write(*,1) 'XSi_cntr_old', XSi_cntr_old
1884 0 : write(*,1) 'delta', XSi_cntr - XSi_cntr_old
1885 : end if
1886 : end function check_XSi_cntr
1887 :
1888 :
1889 10 : integer function check_delta_mdot(s, skip_hard_limit, dt, dt_limit_ratio)
1890 : type (star_info), pointer :: s
1891 : logical, intent(in) :: skip_hard_limit
1892 : real(dp), intent(in) :: dt
1893 : real(dp), intent(inout) :: dt_limit_ratio
1894 : real(dp) :: mdot, mdot_old, delta_mdot, lim, hard_lim
1895 10 : check_delta_mdot = keep_going
1896 10 : mdot = s% mstar_dot
1897 10 : mdot_old = s% mstar_dot_old
1898 : delta_mdot = abs(mdot - mdot_old)/ &
1899 : (s% delta_mdot_atol*Msun/secyer + &
1900 10 : s% delta_mdot_rtol*max(abs(mdot),abs(mdot_old)))
1901 10 : if (delta_mdot == 0) return
1902 0 : lim = s% delta_mdot_limit*s% time_delta_coeff
1903 0 : hard_lim = s% delta_mdot_hard_limit*s% time_delta_coeff
1904 0 : if (hard_lim > 0 .and. (.not. skip_hard_limit)) then
1905 0 : if (delta_mdot > hard_lim) then
1906 0 : if (s% report_dt_hard_limit_retries) &
1907 0 : write(*, '(a30, f20.10, 99e20.10)') 'delta_mdot_hard_limit', &
1908 0 : delta_mdot, hard_lim
1909 0 : s% retry_message = 'delta_mdot_hard_limit'
1910 0 : check_delta_mdot = retry
1911 0 : return
1912 : end if
1913 : end if
1914 0 : if (lim <= 0) return
1915 0 : dt_limit_ratio = delta_mdot/lim
1916 0 : if (dt_limit_ratio <= 1d0) dt_limit_ratio = 0
1917 : end function check_delta_mdot
1918 :
1919 :
1920 10 : integer function check_delta_mstar(s, skip_hard_limit, dt, dt_limit_ratio)
1921 : type (star_info), pointer :: s
1922 : logical, intent(in) :: skip_hard_limit
1923 : real(dp), intent(in) :: dt
1924 : real(dp), intent(inout) :: dt_limit_ratio
1925 : real(dp) :: delta_lg_star_mass, lim, hard_lim
1926 10 : check_delta_mstar = keep_going
1927 10 : delta_lg_star_mass = abs(log10(s% mstar/s% mstar_old))
1928 10 : lim = s% delta_lg_star_mass_limit*s% time_delta_coeff
1929 10 : hard_lim = s% delta_lg_star_mass_hard_limit*s% time_delta_coeff
1930 10 : if (hard_lim > 0 .and. (.not. skip_hard_limit)) then
1931 0 : if (delta_lg_star_mass > hard_lim) then
1932 0 : if (s% report_dt_hard_limit_retries) &
1933 0 : write(*, '(a30, f20.10, 99e20.10)') 'delta_lg_star_mass_hard_limit', &
1934 0 : delta_lg_star_mass, hard_lim
1935 0 : check_delta_mstar = retry
1936 0 : s% retry_message = 'delta_lg_star_mass_hard_limit'
1937 0 : return
1938 : end if
1939 : end if
1940 10 : if (lim <= 0) return
1941 10 : dt_limit_ratio = delta_lg_star_mass/lim
1942 10 : if (dt_limit_ratio <= 1d0) dt_limit_ratio = 0
1943 : end function check_delta_mstar
1944 :
1945 :
1946 10 : integer function check_adjust_J_q(s, skip_hard_limit, dt_limit_ratio)
1947 : type (star_info), pointer :: s
1948 : logical, intent(in) :: skip_hard_limit
1949 : real(dp), intent(inout) :: dt_limit_ratio
1950 : real(dp) :: relative_excess
1951 : include 'formats'
1952 10 : check_adjust_J_q = keep_going
1953 10 : dt_limit_ratio = 0d0
1954 10 : if (s% doing_relax) return
1955 10 : if (.not. (s% rotation_flag .and. s% do_adjust_J_lost .and. s% mstar_dot < 0d0)) return
1956 : ! we care about s% adjust_J_q remaining above a given limit
1957 : ! so we use 1-S% adjust_J_q
1958 : check_adjust_J_q = check_change(s, 1-s% adjust_J_q, &
1959 : 1-s% adjust_J_q_limit, &
1960 : 1-s% adjust_J_q_hard_limit, &
1961 : 0, 'check_adjust_J_q', &
1962 0 : .false., dt_limit_ratio, relative_excess)
1963 0 : end function check_adjust_J_q
1964 :
1965 :
1966 10 : integer function check_delta_lgL(s, skip_hard_limit, dt_limit_ratio)
1967 : type (star_info), pointer :: s
1968 : logical, intent(in) :: skip_hard_limit
1969 : real(dp), intent(inout) :: dt_limit_ratio
1970 : real(dp) :: dlgL, relative_excess
1971 : include 'formats'
1972 10 : check_delta_lgL = keep_going
1973 10 : dt_limit_ratio = 0d0
1974 10 : if (s% doing_relax) return
1975 10 : if (s% L_surf < s% delta_lgL_limit_L_min .or. s% L_surf_old <= 0d0) then
1976 0 : dlgL = 0
1977 : else
1978 10 : dlgL = log10(s% L_surf/s% L_surf_old)
1979 : end if
1980 10 : if (is_bad(dlgL)) then
1981 0 : write(*,2) 's% L_surf', s% model_number, s% L_surf
1982 0 : write(*,2) 's% L_surf_old', s% model_number, s% L_surf_old
1983 0 : call mesa_error(__FILE__,__LINE__,'check_delta_lgL')
1984 : end if
1985 : check_delta_lgL = check_change(s, dlgL, &
1986 : s% delta_lgL_limit, s% delta_lgL_hard_limit, &
1987 10 : 0, 'delta_lgL', skip_hard_limit, dt_limit_ratio, relative_excess)
1988 10 : if (check_delta_lgL /= keep_going .and. s% report_dt_hard_limit_retries) then
1989 0 : write(*,1) 'L_surf', s% L_surf
1990 0 : write(*,1) 'L_surf_old', s% L_surf_old
1991 : end if
1992 : end function check_delta_lgL
1993 :
1994 :
1995 10 : integer function check_delta_HR(s, skip_hard_limit, dt_limit_ratio)
1996 : type (star_info), pointer :: s
1997 : logical, intent(in) :: skip_hard_limit
1998 : real(dp), intent(inout) :: dt_limit_ratio
1999 : real(dp) :: dlgL, dlgTeff, dHR, relative_excess
2000 : include 'formats'
2001 10 : check_delta_HR = keep_going
2002 : if (s% L_phot_old <= 0 .or. s% Teff_old <= 0 .or. &
2003 10 : s% L_phot <= 0 .or. s% Teff <= 0) return
2004 10 : dlgL = log10(s% L_phot/s% L_phot_old)
2005 10 : dlgTeff = log10(s% Teff/s% Teff_old)
2006 10 : dHR = sqrt(pow2(s% delta_HR_ds_L*dlgL) + pow2(s% delta_HR_ds_Teff*dlgTeff))
2007 10 : if (is_bad(dHR)) then
2008 0 : write(*,1) 's% L_phot_old', s% L_phot_old
2009 0 : write(*,1) 's% L_phot', s% L_phot
2010 0 : write(*,1) 's% Teff_old', s% Teff_old
2011 0 : write(*,1) 's% Teff', s% Teff
2012 0 : write(*,1) 'dHR', dHR
2013 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'check_delta_HR')
2014 : end if
2015 : check_delta_HR = check_change(s, dHR, &
2016 : s% delta_HR_limit, s% delta_HR_hard_limit, &
2017 10 : 0, 'delta_HR', skip_hard_limit, dt_limit_ratio, relative_excess)
2018 10 : if (check_delta_HR /= keep_going .and. s% report_dt_hard_limit_retries) then
2019 0 : write(*,1) 's% L_phot_old', s% L_phot_old
2020 0 : write(*,1) 's% L_phot', s% L_phot
2021 0 : write(*,1) 's% Teff_old', s% Teff_old
2022 0 : write(*,1) 's% Teff', s% Teff
2023 0 : write(*,1) 'dHR', dHR
2024 : end if
2025 : end function check_delta_HR
2026 :
2027 :
2028 10 : integer function check_rel_error_in_energy(s, skip_hard_limit, dt_limit_ratio)
2029 : type (star_info), pointer :: s
2030 : logical, intent(in) :: skip_hard_limit
2031 : real(dp), intent(inout) :: dt_limit_ratio
2032 : real(dp) :: rel_error, relative_excess
2033 : include 'formats'
2034 10 : check_rel_error_in_energy = keep_going
2035 10 : dt_limit_ratio = 0d0
2036 10 : if (s% doing_relax) return
2037 10 : rel_error = abs(s% error_in_energy_conservation/s% total_energy_end)
2038 : check_rel_error_in_energy = check_change(s, rel_error, &
2039 : s% limit_for_rel_error_in_energy_conservation, &
2040 : s% hard_limit_for_rel_error_in_energy_conservation, &
2041 : 0, 'hard_limit_for_rel_error_in_energy_conservation', &
2042 10 : .false., dt_limit_ratio, relative_excess)
2043 10 : if (check_rel_error_in_energy /= keep_going .and. s% report_dt_hard_limit_retries) then
2044 0 : write(*,1) 'error_in_energy_conservation', s% error_in_energy_conservation
2045 0 : write(*,1) 'total_energy_end', s% total_energy_end
2046 0 : write(*,1) 'rel_error', s% error_in_energy_conservation/s% total_energy_end
2047 : end if
2048 : end function check_rel_error_in_energy
2049 :
2050 :
2051 10 : integer function check_dt_div_dt_cell_collapse(s, skip_hard_limit, dt, dt_limit_ratio)
2052 : use star_utils, only: eval_min_cell_collapse_time
2053 : type (star_info), pointer :: s
2054 : logical, intent(in) :: skip_hard_limit
2055 : real(dp), intent(in) :: dt
2056 : real(dp), intent(inout) :: dt_limit_ratio
2057 : real(dp) :: ratio, dt_timescale, relative_excess
2058 : integer :: min_collapse_k, ierr
2059 : include 'formats'
2060 10 : check_dt_div_dt_cell_collapse = keep_going
2061 10 : if (s% doing_relax) return
2062 : dt_timescale = eval_min_cell_collapse_time( &
2063 10 : s, 2, s% nz, min_collapse_k, ierr)
2064 10 : if (ierr /= 0) return
2065 10 : if (dt_timescale < 1d-30) return
2066 10 : ratio = dt/dt_timescale
2067 : !write(*,2) 'dt dt_cell_collapse ratio', min_collapse_k, dt, dt_timescale, ratio
2068 : check_dt_div_dt_cell_collapse = check_change(s, ratio, &
2069 : s% dt_div_dt_cell_collapse_limit, s% dt_div_dt_cell_collapse_hard_limit, &
2070 10 : min_collapse_k, 'dt_div_dt_cell_collapse', skip_hard_limit, dt_limit_ratio, relative_excess)
2071 10 : if (check_dt_div_dt_cell_collapse /= keep_going .and. s% report_dt_hard_limit_retries) then
2072 0 : write(*,2) 'min dt_cell_collapse', min_collapse_k, dt_timescale
2073 : end if
2074 : end function check_dt_div_dt_cell_collapse
2075 :
2076 :
2077 10 : integer function check_dt_div_min_dr_div_cs(s, skip_hard_limit, dt, dt_limit_ratio)
2078 : use star_utils, only: min_dr_div_cs
2079 : type (star_info), pointer :: s
2080 : logical, intent(in) :: skip_hard_limit
2081 : real(dp), intent(in) :: dt
2082 : real(dp), intent(inout) :: dt_limit_ratio
2083 : real(dp) :: ratio, dt_x, relative_excess
2084 : include 'formats'
2085 10 : check_dt_div_min_dr_div_cs = keep_going
2086 10 : dt_limit_ratio = 0d0
2087 10 : if (s% doing_relax) return
2088 10 : if (s% dt_div_min_dr_div_cs_limit <= 0d0) return
2089 0 : dt_x = min_dr_div_cs(s, s% Tlim_dt_div_min_dr_div_cs_cell)
2090 : !write(*,2) 'log min_dr_div_cs, q, m', s% Tlim_dt_div_min_dr_div_cs_cell, &
2091 : ! safe_log10(dt_x), s% q(s% Tlim_dt_div_min_dr_div_cs_cell), s% m(s% Tlim_dt_div_min_dr_div_cs_cell)/Msun
2092 0 : ratio = dt/dt_x
2093 : !write(*,3) 'dt/dx_x log', s% Tlim_dt_div_min_dr_div_cs_cell, s% model_number, ratio, safe_log10(ratio)
2094 : check_dt_div_min_dr_div_cs = check_change(s, ratio, &
2095 : s% dt_div_min_dr_div_cs_limit, s% dt_div_min_dr_div_cs_hard_limit, &
2096 : s% Tlim_dt_div_min_dr_div_cs_cell, 'dt_div_min_dr_div_cs', &
2097 0 : skip_hard_limit, dt_limit_ratio, relative_excess)
2098 0 : if ((check_dt_div_min_dr_div_cs /= keep_going .and. s% report_dt_hard_limit_retries) .or. &
2099 : (ratio > 1d0 .and. s% report_min_dr_div_cs)) then
2100 0 : write(*,2) 'min_dr_div_cs', s% Tlim_dt_div_min_dr_div_cs_cell, dt_x
2101 : end if
2102 : end function check_dt_div_min_dr_div_cs
2103 :
2104 :
2105 10 : integer function check_dX_nuc_drop(s, skip_hard_limit, dt, dt_limit_ratio)
2106 : use rates_def, only: i_rate
2107 : type (star_info), pointer :: s
2108 : logical, intent(in) :: skip_hard_limit
2109 : real(dp), intent(in) :: dt
2110 : real(dp), intent(inout) :: dt_limit_ratio
2111 :
2112 : integer :: k, nz, max_k, max_j
2113 10 : real(dp), pointer, dimension(:) :: sig
2114 : real(dp) :: max_dx_nuc_drop, X_limit, A_limit, min_dt, limit, hard_limit
2115 :
2116 : logical, parameter :: dbg = .false.
2117 :
2118 : include 'formats'
2119 :
2120 10 : check_dX_nuc_drop = keep_going
2121 10 : if (s% mix_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) return
2122 :
2123 10 : X_limit = s% dX_nuc_drop_min_X_limit
2124 10 : A_limit = s% dX_nuc_drop_max_A_limit
2125 10 : nz = s% nz
2126 10 : sig => s% sig
2127 :
2128 10 : max_dx_nuc_drop = 0
2129 10 : max_k = 0
2130 10 : max_j = 0
2131 10633 : do k = 1, nz
2132 10633 : call do1(k)
2133 : end do
2134 :
2135 10 : s% Tlim_dXnuc_drop_cell = max_k
2136 10 : s% Tlim_dXnuc_drop_species = max_j
2137 :
2138 10 : hard_limit = s% dX_nuc_drop_hard_limit*s% time_delta_coeff
2139 10 : if (hard_limit > 0 .and. (.not. skip_hard_limit) .and. &
2140 : max_dx_nuc_drop > hard_limit) then
2141 0 : if (s% report_dt_hard_limit_retries) then
2142 0 : write(*,2) trim(chem_isos% name(s% chem_id(max_j))), max_k, s% xa(max_j,max_k)
2143 0 : write(*,2) trim(chem_isos% name(s% chem_id(max_j))) // ' old', max_k, s% xa_old(max_j,max_k)
2144 0 : write(*,2) 'drop', max_k, max_dx_nuc_drop
2145 : end if
2146 0 : s% retry_message = 'dX_nuc_drop_hard_limit'
2147 0 : s% retry_message_k = max_k
2148 0 : check_dX_nuc_drop = retry
2149 0 : return
2150 : end if
2151 :
2152 10 : limit = s% dX_nuc_drop_limit*s% time_delta_coeff
2153 10 : if (s% log_max_temperature >= 9.45d0 .and. s% dX_nuc_drop_limit_at_high_T > 0) &
2154 10 : limit = s% dX_nuc_drop_limit_at_high_T
2155 10 : if (limit <= 0 .or. max_dx_nuc_drop <= 0) return
2156 :
2157 10 : if (dt < secyer*s% dX_nuc_drop_min_yrs_for_dt) return
2158 10 : min_dt = secyer*s% dX_nuc_drop_min_yrs_for_dt
2159 : dt_limit_ratio = min( &
2160 : max_dx_nuc_drop/limit, &
2161 10 : 1d0*dt/min_dt)
2162 10 : if (dt_limit_ratio <= 1d0) dt_limit_ratio = 0
2163 :
2164 10 : s% dX_nuc_drop_max_j = max_j
2165 10 : s% dX_nuc_drop_max_k = max_k
2166 :
2167 : contains
2168 :
2169 10623 : subroutine do1(k)
2170 : integer, intent(in) :: k
2171 :
2172 : integer :: j
2173 : real(dp) :: dx, dx_drop, dm, dt_dm, dx_burning, dx_inflow, dxdt_nuc
2174 : real(dp) :: dx00, dxp1, sig00, sigp1, flux00, fluxp1
2175 :
2176 : include 'formats'
2177 :
2178 10623 : dm = s% dm(k)
2179 10623 : dt_dm = dt/dm
2180 :
2181 95607 : do j=1, s% species
2182 :
2183 84984 : if (chem_isos% W(s% chem_id(j)) > A_limit) then
2184 : if (dbg .and. k == 1387) &
2185 : write(*,2) 'dX_nuc_max_A_limit ' // trim(chem_isos% name(s% chem_id(j))), &
2186 : k, s% xa(j,k), chem_isos% W(s% chem_id(j)), A_limit
2187 : cycle
2188 : end if
2189 :
2190 84984 : if (chem_isos% Z(s% chem_id(j)) <= 2) then
2191 : cycle ! skip the little guys
2192 : end if
2193 :
2194 53115 : if (s% xa(j,k) < X_limit) then
2195 : if (dbg .and. k == 1387) &
2196 : write(*,2) &
2197 : 'dX_nuc_drop_min_X_limit ' // trim(chem_isos% name(s% chem_id(j))), &
2198 : k, s% xa(j,k), X_limit
2199 : cycle
2200 : end if
2201 :
2202 51836 : dxdt_nuc = s% dxdt_nuc(j,k)
2203 51836 : if (dxdt_nuc >= 0) cycle
2204 :
2205 7651 : dx_burning = dxdt_nuc*dt
2206 7651 : sig00 = sig(k)
2207 :
2208 7651 : if (k < s% nz) then
2209 7631 : sigp1 = sig(k+1)
2210 : else
2211 : sigp1 = 0
2212 : end if
2213 :
2214 7651 : if (k > 1) then
2215 7651 : dx00 = s% xa(j,k-1) - s% xa(j,k)
2216 7651 : flux00 = -sig00*dx00
2217 : else
2218 : flux00 = 0
2219 : end if
2220 :
2221 7651 : if (k < s% nz) then
2222 7631 : dxp1 = s% xa(j,k) - s% xa(j,k+1)
2223 7631 : fluxp1 = -sigp1*dxp1
2224 : else
2225 : fluxp1 = 0
2226 : end if
2227 :
2228 7651 : dx_inflow = max(0d0, fluxp1, -flux00)*dt_dm
2229 :
2230 7651 : dx_drop = -(dx_burning + dx_inflow) ! dx_burning < 0 for drop
2231 :
2232 7651 : dx = s% xa_old(j,k) - s% xa(j,k) ! the actual drop
2233 : if (dx < dx_drop) dx_drop = dx
2234 :
2235 18274 : if (dx_drop > max_dx_nuc_drop) then
2236 2607 : max_dx_nuc_drop = dx_drop
2237 2607 : max_k = k
2238 2607 : max_j = j
2239 : end if
2240 :
2241 : end do
2242 :
2243 10623 : end subroutine do1
2244 :
2245 :
2246 : end function check_dX_nuc_drop
2247 :
2248 :
2249 11 : integer function check_varcontrol_limit(s, dt_limit_ratio)
2250 : type (star_info), pointer :: s
2251 : real(dp), intent(inout) :: dt_limit_ratio
2252 :
2253 : real(dp) :: varcontrol, vc_target
2254 : integer :: ierr
2255 :
2256 : include 'formats'
2257 :
2258 : ierr = 0
2259 11 : check_varcontrol_limit = keep_going
2260 :
2261 11 : varcontrol = eval_varcontrol(s, ierr)
2262 11 : if (ierr /= 0) then
2263 0 : check_varcontrol_limit = retry
2264 0 : s% retry_message = 'varcontrol hard limit'
2265 0 : if (s% report_ierr) write(*, *) 'check_varcontrol_limit: eval_varcontrol ierr', ierr
2266 0 : s% result_reason = nonzero_ierr
2267 0 : return
2268 : end if
2269 :
2270 11 : if (s% varcontrol_target < s% min_allowed_varcontrol_target) then
2271 0 : check_varcontrol_limit = terminate
2272 0 : write(*, *) 'ERROR: timestep varcontrol_target < min_allowed_varcontrol_target'
2273 0 : s% result_reason = nonzero_ierr
2274 0 : return
2275 : end if
2276 :
2277 11 : vc_target = max(1d-99,s% varcontrol_target)
2278 11 : dt_limit_ratio = varcontrol/vc_target
2279 :
2280 11 : if (s% report_solver_dt_info) then
2281 0 : write(*, 1) 's% varcontrol_target', s% varcontrol_target
2282 0 : write(*, 1) 'vc_target', vc_target
2283 0 : write(*, 1) 'varcontrol', varcontrol
2284 0 : write(*, 1) 'dt_limit_ratio', dt_limit_ratio
2285 0 : write(*, *)
2286 : end if
2287 :
2288 11 : if (dt_limit_ratio > s% varcontrol_dt_limit_ratio_hard_max) then
2289 0 : write(*, '(a50, f20.10, 99e20.10)') 'varcontrol_dt_limit_ratio too large', &
2290 0 : dt_limit_ratio, varcontrol, vc_target
2291 0 : check_varcontrol_limit = retry
2292 0 : s% retry_message = 'varcontrol_dt_limit_ratio_hard_max'
2293 0 : return
2294 : end if
2295 :
2296 : end function check_varcontrol_limit
2297 :
2298 :
2299 11 : real(dp) function eval_varcontrol(s, ierr) result(varcontrol)
2300 : type (star_info), pointer :: s
2301 : integer, intent(out) :: ierr
2302 :
2303 : integer :: j, nterms, nvar_hydro, nz, k
2304 22 : real(dp) :: sumj, sumvar, sumscales, sumterm(s% nvar_total)
2305 : real(dp), parameter :: xscale_min = 1
2306 :
2307 : include 'formats'
2308 11 : ierr = 0
2309 :
2310 11 : varcontrol = 1d99
2311 11 : nvar_hydro = s% nvar_hydro
2312 11 : nz = s% nz
2313 :
2314 11 : nterms = 0
2315 11 : sumvar = 0
2316 11 : sumscales = 0
2317 143 : sumterm(:) = 0
2318 :
2319 11 : if (.not. associated(s% xh_old)) then
2320 0 : call mesa_error(__FILE__,__LINE__,'not associated xh_old')
2321 : end if
2322 :
2323 : ! use differences in smoothed old and new to filter out high frequency noise.
2324 55 : do j = 1, nvar_hydro
2325 :
2326 : if (j == s% i_lum .or. &
2327 : j == s% i_u .or. &
2328 : j == s% i_v .or. &
2329 : j == s% i_w .or. &
2330 : j == s% i_Hp .or. &
2331 : j == s% i_j_rot .or. &
2332 : j == s% i_w_div_wc .or. &
2333 44 : j == s% i_alpha_RTI .or. &
2334 : j == s% i_Et_RSP) cycle
2335 :
2336 33 : nterms = nterms + nz
2337 39162 : do k = 3, nz-2
2338 430419 : sumj = abs(sum(s% xh(j,k-2:k+2)) - sum(s% xh_old(j,k-2:k+2)))/5
2339 39162 : sumterm(j) = sumterm(j) + sumj
2340 : end do
2341 : sumterm(j) = sumterm(j) + &
2342 : abs((2*s% xh(j,1) + s% xh(j,2)) - (2*s% xh_old(j,1) + s% xh_old(j,2)))/3 + &
2343 33 : abs((2*s% xh(j,nz) + s% xh(j,nz-1)) - (2*s% xh_old(j,nz) + s% xh_old(j,nz-1)))/3
2344 33 : k = 2
2345 231 : sumj = abs(sum(s% xh(j,k-1:k+1)) - sum(s% xh_old(j,k-1:k+1)))/3
2346 33 : sumterm(j) = sumterm(j) + sumj
2347 33 : k = nz-1
2348 33 : sumj = abs(sum(s% xh(j,k-1:k+1)) - sum(s% xh_old(j,k-1:k+1)))/3
2349 :
2350 33 : if (j == s% i_lnd) then
2351 11 : sumterm(j) = sumterm(j)/3 ! Seems to help. from Eggleton.
2352 : end if
2353 :
2354 33 : sumvar = sumvar + sumterm(j)
2355 55 : sumscales = sumscales + max(xscale_min, abs(s% xh_old(j,1)))
2356 :
2357 : end do
2358 :
2359 : sumterm(:) = sumterm(:)/sumscales
2360 11 : sumvar = sumvar/sumscales
2361 :
2362 11 : varcontrol = sumvar/nterms
2363 :
2364 11 : end function eval_varcontrol
2365 :
2366 :
2367 11 : subroutine filter_dt_next(s, order, dt_limit_ratio_in)
2368 : ! H211b "low pass" controller.
2369 : ! Soderlind & Wang, J of Computational and Applied Math 185 (2006) 225 – 243.
2370 : use num_def
2371 : type (star_info), pointer :: s
2372 : real(dp), intent(in) :: order, dt_limit_ratio_in
2373 :
2374 : real(dp) :: ratio, ratio_prev, limtr, dt_limit_ratio_target, &
2375 : dt_limit_ratio, beta1, beta2, alpha2
2376 :
2377 : include 'formats'
2378 :
2379 11 : beta1 = 0.25d0/order
2380 11 : beta2 = 0.25d0/order
2381 11 : alpha2 = 0.25d0
2382 :
2383 11 : dt_limit_ratio = max(1d-10, dt_limit_ratio_in)
2384 11 : s% dt_limit_ratio = dt_limit_ratio
2385 11 : dt_limit_ratio_target = 1d0
2386 :
2387 : if (s% use_dt_low_pass_controller .and. &
2388 11 : s% dt_limit_ratio_old > 0 .and. s% dt_old > 0) then ! use 2 values to do "low pass" controller
2389 10 : ratio = limiter(dt_limit_ratio_target/dt_limit_ratio)
2390 10 : ratio_prev = limiter(dt_limit_ratio_target/s% dt_limit_ratio_old)
2391 : limtr = limiter( &
2392 10 : pow(ratio,beta1) * pow(ratio_prev,beta2) * pow(s% dt/s% dt_old,-alpha2))
2393 10 : s% dt_next = s% dt*limtr
2394 :
2395 10 : if (s% report_solver_dt_info) then
2396 0 : write(*,2) 'dt_limit_ratio_target', s% model_number, dt_limit_ratio_target
2397 0 : write(*,2) 'dt_limit_ratio', s% model_number, dt_limit_ratio
2398 0 : write(*,2) 's% dt_limit_ratio_old', s% model_number, s% dt_limit_ratio_old
2399 0 : write(*,2) 'order', s% model_number, order
2400 0 : write(*,2) 'ratio', s% model_number, ratio
2401 0 : write(*,2) 'ratio_prev', s% model_number, ratio_prev
2402 0 : write(*,2) 'limtr', s% model_number, limtr
2403 0 : write(*,2) 's% dt_next', s% model_number, s% dt_next
2404 0 : write(*,'(A)')
2405 : end if
2406 :
2407 : else ! no history available, so fall back to the 1st order controller
2408 1 : s% dt_next = s% dt*dt_limit_ratio_target/dt_limit_ratio
2409 1 : if (s% report_solver_dt_info) then
2410 0 : write(*,2) 's% dt', s% model_number, s% dt
2411 0 : write(*,2) 'dt_limit_ratio_target', s% model_number, dt_limit_ratio_target
2412 0 : write(*,2) 'dt_limit_ratio', s% model_number, dt_limit_ratio
2413 0 : write(*,2) 'filter_dt_next', s% model_number, s% dt_next
2414 : end if
2415 1 : if (s% dt_next == 0d0) call mesa_error(__FILE__,__LINE__,'filter_dt_next')
2416 : end if
2417 :
2418 :
2419 : contains
2420 :
2421 :
2422 30 : real(dp) function limiter(x)
2423 : real(dp), intent(in) :: x
2424 : real(dp), parameter :: kappa = 10 ! 2
2425 : ! for x >= 0 and kappa = 2, limiter value is between 0.07 and 4.14
2426 : ! for x >= 0 and kappa = 10, limiter value is between 0.003 and 16.7
2427 : ! for x = 1, limiter = 1
2428 30 : limiter = 1 + kappa*atan((x-1)/kappa)
2429 30 : end function limiter
2430 :
2431 :
2432 : end subroutine filter_dt_next
2433 :
2434 :
2435 : end module timestep
2436 :
2437 :
|
