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