Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2013-2022 Pablo Marchant, Matthias Fabry & 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 run_binary_support
21 :
22 : use star_lib
23 : use star_def
24 : use const_def, only: dp, i8, secday
25 : use utils_lib
26 : use binary_def
27 : use binary_private_def
28 : use binary_ctrls_io, only : do_one_binary_setup
29 : use binary_do_one_utils
30 : use binary_ce
31 : use binary_photos
32 :
33 : implicit none
34 :
35 : contains
36 :
37 0 : subroutine do_run1_binary(tst, &
38 : ! star extras
39 : extras_controls, &
40 : ! binary extras
41 : extras_binary_controls, &
42 : ierr, &
43 : inlist_fname_arg)
44 :
45 : use binary_job_ctrls_io
46 : use binary_mdot, only : adjust_mdots, set_accretion_composition
47 : use binary_tides, only : sync_spin_orbit_torque
48 : use binary_evolve
49 : use mod_other_rlo_mdot
50 : use mod_other_implicit_rlo
51 : use mod_other_tsync
52 : use mod_other_sync_spin_to_orbit
53 : use mod_other_mdot_edd
54 : use mod_other_adjust_mdots
55 : use mod_other_accreted_material_j
56 : use mod_other_binary_jdot
57 : use mod_other_binary_wind_transfer
58 : use mod_other_binary_edot
59 : use mod_other_binary_ce
60 : use mod_other_binary_extras
61 : use mod_other_binary_photos
62 : use mod_other_e2
63 : use mod_other_pgbinary_plots
64 : use binary_timestep
65 : use binary_history
66 : use binary_history_specs
67 : use run_star_support
68 : use pgbinary, only : read_pgbinary_inlist, update_pgbinary_plots, &
69 : start_new_run_for_pgbinary, restart_run_for_pgbinary
70 :
71 : logical, intent(in) :: tst
72 :
73 : interface
74 :
75 : subroutine extras_controls(id, ierr)
76 : implicit none
77 : integer, intent(in) :: id
78 : integer, intent(out) :: ierr
79 : end subroutine extras_controls
80 :
81 : subroutine extras_binary_controls(binary_id, ierr)
82 : implicit none
83 : integer :: binary_id
84 : integer, intent(out) :: ierr
85 : end subroutine extras_binary_controls
86 :
87 : end interface
88 :
89 : integer, intent(out) :: ierr
90 : character (len = *) :: inlist_fname_arg
91 : optional inlist_fname_arg
92 :
93 : integer :: id, binary_id, i, j, k, l, i_prev, result, partial_result, &
94 : result_reason, model_number, iounit, binary_startup, model, num_stars
95 : type (star_info), pointer :: s
96 : character (len = 256) :: restart_filename, photo_filename
97 : integer(i8) :: time0, clock_rate
98 : logical :: doing_restart, first_try, continue_evolve_loop, &
99 : get_history_info, write_history, write_terminal, will_read_pgbinary_inlist
100 : type (binary_info), pointer :: b
101 : character (len = strlen) :: inlist_fname
102 :
103 : include 'formats'
104 :
105 0 : ierr = 0
106 0 : call system_clock(time0, clock_rate)
107 :
108 0 : call resolve_inlist_fname(inlist_fname, inlist_fname_arg)
109 0 : MESA_INLIST_RESOLVED = .true. ! Now any call to resolve_inlist_fname will only return inlist_fname_arg or 'inlist'
110 :
111 : ! Find out if this is a restart
112 0 : open(newunit = iounit, file = '.restart', status = 'old', action = 'read', iostat = ierr)
113 0 : doing_restart = (ierr == 0)
114 0 : if (doing_restart) then
115 : ! photo_filename will be like 'x100'
116 : ! the photo for star 1 is '1_x100', for star 2 '2_x100' and for binary 'b_x100'
117 0 : read(iounit, '(a)', iostat = ierr) photo_filename
118 0 : if (ierr /= 0) then
119 0 : stop "Problem while reading restart info"
120 : end if
121 : else
122 : ierr = 0
123 : end if
124 :
125 0 : binary_id = alloc_binary(ierr)
126 0 : if (ierr /= 0) then
127 0 : write(*, *) 'failed in alloc_binary'
128 0 : return
129 : end if
130 :
131 0 : call binary_ptr(binary_id, b, ierr)
132 0 : if (ierr /= 0) then
133 0 : write(*, *) 'failed in binary_ptr'
134 0 : return
135 : end if
136 :
137 0 : if (.not. doing_restart) then
138 0 : b% model_number = 0
139 0 : b% model_number_old = 0
140 0 : b% binary_age = 0
141 0 : b% binary_age_old = 0
142 : end if
143 :
144 0 : call do_read_binary_job(b, inlist_fname, ierr)
145 0 : if (ierr /= 0) then
146 0 : write(*, *) 'failed in read_binary_job'
147 0 : return
148 : end if
149 :
150 0 : if (.not. b% job% evolve_both_stars .or. &
151 : ((b% job% change_initial_model_twins_flag .or. b% job% change_model_twins_flag) &
152 : .and. b% job% new_model_twins_flag)) then
153 0 : b% have_star_1 = .true.
154 0 : b% have_star_2 = .false.
155 : else
156 0 : b% have_star_1 = .true.
157 0 : b% have_star_2 = .true.
158 : end if
159 :
160 0 : write(*, '(A)')
161 0 : write(*, '(A)')
162 :
163 0 : result_reason = 0
164 :
165 : ! Setup null hooks
166 0 : b% other_rlo_mdot => null_other_rlo_mdot
167 0 : b% other_implicit_function_to_solve => null_other_implicit_function_to_solve
168 0 : b% other_check_implicit_rlo => null_other_check_implicit_rlo
169 0 : b% other_tsync => null_other_tsync
170 0 : b% other_sync_spin_to_orbit => null_other_sync_spin_to_orbit
171 0 : b% other_mdot_edd => null_other_mdot_edd
172 0 : b% other_adjust_mdots => null_other_adjust_mdots
173 0 : b% other_accreted_material_j => null_other_accreted_material_j
174 0 : b% other_jdot_gr => null_other_jdot_gr
175 0 : b% other_jdot_ml => null_other_jdot_ml
176 0 : b% other_jdot_ls => null_other_jdot_ls
177 0 : b% other_jdot_missing_wind => null_other_jdot_missing_wind
178 0 : b% other_jdot_mb => null_other_jdot_mb
179 0 : b% other_extra_jdot => null_other_extra_jdot
180 0 : b% other_binary_wind_transfer => null_other_binary_wind_transfer
181 0 : b% other_edot_tidal => null_other_edot_tidal
182 0 : b% other_edot_enhance => null_other_edot_enhance
183 0 : b% other_extra_edot => null_other_extra_edot
184 0 : b% other_CE_init => null_other_CE_init
185 0 : b% other_CE_rlo_mdot => null_other_CE_rlo_mdot
186 0 : b% other_CE_binary_evolve_step => null_other_CE_binary_evolve_step
187 0 : b% other_CE_binary_finish_step => null_other_CE_binary_finish_step
188 0 : b% other_e2 => null_other_e2
189 0 : b% other_pgbinary_plots_info => null_other_pgbinary_plots_info
190 :
191 0 : b% extras_binary_startup => null_extras_binary_startup
192 0 : b% extras_binary_start_step => null_extras_binary_start_step
193 0 : b% extras_binary_check_model => null_extras_binary_check_model
194 0 : b% extras_binary_finish_step => null_extras_binary_finish_step
195 0 : b% extras_binary_after_evolve => null_extras_binary_after_evolve
196 0 : b% how_many_extra_binary_history_columns => null_how_many_extra_binary_history_columns
197 0 : b% data_for_extra_binary_history_columns => null_data_for_extra_binary_history_columns
198 0 : b% how_many_extra_binary_history_header_items => null_how_many_extra_binary_history_header_items
199 0 : b% data_for_extra_binary_history_header_items => null_data_for_extra_binary_history_header_items
200 :
201 0 : b% other_binary_photo_read => default_other_binary_photo_read
202 0 : b% other_binary_photo_write => default_other_binary_photo_write
203 :
204 0 : b% ignore_hard_limits_this_step = .false.
205 :
206 0 : call do_one_binary_setup(b, inlist_fname, ierr)
207 : ! extras_binary_controls is defined in run_binary_extras.f and hooks can
208 : ! be specified there
209 0 : call extras_binary_controls(b% binary_id, ierr)
210 :
211 : ! load binary photo if this is a restart
212 0 : if (doing_restart) then
213 0 : restart_filename = trim(trim(b% photo_directory) // '/b_' // photo_filename)
214 0 : call binary_load_photo(b, restart_filename, ierr)
215 0 : if (failed('binary_load_photo', ierr)) return
216 : end if
217 :
218 0 : b% donor_id = -1
219 0 : b% accretor_id = -1
220 0 : do i = 1, 2
221 :
222 0 : if (i==1 .and. .not. b% have_star_1) then
223 : cycle
224 0 : else if (i==2 .and. .not. b% have_star_2) then
225 : cycle
226 : end if
227 :
228 0 : call do_read_star_job(b% job% inlist_names(i), ierr)
229 0 : if (failed('do_read_star_job', ierr)) return
230 :
231 0 : if (i==1) then
232 0 : restart_filename = trim('1_' // photo_filename)
233 : else
234 0 : restart_filename = trim('2_' // photo_filename)
235 : end if
236 :
237 : ! the star is initialized in this call
238 : call before_evolve_loop(.true., doing_restart, doing_restart, &
239 : binary_controls, extras_controls, &
240 : id_from_read_star_job, b% job% inlist_names(i), restart_filename, &
241 0 : .false., binary_id, id, ierr)
242 :
243 0 : if (ierr /= 0) return
244 :
245 0 : call star_ptr(id, s, ierr)
246 0 : if (failed('star_ptr', ierr)) return
247 0 : b% star_ids(i) = id
248 :
249 0 : s% include_binary_history_in_log_file = b% append_to_star_history
250 :
251 0 : s% how_many_binary_history_columns => how_many_binary_history_columns
252 0 : s% data_for_binary_history_columns => data_for_binary_history_columns
253 :
254 0 : s% how_many_extra_binary_history_columns => b% how_many_extra_binary_history_columns
255 0 : s% data_for_extra_binary_history_columns => b% data_for_extra_binary_history_columns
256 :
257 : ! additional settings for mass transfer and tides
258 0 : if (b% do_j_accretion) then
259 0 : s% use_accreted_material_j = .true.
260 : end if
261 0 : s% accrete_given_mass_fractions = .true.
262 0 : s% accrete_same_as_surface = .false.
263 0 : s% binary_other_torque => sync_spin_orbit_torque
264 :
265 0 : s% doing_timing = .false.
266 :
267 0 : write(*, '(A)')
268 0 : write(*, '(A)')
269 :
270 : end do
271 :
272 : ! Error if users set the saved model to the same filename
273 0 : if(b% have_star_1 .and. b% have_star_2) then
274 0 : if(b% s1% job% save_model_when_terminate .and. b% s2% job% save_model_when_terminate) then
275 0 : if(len_trim(b% s1% job% save_model_filename) > 0 .and. len_trim(b% s2% job% save_model_filename) >0) then
276 0 : if(trim(b% s1% job% save_model_filename) == trim(b% s2% job% save_model_filename)) then
277 0 : write(*, *) "ERROR: Both stars are set to write save_model_filename to the same file"
278 0 : call mesa_error(__FILE__, __LINE__)
279 : end if
280 : end if
281 : end if
282 : end if
283 :
284 :
285 : ! binary data must be initiated after stars, such that masses are available
286 : ! if using saved models
287 0 : call binarydata_init(b, doing_restart)
288 0 : call binary_private_def_init
289 0 : call binary_history_column_names_init(ierr)
290 0 : call set_binary_history_columns(b, b% job% binary_history_columns_file, .true., ierr)
291 :
292 : ! setup pgbinary
293 0 : if (.not. doing_restart) then
294 0 : call start_new_run_for_pgbinary(b, ierr)
295 0 : if (failed('start_new_run_for_pgbinary', ierr)) return
296 : else
297 0 : call restart_run_for_pgbinary(b, ierr)
298 0 : if (failed('restart_run_for_pgbinary', ierr)) return
299 : end if
300 :
301 0 : if (b% job% show_binary_log_description_at_start .and. .not. doing_restart) then
302 0 : write(*, '(A)')
303 0 : call do_show_binary_log_description(id, ierr)
304 0 : if (failed('show_log_description', ierr)) return
305 : end if
306 :
307 0 : if (b% point_mass_i /= 1 .and. b% do_initial_orbit_sync_1 .and. &
308 : .not. doing_restart) then
309 : call star_relax_uniform_omega(&
310 : b% s1% id, 0, (2 * pi) / b% period, b% s1% job% num_steps_to_relax_rotation, &
311 0 : b% s1% job% relax_omega_max_yrs_dt, ierr)
312 0 : if (ierr /= 0) then
313 0 : write(*, *) 'failed in initial orbital sync'
314 0 : return
315 : end if
316 : end if
317 :
318 0 : if (b% point_mass_i /= 2 .and. b% do_initial_orbit_sync_2 .and. &
319 : .not. doing_restart) then
320 : call star_relax_uniform_omega(&
321 : b% s2% id, 0, (2 * pi) / b% period, b% s2% job% num_steps_to_relax_rotation, &
322 0 : b% s2% job% relax_omega_max_yrs_dt, ierr)
323 0 : if (ierr /= 0) then
324 0 : write(*, *) 'failed in initial orbital sync'
325 0 : return
326 : end if
327 : end if
328 :
329 0 : continue_evolve_loop = .true.
330 0 : s% doing_timing = .false.
331 0 : i_prev = 0
332 :
333 0 : binary_startup = b% extras_binary_startup(b% binary_id, doing_restart, ierr)
334 0 : if (ierr /= 0) return
335 :
336 : ! perform changes specified by binary_job variables
337 0 : call do_binary_job_controls_after(b% binary_id, b, doing_restart, ierr)
338 0 : if (ierr /= 0) return
339 :
340 : ! setup things for common envelope in case of a restart
341 0 : if (doing_restart .and. b% CE_flag .and. b% CE_init) then
342 0 : call ce_init(b, doing_restart, ierr)
343 : end if
344 :
345 0 : evolve_loop : do while(continue_evolve_loop) ! evolve one step per loop
346 :
347 0 : if (b% point_mass_i /= 0) then
348 : num_stars = 1
349 : else
350 0 : num_stars = 2
351 : end if
352 :
353 0 : do l = 1, num_stars
354 0 : if (l == 1 .and. b% point_mass_i == 1) then
355 : i = 2
356 : else
357 0 : i = l
358 : end if
359 :
360 0 : id = b% star_ids(i)
361 0 : call star_ptr(id, s, ierr)
362 0 : if (failed('star_ptr', ierr)) return
363 0 : call before_step_loop(id, ierr)
364 0 : if (ierr /= 0) return
365 :
366 0 : result = s% extras_start_step(id)
367 0 : if (result /= keep_going) then
368 0 : continue_evolve_loop = .false.
369 : exit evolve_loop
370 : end if
371 : end do
372 :
373 0 : first_try = .true.
374 0 : b% donor_started_implicit_wind = .false.
375 0 : b% num_tries = 0
376 :
377 0 : if (b% CE_flag .and. .not. b% CE_init) then
378 0 : call ce_init(b, .false., ierr)
379 0 : write(*, *) "CE flag is on!!"
380 : end if
381 :
382 0 : step_loop : do ! may need to repeat this loop
383 :
384 0 : result = b% extras_binary_start_step(b% binary_id, ierr)
385 0 : if (ierr /= 0) then
386 0 : write(*, *) "error in extras_binary_start_step"
387 0 : return
388 : end if
389 0 : if (result /= keep_going) exit evolve_loop
390 : !update num_stars in case point_mass_i has changed
391 0 : if (b% point_mass_i /= 0) then
392 : num_stars = 1
393 : else
394 0 : num_stars = 2
395 : end if
396 :
397 0 : call set_donor_star(b)
398 0 : call set_star_timesteps(b)
399 0 : result = keep_going
400 :
401 : ! Store mtransfer_rate used in a step, as it is rewritten by binary_check_model and
402 : ! that can produce inconsistent output.
403 0 : b% step_mtransfer_rate = b% mtransfer_rate
404 :
405 0 : if (.not. b% CE_flag) then
406 : ! if donor reaches implicit wind limit during mass transfer, set was_in_implicit_wind_limit = .true.
407 : ! so that further iterations of the implicit wind do not start far off from the required mdot.
408 : ! system will likely detach suddenly, so set change_factor to double its maximum
409 0 : if (b% donor_started_implicit_wind) then
410 0 : b% s_donor% was_in_implicit_wind_limit = .true.
411 0 : b% change_factor = 2 * b% max_change_factor
412 : end if
413 : end if
414 :
415 0 : do l = 1, num_stars
416 :
417 0 : if (b% point_mass_i /= 0) then ! if any point mass, evolve the other one
418 0 : j = 3 - b% point_mass_i
419 : else ! both stars
420 0 : if (l == 1) then
421 0 : j = b% d_i ! evolve the donor first
422 : else
423 0 : j = b% a_i
424 : end if
425 : end if
426 :
427 0 : id = b% star_ids(j)
428 :
429 : ! Avoid repeating the accretor when using the implicit scheme plus
430 : ! rotation and implicit winds. When this happens the accretor won't
431 : ! usually care about the result of the evolution of the donor.
432 0 : if (j == b% a_i .and. b% num_tries >0 .and. s% was_in_implicit_wind_limit) &
433 : cycle
434 :
435 : ! fix sync timescales to zero. If synching stars these will be
436 : ! updated at each star_evolve_step
437 0 : if (j == 1) then
438 0 : b% t_sync_1 = 0
439 0 : else if (j == 2) then
440 0 : b% t_sync_2 = 0
441 : end if
442 :
443 : result = worst_result(result, &
444 0 : star_evolve_step_part1(id, first_try))
445 :
446 : end do
447 :
448 : ! modify mdots to account for mass transfer
449 0 : call adjust_mdots(b)
450 : ! if both stars are accreting mass at this point, then there is
451 : ! something very wrong! If one star loses and the other gains mass,
452 : ! then the mass losing star must be evolved first
453 0 : k = b% d_i
454 0 : if (b% point_mass_i == 0 .and. b% s_donor% mstar_dot > 0) then
455 0 : if (b% s_accretor% mstar_dot > 0) then
456 0 : write(*, *) "ERROR: both stars accreting, terminating evolution"
457 0 : result = terminate
458 0 : exit step_loop
459 : end if
460 0 : k = b% a_i ! donor is gaining while accretor is losing, accretor plays donor in this step
461 0 : else if (b% point_mass_i /= 0 .and. b% s_donor% mstar_dot > 0) then
462 0 : write(*, *) "ERROR: donor accreting, terminating evolution"
463 0 : result = terminate
464 0 : exit step_loop
465 : end if
466 :
467 0 : if (result == keep_going) then
468 0 : do l = 1, num_stars
469 : ! if there's any point mass, evolve the non-point mass, num_stars should be 1 here
470 : ! so no risk of evolving things twice
471 0 : if (b% point_mass_i /= 0) then
472 0 : k = 3 - b% point_mass_i
473 : else ! two stars present
474 0 : if (l == 1) then ! evolve donor in first loop pass
475 0 : j = k
476 : else ! evolve acc in second loop pass
477 0 : j = 3 - k
478 : end if
479 : end if
480 0 : id = b% star_ids(j)
481 :
482 0 : if (j == b% a_i .and. b% num_tries >0 .and. s% was_in_implicit_wind_limit) &
483 : cycle
484 :
485 : ! set accretion composition
486 0 : if (.not. b% CE_flag) then
487 0 : if (i == 2) call set_accretion_composition(b, j)
488 : end if
489 :
490 : result = worst_result(result, &
491 0 : star_evolve_step_part2(id, first_try))
492 :
493 : end do
494 : end if
495 :
496 : ! do not evolve binary step on failure, its unnecessary and some variables are not properly
497 : ! set when the newton solver fails.
498 0 : if (result == keep_going) then
499 0 : if (.not. b% CE_flag) then
500 0 : result = worst_result(result, binary_evolve_step(b))
501 : else
502 0 : result = worst_result(result, CE_binary_evolve_step(b))
503 : end if
504 : end if
505 :
506 0 : if (result == keep_going) then
507 0 : result = worst_result(result, binary_check_model(b))
508 : end if
509 :
510 0 : do l = 1, num_stars
511 0 : if (l == 1 .and. b% point_mass_i == 1) then
512 : i = 2
513 : else
514 : i = l
515 : end if
516 : if (i == 1 .and. b% point_mass_i == 1) i = 2
517 0 : if (result == keep_going) then
518 0 : id = b% star_ids(i)
519 0 : call star_ptr(id, s, ierr)
520 0 : result = worst_result(result, s% extras_check_model(id))
521 0 : result = worst_result(result, star_check_model(id))
522 : end if
523 : end do
524 :
525 0 : partial_result = b% extras_binary_check_model(b% binary_id)
526 0 : result = worst_result(result, partial_result)
527 :
528 : ! solve first binary timestep limit because star_pick_next_timestep needs it
529 : ! check redos as well to stop the implicit solver if a hard limit is hit
530 0 : if (result == keep_going .or. result == redo) then
531 0 : result = worst_result(result, binary_pick_next_timestep(b))
532 : end if
533 :
534 0 : if (result == keep_going) then
535 0 : do l = 1, num_stars
536 0 : if (l == 1 .and. b% point_mass_i == 1) then
537 : i = 2
538 : else
539 0 : i = l
540 : end if
541 0 : id = b% star_ids(i)
542 0 : call star_ptr(id, s, ierr)
543 0 : if (failed('star_ptr', ierr)) return
544 0 : result = worst_result(result, star_pick_next_timestep(id))
545 : end do
546 : end if
547 0 : if (result == keep_going) then
548 : exit step_loop
549 : end if
550 :
551 0 : do l = 1, num_stars
552 0 : if (l == 1 .and. b% point_mass_i == 1) then
553 : i = 2
554 : else
555 0 : i = l
556 : end if
557 :
558 0 : id = b% star_ids(i)
559 0 : model_number = get_model_number(id, ierr)
560 0 : if (failed('get_model_number', ierr)) return
561 :
562 0 : result_reason = get_result_reason(id, ierr)
563 0 : if (result == retry) then
564 0 : if (failed('get_result_reason', ierr)) return
565 0 : if (s% job% report_retries) &
566 0 : write(*, '(i6,3x,a,/)') model_number, &
567 0 : 'retry reason ' // trim(result_reason_str(result_reason))
568 : end if
569 :
570 : end do
571 :
572 0 : b% generations = b% s_donor% generations
573 0 : if (result == redo) then
574 0 : do l = 1, num_stars
575 0 : if (l == 1 .and. b% point_mass_i == 1) then
576 : i = 2
577 : else
578 0 : i = l
579 : end if
580 0 : id = b% star_ids(i)
581 :
582 : ! Avoid repeating the accretor when using the implicit scheme plus
583 : ! rotation and implicit winds. When this happens the accretor won't
584 : ! usually care about the result of the evolution of the donor.
585 :
586 0 : if (i == b% a_i .and. b% num_tries >0 .and. s% was_in_implicit_wind_limit) &
587 : cycle
588 0 : result = worst_result(result, star_prepare_to_redo(id))
589 : end do
590 0 : result = worst_result(result, binary_prepare_to_redo(b))
591 : end if
592 0 : if (result == retry) then
593 0 : do l = 1, num_stars
594 0 : if (l == 1 .and. b% point_mass_i == 1) then
595 : i = 2
596 : else
597 0 : i = l
598 : end if
599 0 : id = b% star_ids(i)
600 0 : result = worst_result(result, star_prepare_to_retry(id))
601 : end do
602 0 : result = worst_result(result, binary_prepare_to_retry(b))
603 : end if
604 :
605 : !correct pointers to stars after retry
606 0 : if (b% d_i == 1) then
607 0 : if (b% have_star_1) then
608 0 : b% s_donor => b% s1
609 0 : if (b% point_mass_i == 0) then
610 0 : if (b% have_star_2) then
611 0 : b% s_accretor => b% s2
612 : else
613 0 : call mesa_error(__FILE__, __LINE__, 'ERROR: missing star pointer for accretor')
614 : end if
615 : end if
616 : else
617 0 : call mesa_error(__FILE__, __LINE__, 'ERROR: missing star pointer for donor')
618 : end if
619 : else
620 0 : if (b% have_star_2) then
621 0 : b% s_donor => b% s2
622 0 : if (b% point_mass_i == 0) then
623 0 : if (b% have_star_1) then
624 0 : b% s_accretor => b% s1
625 : else
626 0 : call mesa_error(__FILE__, __LINE__, 'ERROR: missing star pointer for accretor')
627 : end if
628 : end if
629 : else
630 0 : call mesa_error(__FILE__, __LINE__, 'ERROR: missing star pointer for donor')
631 : end if
632 : end if
633 :
634 0 : if (result == terminate) then
635 : continue_evolve_loop = .false.
636 : exit step_loop
637 : end if
638 0 : first_try = .false.
639 :
640 : end do step_loop
641 :
642 0 : partial_result = b% extras_binary_finish_step(b% binary_id)
643 0 : result = worst_result(result, partial_result)
644 :
645 0 : if(result == keep_going) result = binary_finish_step(b)
646 0 : if (b% CE_flag .and. b% CE_init .and. result == keep_going) then
647 0 : result = worst_result(result, CE_binary_finish_step(b))
648 : end if
649 :
650 0 : if (result == keep_going) then
651 : ! write terminal info
652 0 : model = b% model_number
653 0 : if (b% history_interval > 0) then
654 0 : write_history = (mod(model, b% history_interval) == 0)
655 : else
656 : write_history = .false.
657 : end if
658 0 : if (s% terminal_interval > 0) then
659 0 : write_terminal = (mod(model, b% terminal_interval) == 0)
660 : else
661 : write_terminal = .false.
662 : end if
663 0 : if (write_history) b% need_to_update_binary_history_now = .true.
664 0 : get_history_info = b% need_to_update_binary_history_now .or. write_terminal
665 : if (get_history_info) then
666 0 : if (b% write_header_frequency * b% terminal_interval > 0) then
667 : if (mod(model, b% write_header_frequency * b% terminal_interval) == 0 &
668 0 : .and. .not. b% doing_first_model_of_run) then
669 0 : write(*, '(A)')
670 0 : call write_binary_terminal_header(b)
671 : end if
672 : end if
673 0 : if (write_terminal) call do_binary_terminal_summary(b)
674 0 : if (b% need_to_update_binary_history_now) call write_binary_history_info(b, ierr)
675 0 : b% need_to_update_binary_history_now = .false.
676 : end if
677 : end if
678 :
679 0 : do l = 1, num_stars
680 0 : if (l == 1 .and. b% point_mass_i == 1) then
681 : i = 2
682 : else
683 0 : i = l
684 : end if
685 0 : id = b% star_ids(i)
686 0 : call star_ptr(id, s, ierr)
687 0 : partial_result = result
688 : call after_step_loop(id, b% job% inlist_names(i), &
689 0 : .false., partial_result, ierr)
690 0 : if (ierr /= 0) return
691 0 : result = worst_result(result, partial_result)
692 : end do
693 :
694 : ! do pgbinary stuff
695 0 : if (result == keep_going .and. b% job% pgbinary_flag) then
696 0 : will_read_pgbinary_inlist = .false.
697 0 : if (b% pg% pgbinary_interval <= 0) then
698 : will_read_pgbinary_inlist = .true.
699 0 : else if(mod(b% model_number, b% pg% pgbinary_interval) == 0) then
700 : will_read_pgbinary_inlist = .true.
701 : end if
702 0 : if (will_read_pgbinary_inlist) then
703 0 : call read_pgbinary_inlist(b, inlist_fname, ierr)
704 0 : if (failed('read_pgbinary_controls', ierr)) return
705 : end if
706 : end if
707 0 : if (result == keep_going .and. b% job% pgbinary_flag) then
708 0 : call update_pgbinary_plots(b, .false., ierr)
709 0 : if (failed('update_pgbinary_plots', ierr)) return
710 : end if
711 :
712 0 : if (result /= keep_going) then
713 0 : if (result /= terminate) then
714 0 : write(*, 2) 'ERROR in result value in run_star_extras: model', &
715 0 : s% model_number
716 0 : write(*, 2) 'result', result
717 0 : exit evolve_loop
718 : end if
719 0 : do l = 1, num_stars
720 0 : if (l == 1 .and. b% point_mass_i == 1) then
721 : i = 2
722 : else
723 0 : i = l
724 : end if
725 0 : id = b% star_ids(i)
726 0 : call star_ptr(id, s, ierr)
727 0 : if (s% result_reason == result_reason_normal) then
728 :
729 : partial_result = result
730 : call terminate_normal_evolve_loop(id, &
731 0 : .false., partial_result, ierr)
732 0 : if (ierr /= 0) return
733 0 : result = worst_result(result, partial_result)
734 :
735 : end if
736 : end do
737 0 : call write_binary_history_info(b, ierr)
738 0 : call do_binary_terminal_summary(b)
739 0 : exit evolve_loop
740 : end if
741 :
742 0 : do l = 1, num_stars
743 0 : if (l == 1 .and. b% point_mass_i == 1) then
744 : i = 2
745 : else
746 0 : i = l
747 : end if
748 0 : id = b% star_ids(i)
749 0 : call star_ptr(id, s, ierr)
750 :
751 : call do_saves(&
752 0 : id, ierr)
753 0 : if (ierr /= 0) return
754 :
755 0 : if (s% doing_timing) then
756 0 : call system_clock(s% job% time1_extra, s% job% clock_rate)
757 : s% job% after_step_timing = s% job% after_step_timing + &
758 0 : dble(s% job% time1_extra - s% job% time0_extra) / s% job% clock_rate
759 0 : s% job% check_time_end = eval_total_times(s% id, ierr)
760 : s% job% check_after_step_timing = s% job% check_after_step_timing + &
761 0 : (s% job% check_time_end - s% job% check_time_start)
762 : end if
763 : end do
764 0 : if (b% photo_interval > 0 .and. mod(b% model_number, b% photo_interval) == 0) then
765 0 : call do_saves_for_binary(b, ierr)
766 0 : if (ierr /= 0) return
767 : end if
768 :
769 0 : if (b% doing_first_model_of_run) b% doing_first_model_of_run = .false.
770 :
771 : end do evolve_loop
772 :
773 : ! save photos at end of the run
774 0 : call do_saves_for_binary(b, ierr)
775 0 : if (ierr /= 0) return
776 :
777 : ! deallocate arrays used for calculation of phase dependent variables
778 0 : deallocate(b% theta_co, b% time_co, b% mdot_donor_theta)
779 0 : deallocate(b% edot_theta, b% e1, b% e2, b% e3)
780 :
781 0 : ierr = binary_after_evolve(b)
782 0 : if (ierr /= 0) then
783 0 : write(*, *) "Error in binary_after_evolve"
784 0 : return
785 : end if
786 0 : call b% extras_binary_after_evolve(b% binary_id, ierr)
787 0 : if (ierr /= 0) then
788 0 : write(*, *) "Error in extras_binary_after_evolve"
789 0 : return
790 : end if
791 :
792 0 : do i = 1, num_stars
793 0 : if (i == 1 .and. b% point_mass_i == 1) then
794 : l = 2
795 : else
796 0 : l = i
797 : end if
798 :
799 0 : id = b% star_ids(l)
800 :
801 0 : call star_ptr(id, s, ierr)
802 0 : if (failed('star_ptr', ierr)) then
803 0 : ierr = 0
804 0 : cycle
805 : end if
806 :
807 0 : call after_evolve_loop(id, .true., ierr)
808 :
809 0 : if (s% doing_timing) then
810 0 : call system_clock(s% job% time1_extra, s% job% clock_rate)
811 : s% job% after_step_timing = s% job% after_step_timing + &
812 0 : dble(s% job% time1_extra - s% job% time0_extra) / s% job% clock_rate
813 0 : s% job% check_time_end = eval_total_times(s% id, ierr)
814 : s% job% check_after_step_timing = s% job% check_after_step_timing + &
815 0 : (s% job% check_time_end - s% job% check_time_start)
816 : end if
817 :
818 : end do
819 :
820 0 : call starlib_shutdown
821 :
822 0 : end subroutine do_run1_binary
823 :
824 0 : integer function worst_result(result1, result2)
825 : integer, intent(in) :: result1, result2
826 :
827 0 : if(result1 == terminate .or. result2 == terminate) then
828 0 : worst_result = terminate
829 : return
830 : end if
831 :
832 0 : if(result1 == retry .or. result2 == retry) then
833 0 : worst_result = retry
834 : return
835 : end if
836 :
837 0 : if(result1 == redo .or. result2 == redo) then
838 0 : worst_result = redo
839 0 : return
840 : end if
841 :
842 0 : worst_result = keep_going
843 : return
844 :
845 0 : end function worst_result
846 :
847 0 : subroutine binary_controls(id, binary_id, ierr)
848 : integer, intent(in) :: id, binary_id
849 : integer, intent(out) :: ierr
850 : type (star_info), pointer :: s
851 : type (binary_info), pointer :: b
852 : ierr = 0
853 :
854 0 : call star_ptr(id, s, ierr)
855 0 : if (ierr /= 0) then
856 0 : write(*, *) 'failed in star_ptr'
857 0 : return
858 : end if
859 0 : call binary_ptr(binary_id, b, ierr)
860 0 : if (ierr /= 0) then
861 0 : write(*, *) 'failed in binary_ptr'
862 0 : return
863 : end if
864 :
865 0 : s% binary_id = binary_id
866 : ! binary takes care of writing photos
867 0 : s% photo_interval = -1
868 0 : s% job% save_photo_when_terminate = .false.
869 0 : s% photo_digits = b% photo_digits
870 :
871 0 : if (b% donor_id == -1) then
872 0 : b% donor_id = id
873 0 : b% s1 => s
874 0 : s% initial_mass = b% m1
875 : else
876 0 : b% accretor_id = id
877 0 : b% s2 => s
878 0 : s% initial_mass = b% m2
879 : end if
880 : end subroutine binary_controls
881 :
882 0 : subroutine do_binary_job_controls_after(binary_id, b, restart, ierr)
883 : use binary_utils
884 : use binary_evolve, only : binary_finish_step
885 : integer, intent(in) :: binary_id
886 : type (binary_info), pointer :: b
887 : logical, intent(in) :: restart
888 : integer, intent(out) :: ierr
889 : logical :: need_to_fix_generations
890 :
891 : ! If something is changed in here, fix generations to zero so a retry wont
892 : ! erase change
893 0 : need_to_fix_generations = .false.
894 :
895 0 : if (b% job% change_ignore_rlof_flag .or. &
896 : (b% job% change_initial_ignore_rlof_flag .and. .not. restart)) then
897 0 : if (b% job% new_ignore_rlof_flag .neqv. b% ignore_rlof_flag) then
898 0 : call set_ignore_rlof_flag(binary_id, b% job% new_ignore_rlof_flag, ierr)
899 0 : if (ierr /= 0) then
900 0 : write(*, *) "error in set_ignore_rlof_flag"
901 0 : return
902 : end if
903 : need_to_fix_generations = .true.
904 : end if
905 : end if
906 :
907 0 : if (b% job% change_model_twins_flag .or. &
908 : (b% job% change_initial_model_twins_flag .and. .not. restart)) then
909 0 : if (b% job% new_model_twins_flag .neqv. b% model_twins_flag) then
910 0 : call set_model_twins_flag(binary_id, b% job% new_model_twins_flag, ierr)
911 0 : if (ierr /= 0) then
912 0 : write(*, *) "error in set_model_twins_flag"
913 0 : return
914 : end if
915 : need_to_fix_generations = .true.
916 : end if
917 : end if
918 :
919 0 : if (b% job% change_point_mass_i .or. &
920 : (b% job% change_initial_point_mass_i .and. .not. restart)) then
921 0 : if (b% job% new_point_mass_i /= b% point_mass_i) then
922 0 : call set_point_mass_i(binary_id, b% job% new_point_mass_i, ierr)
923 0 : if (ierr /= 0) then
924 0 : write(*, *) "error in set_point_mass_i"
925 0 : return
926 : end if
927 : need_to_fix_generations = .true.
928 : end if
929 : end if
930 :
931 0 : if (b% job% change_m1 .or. &
932 : (b% job% change_initial_m1 .and. .not. restart)) then
933 0 : if (b% m(1) /= b% job% new_m1) then
934 0 : call set_m1(binary_id, b% job% new_m1, ierr)
935 0 : if (ierr /= 0) then
936 0 : write(*, *) "error in set_m1"
937 0 : return
938 : end if
939 : need_to_fix_generations = .true.
940 : end if
941 : end if
942 :
943 0 : if (b% job% change_m2 .or. &
944 : (b% job% change_initial_m2 .and. .not. restart)) then
945 0 : if (b% m(2) /= b% job% new_m2) then
946 0 : call set_m2(binary_id, b% job% new_m2, ierr)
947 0 : if (ierr /= 0) then
948 0 : write(*, *) "error in set_m2"
949 0 : return
950 : end if
951 : need_to_fix_generations = .true.
952 : end if
953 : end if
954 :
955 0 : if (b% job% change_period_eccentricity .or. &
956 : (b% job% change_initial_period_eccentricity .and. .not. restart)) then
957 0 : if (b% period /= b% job% new_period * secday .or. &
958 : b% eccentricity /= b% job% new_eccentricity) then
959 : call set_period_eccentricity(binary_id, &
960 0 : b% job% new_period * secday, b% job% new_eccentricity, ierr)
961 0 : if (ierr /= 0) then
962 0 : write(*, *) "error in set_period_eccentricity"
963 0 : return
964 : end if
965 : need_to_fix_generations = .true.
966 : end if
967 : end if
968 :
969 0 : if (b% job% change_separation_eccentricity .or. &
970 : (b% job% change_initial_separation_eccentricity .and. .not. restart)) then
971 0 : if (b% separation /= b% job% new_separation * Rsun .or. &
972 : b% eccentricity /= b% job% new_eccentricity) then
973 : call set_separation_eccentricity(binary_id, &
974 0 : b% job% new_separation * Rsun, b% job% new_eccentricity, ierr)
975 0 : if (ierr /= 0) then
976 0 : write(*, *) "error in set_separation_eccentricity"
977 0 : return
978 : end if
979 : need_to_fix_generations = .true.
980 : end if
981 : end if
982 :
983 0 : if (need_to_fix_generations) then
984 0 : if (b% point_mass_i /= 1) then
985 0 : b% s1% generations = 0
986 : end if
987 0 : if (b% point_mass_i /= 2) then
988 0 : b% s2% generations = 0
989 : end if
990 0 : b% generations = 0
991 :
992 : ! this updates 'old' values in case there is a retry on the first step
993 0 : ierr = binary_finish_step(b)
994 :
995 0 : if (ierr == keep_going) then
996 0 : ierr = 0
997 : else
998 0 : ierr = -1
999 : end if
1000 : end if
1001 :
1002 0 : end subroutine do_binary_job_controls_after
1003 :
1004 : end module run_binary_support
|
