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