Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 Pablo Marchant & 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 :
21 : module binary_timestep
22 :
23 : use const_def, only: dp, msun, secyer
24 : use math_lib
25 : use star_lib
26 : use star_def
27 : use binary_def
28 :
29 : implicit none
30 :
31 : contains
32 :
33 0 : subroutine set_star_timesteps(b) ! sets the smallest next timestep for all stars
34 : type (binary_info), pointer :: b
35 : integer :: i, l
36 0 : real(dp) :: dt_min
37 : type (star_info), pointer :: s
38 : integer :: ierr, num_stars
39 0 : ierr = 0
40 0 : dt_min = 1d99
41 0 : num_stars = 2
42 0 : if (b% point_mass_i > 0) num_stars = 1
43 0 : do l = 1, num_stars
44 0 : if (l == 1 .and. b% point_mass_i == 1) then
45 : i = 2
46 : else
47 0 : i = l
48 : end if
49 0 : call star_ptr(b% star_ids(i), s, ierr)
50 0 : if (ierr /= 0) then
51 0 : write(*, *) trim('star_ptr') // ' ierr', ierr
52 0 : return
53 : end if
54 0 : if (s% dt_next < dt_min) then
55 0 : dt_min = s% dt_next
56 : end if
57 : end do
58 0 : if (b% max_timestep <= dt_min) then
59 : dt_min = b% max_timestep
60 : else
61 0 : b% dt_why_reason = b_Tlim_comp
62 : end if
63 : ! just to be sure we dont cause a segfault
64 0 : if (b% dt_why_reason < 1 .or. b% dt_why_reason > b_numTlim) then
65 0 : dt_why_str(Tlim_binary) = " "
66 : else
67 0 : dt_why_str(Tlim_binary) = binary_dt_why_str(b% dt_why_reason)
68 : end if
69 0 : do l = 1, num_stars
70 0 : if (l == 1 .and. b% point_mass_i == 1) then
71 : i = 2
72 : else
73 0 : i = l
74 : end if
75 0 : call star_ptr(b% star_ids(i), s, ierr)
76 0 : if (ierr /= 0) then
77 0 : write(*, *) trim('star_ptr') // ' ierr', ierr
78 0 : return
79 : end if
80 0 : if (s% dt_next > dt_min) then
81 0 : s% dt_next = dt_min
82 0 : s% why_Tlim = Tlim_binary
83 : end if
84 : end do
85 :
86 0 : if (b% have_to_reduce_timestep_due_to_j) then
87 : ! lower timesteps after retries due to large changes in angular momentum
88 0 : if (b% point_mass_i /= 1) then
89 0 : b% s1% dt = b% s1% dt*b% dt_reduction_factor_for_j
90 : end if
91 0 : if (b% point_mass_i /= 2) then
92 0 : b% s2% dt = b% s2% dt*b% dt_reduction_factor_for_j
93 : end if
94 0 : b% have_to_reduce_timestep_due_to_j = .false.
95 : end if
96 :
97 : end subroutine set_star_timesteps
98 :
99 2 : integer function binary_pick_next_timestep(b)
100 : type (binary_info), pointer :: b
101 : type (star_info), pointer :: s
102 :
103 : real(dp) :: &
104 2 : env_change, dtm, dtj, dta, dtr, dte, dtdm, &
105 2 : j_change, sep_change, rel_gap_change, e_change, set_dt, &
106 2 : rel_change
107 :
108 : include 'formats'
109 :
110 2 : dtm = 1d99
111 2 : dtj = 1d99
112 2 : dta = 1d99
113 2 : dtr = 1d99
114 2 : dte = 1d99
115 2 : dtdm = 1d99
116 :
117 2 : binary_pick_next_timestep = keep_going
118 :
119 2 : s => b% s_donor
120 :
121 2 : if (b% max_timestep < 0) b% max_timestep = b% s_donor% dt
122 :
123 2 : b% env(b% d_i) = s% star_mass - s% he_core_mass
124 2 : if (b% point_mass_i == 0) then
125 0 : b% env(b% a_i) = b% s_accretor% star_mass - b% s_accretor% he_core_mass
126 : end if
127 :
128 2 : if(.not. b% doing_first_model_of_run) then
129 2 : if (b% env_old(b% d_i) /= 0) then
130 0 : env_change = b% env(b% d_i) - b% env_old(b% d_i)
131 : else
132 : env_change = 0
133 : end if
134 :
135 2 : if (b% rl_relative_gap_old(b% d_i) /= 0) then
136 2 : rel_gap_change = b% rl_relative_gap_old(b% d_i) - b% rl_relative_gap(b% d_i)
137 : else
138 : rel_gap_change = 0
139 : end if
140 :
141 2 : if (b% angular_momentum_j_old /= 0) then
142 0 : j_change = b% angular_momentum_j - b% angular_momentum_j_old
143 : else
144 : j_change = 0
145 : end if
146 :
147 2 : if (b% separation_old /= 0) then
148 0 : sep_change = b% separation - b% separation_old
149 : else
150 : sep_change = 0
151 : end if
152 2 : if (b% eccentricity_old /= 0) then
153 0 : e_change = b% eccentricity - b% eccentricity_old
154 : else
155 : e_change = 0
156 : end if
157 : else
158 : env_change = 0
159 : rel_gap_change = 0
160 : j_change = 0
161 : sep_change = 0
162 : e_change = 0
163 : end if
164 :
165 : ! get limits for dt based on relative changes
166 2 : if (b% fj > 0) then
167 0 : rel_change = abs(j_change/b% angular_momentum_j)
168 : if (.not. b% ignore_hard_limits_this_step .and. &
169 0 : b% fj_hard > 0d0 .and. rel_change > b% fj_hard * b% time_delta_coeff) then
170 0 : write(*,*) "retry because of fj_hard limit,", &
171 0 : "fj_hard:", b% fj_hard, "rel_change:", rel_change
172 0 : binary_pick_next_timestep = retry
173 0 : b% have_to_reduce_timestep_due_to_j = .true.
174 0 : return
175 : end if
176 0 : dtj = s% time_step/(rel_change/(b% fj * b% time_delta_coeff)+1d-99)
177 : end if
178 :
179 2 : if (b% fm > 0) then
180 0 : rel_change = abs(env_change/max(b% env(b% d_i), b% fm_limit))
181 : if (.not. b% ignore_hard_limits_this_step .and. &
182 0 : b% fm_hard > 0d0 .and. rel_change > b% fm_hard * b% time_delta_coeff) then
183 0 : write(*,*) "retry because of fm_hard limit,", &
184 0 : "fm_hard:", b% fm_hard, "rel_change:", rel_change
185 0 : binary_pick_next_timestep = retry
186 0 : return
187 : end if
188 0 : dtm = s% time_step/(rel_change/(b% fm * b% time_delta_coeff)+1d-99)
189 : end if
190 :
191 2 : if (b% fr > 0) then
192 2 : rel_change = abs(rel_gap_change/max(abs(b% rl_relative_gap(b% d_i)), b% fr_limit))
193 : if (.not. b% ignore_hard_limits_this_step .and. &
194 2 : b% fr_hard > 0d0 .and. rel_change > b% fr_hard * b% time_delta_coeff) then
195 0 : write(*,*) "retry because of fr_hard limit for donor,", &
196 0 : "fr_hard:", b% fr_hard, "rel_change:", rel_change
197 0 : binary_pick_next_timestep = retry
198 0 : return
199 : end if
200 2 : dtr = s% time_step/(rel_change/(b% fr * b% time_delta_coeff)+1d-99)
201 :
202 : ! Check for accretor as well
203 2 : if (b% rl_relative_gap_old(b% a_i) /= 0) then
204 0 : rel_gap_change = b% rl_relative_gap_old(b% a_i) - b% rl_relative_gap(b% a_i)
205 : else
206 : rel_gap_change = 0
207 : end if
208 2 : rel_change = abs(rel_gap_change/max(abs(b% rl_relative_gap(b% a_i)), b% fr_limit))
209 : if (.not. b% ignore_hard_limits_this_step .and. &
210 2 : b% fr_hard > 0d0 .and. rel_change > b% fr_hard * b% time_delta_coeff) then
211 0 : write(*,*) "retry because of fr_hard limit for accretor,", &
212 0 : "fr_hard:", b% fr_hard, "rel_change:", rel_change
213 0 : binary_pick_next_timestep = retry
214 0 : return
215 : end if
216 2 : dtr = min(dtr, s% time_step/(rel_change/(b% fr * b% time_delta_coeff)+1d-99))
217 : end if
218 2 : if (dtr < b% fr_dt_limit) dtr = b% fr_dt_limit
219 :
220 2 : if (b% fa > 0) then
221 0 : rel_change = abs(sep_change/b% separation)
222 : if (.not. b% ignore_hard_limits_this_step .and. &
223 0 : b% fa_hard > 0d0 .and. rel_change > b% fa_hard * b% time_delta_coeff) then
224 0 : write(*,*) "retry because of fa_hard limit,", &
225 0 : "fa_hard:", b% fa_hard, "rel_change:", rel_change
226 0 : binary_pick_next_timestep = retry
227 0 : return
228 : end if
229 0 : dta = s% time_step/(rel_change/(b% fa * b% time_delta_coeff)+1d-99)
230 : end if
231 :
232 2 : if (b% fe > 0) then
233 0 : rel_change = abs(e_change/ max(b% eccentricity, b% fe_limit))
234 : if (.not. b% ignore_hard_limits_this_step .and. &
235 0 : b% fe_hard > 0d0 .and. rel_change > b% fe_hard * b% time_delta_coeff) then
236 0 : write(*,*) "retry because of fe_hard limit,", &
237 0 : "fe_hard:", b% fe_hard, "rel_change:", rel_change
238 0 : binary_pick_next_timestep = retry
239 0 : return
240 : end if
241 0 : dte = s% time_step/(rel_change/(b% fe * b% time_delta_coeff)+1d-99)
242 : end if
243 :
244 2 : if (b% fdm > 0d0) then
245 0 : rel_change = abs(b% m(b% d_i) - b% m_old(b% d_i))/b% m_old(b% d_i)
246 : if (.not. b% ignore_hard_limits_this_step .and. &
247 0 : b% fdm_hard > 0d0 .and. rel_change > b% fdm_hard * b% time_delta_coeff) then
248 0 : write(*,*) "retry because of fdm_hard limit for donor,", &
249 0 : "fdm_hard:", b% fdm_hard, "rel_change:", rel_change
250 0 : binary_pick_next_timestep = retry
251 0 : return
252 : end if
253 0 : dtdm = s% time_step/(rel_change/(b% fdm * b% time_delta_coeff)+1d-99)
254 :
255 0 : rel_change = abs(b% m(b% a_i) - b% m_old(b% a_i))/b% m_old(b% a_i)
256 : if (.not. b% ignore_hard_limits_this_step .and. &
257 0 : b% fdm_hard > 0d0 .and. rel_change > b% fdm_hard * b% time_delta_coeff) then
258 0 : write(*,*) "retry because of fdm_hard limit for accretor,", &
259 0 : "fdm_hard:", b% fdm_hard, "rel_change:", rel_change
260 0 : binary_pick_next_timestep = retry
261 0 : return
262 : end if
263 0 : dtdm = min(dtdm, s% time_step/(rel_change/(b% fdm * b% time_delta_coeff)+1d-99))
264 : end if
265 :
266 2 : set_dt = min(dtm, dtr, dtj, dta, dte, dtdm)
267 :
268 2 : if (set_dt == dtm) then
269 0 : b% dt_why_reason = b_Tlim_env
270 2 : else if (set_dt == dtr) then
271 2 : b% dt_why_reason = b_Tlim_roche
272 0 : else if (set_dt == dtj) then
273 0 : b% dt_why_reason = b_Tlim_jorb
274 0 : else if (set_dt == dta) then
275 0 : b% dt_why_reason = b_Tlim_sep
276 0 : else if (set_dt == dte) then
277 0 : b% dt_why_reason = b_Tlim_ecc
278 0 : else if (set_dt == dtdm) then
279 0 : b% dt_why_reason = b_Tlim_dm
280 : else
281 0 : call mesa_error(__FILE__,__LINE__,'Something wrong in binary timestep')
282 : end if
283 :
284 2 : if (set_dt < 1d-7) set_dt = 1d-7 ! there's a limit to everything
285 :
286 : b% max_timestep = exp10(b% dt_softening_factor*log10(b% max_timestep) + &
287 2 : (1-b% dt_softening_factor)*log10(set_dt*secyer))
288 :
289 : ! use variable varcontrols for different phases of evolution
290 2 : if (abs(b% mtransfer_rate)/Msun*secyer > 1d-20) then
291 0 : if (b% s_donor% center_h1 > 1d-12 .and. b% varcontrol_case_a > 0d0) then
292 0 : b% s_donor% varcontrol_target = b% varcontrol_case_a
293 0 : if (b% point_mass_i == 0) &
294 0 : b% s_accretor% varcontrol_target = b% varcontrol_case_a
295 0 : else if (b% s_donor% center_h1 < 1d-12 .and. b% varcontrol_case_b > 0d0) then
296 0 : b% s_donor% varcontrol_target = b% varcontrol_case_b
297 0 : if (b% point_mass_i == 0) &
298 0 : b% s_accretor% varcontrol_target = b% varcontrol_case_b
299 : end if
300 : else
301 2 : if (b% s_donor% center_h1 > 1d-12) then
302 0 : if (b% varcontrol_ms > 0d0) &
303 0 : b% s_donor% varcontrol_target = b% varcontrol_ms
304 : else
305 2 : if (b% varcontrol_post_ms > 0d0) &
306 0 : b% s_donor% varcontrol_target = b% varcontrol_post_ms
307 : end if
308 :
309 2 : if (b% point_mass_i == 0) then
310 0 : if (b% s_accretor% center_h1 > 1d-12) then
311 0 : if (b% varcontrol_ms > 0d0) &
312 0 : b% s_accretor% varcontrol_target = b% varcontrol_ms
313 : else
314 0 : if (b% varcontrol_post_ms > 0d0) &
315 0 : b% s_accretor% varcontrol_target = b% varcontrol_post_ms
316 : end if
317 : end if
318 : end if
319 :
320 2 : b% ignore_hard_limits_this_step = .false.
321 :
322 2 : end function binary_pick_next_timestep
323 :
324 :
325 : end module binary_timestep
|
