diff --git a/binary/private/binary_evolve.f90 b/binary/private/binary_evolve.f90 index 796eabcf8..41fdddc34 100644 --- a/binary/private/binary_evolve.f90 +++ b/binary/private/binary_evolve.f90 @@ -206,21 +206,17 @@ subroutine binarydata_init(b, doing_restart) end subroutine binarydata_init subroutine set_donor_star(b) + use binary_mdot, only: switch_donor + type (binary_info), pointer :: b - logical :: switch_donor real(dp) :: mdot_hi_temp include 'formats' - switch_donor = .false. - if (b% keep_donor_fixed .and. b% mdot_scheme /= "contact") return - if (b% mdot_scheme == "roche_lobe" .and. & - abs(b% mtransfer_rate/(Msun/secyer)) < b% mdot_limit_donor_switch .and. & - b% rl_relative_gap_old(b% a_i) > b% rl_relative_gap_old(b% d_i)) then - switch_donor = .true. - else if (b% mtransfer_rate > 0d0) then - switch_donor = .true. + ! main reason to change donor: the implicit solver requires a "positive" MT rate + ! switch donor and flip the signs of the implicit solver values + if (b% mtransfer_rate > 0d0) then b% mtransfer_rate = - b% mtransfer_rate mdot_hi_temp = b% mdot_hi b% mdot_hi = - b% mdot_lo @@ -230,28 +226,25 @@ subroutine set_donor_star(b) end if b% have_mdot_lo = .true. b% fixed_delta_mdot = b% fixed_delta_mdot / 2.0d0 - else if (b% mdot_scheme == "contact" .and. & - b% rl_relative_gap_old(b% a_i) > b% rl_relative_gap_old(b% d_i) .and. & - b% rl_relative_gap_old(b% a_i) < - b% implicit_scheme_tolerance .and. & - b% rl_relative_gap_old(b% d_i) < - b% implicit_scheme_tolerance .and. & - abs(b% mtransfer_rate/(Msun/secyer)) < b% mdot_limit_donor_switch) then - switch_donor = .true. + + call switch_donor(b) + return end if - if (switch_donor) then - if (b% report_rlo_solver_progress) write(*,*) "switching donor" - if (b% d_i == 2) then - b% d_i = 1 - b% a_i = 2 - b% s_donor => b% s1 - b% s_accretor => b% s2 - else - b% d_i = 2 - b% a_i = 1 - b% s_donor => b% s2 - b% s_accretor => b% s1 - end if + ! if accretor is bigger than the donor, switch donor, but only when the MT rate is low to avoid erratic switching + if ((b% mdot_scheme == "roche_lobe" .and. & + abs(b% mtransfer_rate/(Msun/secyer)) < b% mdot_limit_donor_switch .and. & + b% rl_relative_gap_old(b% a_i) > b% rl_relative_gap_old(b% d_i)) & + .or. & ! same in the contact scheme with both stars under the RL + (b% mdot_scheme == "contact" .and. & + b% rl_relative_gap_old(b% a_i) > b% rl_relative_gap_old(b% d_i) .and. & + b% rl_relative_gap_old(b% a_i) < - b% implicit_scheme_tolerance .and. & + b% rl_relative_gap_old(b% d_i) < - b% implicit_scheme_tolerance .and. & + abs(b% mtransfer_rate/(Msun/secyer)) < b% mdot_limit_donor_switch)) then + + call switch_donor(b) end if + end subroutine set_donor_star integer function binary_evolve_step(b) diff --git a/binary/private/binary_mdot.f90 b/binary/private/binary_mdot.f90 index 4af98df65..c2bb7a3e7 100644 --- a/binary/private/binary_mdot.f90 +++ b/binary/private/binary_mdot.f90 @@ -492,19 +492,38 @@ subroutine eval_mdot_edd(binary_id, mdot_edd, mdot_edd_eta, ierr) end subroutine eval_mdot_edd + ! switch donor and accretor pointers + subroutine switch_donor(b) + type (binary_info), pointer :: b + + if (b% d_i == 2) then + b% d_i = 1 + b% a_i = 2 + b% s_donor => b% s1 + b% s_accretor => b% s2 + else + b% d_i = 2 + b% a_i = 1 + b% s_donor => b% s2 + b% s_accretor => b% s1 + end if + end subroutine switch_donor + subroutine adjust_mdots(b) use binary_wind, only: eval_wind_xfer_fractions type (binary_info), pointer :: b - real(dp) :: actual_mtransfer_rate + real(dp) :: actual_mtransfer_rate, mdot_donor, mdot_acc integer :: ierr + include 'formats' + actual_mtransfer_rate = 0d0 if (b% use_other_adjust_mdots) then call b% other_adjust_mdots(b% binary_id, ierr) if (ierr /= 0) then - write(*,*) "Error in other_adjust_mdots" + write(*, 1) "Error in other_adjust_mdots" stop end if return @@ -523,7 +542,7 @@ subroutine adjust_mdots(b) call Tout_enhance_wind(b, b% s_donor) if (b% point_mass_i == 0) then ! do not repeat if using the implicit wind - if (.not. (b% num_tries >0 .and. b% s_accretor% was_in_implicit_wind_limit)) & + if (.not. (b% num_tries > 0 .and. b% s_accretor% was_in_implicit_wind_limit)) & call Tout_enhance_wind(b, b% s_accretor) end if @@ -533,7 +552,7 @@ subroutine adjust_mdots(b) ! accretor. call eval_wind_xfer_fractions(b% binary_id, ierr) if (ierr/=0) then - write(*,*) "Error in eval_wind_xfer_fractions" + write(*, 1) "Error in eval_wind_xfer_fractions" return end if b% mdot_wind_transfer(b% d_i) = b% s_donor% mstar_dot * & @@ -545,20 +564,37 @@ subroutine adjust_mdots(b) b% mdot_wind_transfer(b% a_i) = 0d0 end if - ! Set mdot for the donor + ! compute total mdots in current configuration + mdot_donor = b% s_donor% mstar_dot + b% mtransfer_rate - & + b% mdot_wind_transfer(b% a_i) + if (b% point_mass_i == 0) then + mdot_acc = b% s_accretor% mstar_dot - & + b% mtransfer_rate*b% fixed_xfer_fraction - b% mdot_wind_transfer(b% d_i) + else + mdot_acc = -b% mtransfer_rate*b% fixed_xfer_fraction - b% mdot_wind_transfer(b% d_i) + end if + + ! if accretor ended up losing mass while the donor is gaining it, switch donor.\\ + ! this is likely due to wind-mass transfer, b% mtransfer_rate is probably low + ! we do not switch the sign of mtransfer rate here. + if (mdot_donor > 0d0 .and. mdot_acc < 0d0) then + if (b% report_rlo_solver_progress) write(*, 1) "adjust_mdots: switching donor due to accreting donor and donating accretor" + call switch_donor(b) + end if + + ! now set mdot for the (possibly new) donor b% s_donor% mstar_dot = b% s_donor% mstar_dot + b% mtransfer_rate - & b% mdot_wind_transfer(b% a_i) ! Set mdot for the accretor if (b% point_mass_i == 0 .and. .not. b% CE_flag) then ! do not repeat if using the implicit wind - if (.not. (b% num_tries >0 .and. b% s_accretor% was_in_implicit_wind_limit)) then + if (.not. (b% num_tries > 0 .and. b% s_accretor% was_in_implicit_wind_limit)) then b% accretion_mode = 0 b% acc_am_div_kep_am = 0.0d0 b% s_accretor% mstar_dot = b% s_accretor% mstar_dot - & b% mtransfer_rate*b% fixed_xfer_fraction - b% mdot_wind_transfer(b% d_i) - !set angular momentum accretion as described in A.3.3 of de Mink et al. 2013 if (b% do_j_accretion) then if (.not. b% use_other_accreted_material_j) then call eval_accreted_material_j(b% binary_id, ierr) @@ -566,7 +602,7 @@ subroutine adjust_mdots(b) call b% other_accreted_material_j(b% binary_id, ierr) end if if (ierr /= 0) then - write(*,*) 'error in accreted_material_j' + write(*, 1) 'error in accreted_material_j' return end if end if @@ -947,11 +983,12 @@ subroutine get_info_for_kolb_eccentric(b) end subroutine get_info_for_kolb_eccentric + !Angular momentum accretion as described in A.3.3 of de Mink et al. 2013 subroutine eval_accreted_material_j(binary_id, ierr) integer, intent(in) :: binary_id integer, intent(out) :: ierr type(binary_info), pointer :: b - real(dp) :: qratio, min_r + real(dp) :: q, min_r logical, parameter :: dbg = .false. include 'formats' @@ -961,18 +998,20 @@ subroutine eval_accreted_material_j(binary_id, ierr) write(*,*) 'failed in binary_ptr' return end if - qratio = b% m(b% a_i) / b% m(b% d_i) - qratio = min(max(qratio,0.0667d0),15d0) - min_r = 0.0425d0*b% separation*pow(qratio+qratio*qratio, 0.25d0) - !TODO: MUST USE EQUATORIAL RADIUS - if (dbg) write(*,*) "radius, impact_radius, separation: ", & + ! minimum approach of the gas stream from Ulrich & Burger, 1976, ApJ 206:509-514 + q = b% m(b% a_i) / b% m(b% d_i) + q = min(max(q, 0.0667d0), 15d0) + min_r = 0.0425d0 * b% separation * pow(q + q * q, 0.25d0) + + !TODO: find an appropriate way to define an "equatorial radius" in a binary + if (dbg) write(*, *) "radius, impact_radius, separation: ", & b% r(b% a_i), min_r/rsun, b% separation/rsun - if (b% r(b% a_i) < min_r) then + if (b% r(b% a_i) < min_r) then ! accretion through disk; inner rim has j = sqrt(GMR) b% accretion_mode = 2 b% s_accretor% accreted_material_j = & sqrt(standard_cgrav * b% m(b% a_i) * b% r(b% a_i)) - else + else ! direct impact; stream has j \approx sqrt(GM * 1.7 Rmin) b% accretion_mode = 1 b% s_accretor% accreted_material_j = & sqrt(standard_cgrav * b% m(b% a_i) * 1.7d0*min_r) @@ -980,9 +1019,6 @@ subroutine eval_accreted_material_j(binary_id, ierr) b% acc_am_div_kep_am = b% s_accretor% accreted_material_j / & sqrt(standard_cgrav * b% m(b% a_i) * b% r(b% a_i)) - !TODO: when using wind mass transfer donor star can end up - ! with positive mdot, need to properly set jdot in that case - end subroutine eval_accreted_material_j subroutine set_accretion_composition(b, acc_index) diff --git a/binary/private/binary_wind.f90 b/binary/private/binary_wind.f90 index cab4e254c..db21c12c5 100644 --- a/binary/private/binary_wind.f90 +++ b/binary/private/binary_wind.f90 @@ -44,7 +44,7 @@ subroutine eval_wind_xfer_fractions(binary_id, ierr) if (b% point_mass_i /= 1) then if (.not. b% do_wind_mass_transfer_1 .or. b% model_twins_flag) then b% wind_xfer_fraction(1) = 0d0 - else if(.not. b% use_other_binary_wind_transfer) then + else if (.not. b% use_other_binary_wind_transfer) then call Bondi_Hoyle_wind_transfer(b% binary_id, 1, ierr) if (ierr /=0) then write(*,*) "Error in Bondi_Hoyle_wind_transfer(b% binary_id, 1, ierr)" @@ -52,7 +52,7 @@ subroutine eval_wind_xfer_fractions(binary_id, ierr) end if else call b% other_binary_wind_transfer(b% binary_id, 1, ierr) - if (ierr /=0) then + if (ierr /= 0) then write(*,*) "Error in other_binary_wind_transfer(b% binary_id, 1, ierr)" return end if @@ -63,7 +63,7 @@ subroutine eval_wind_xfer_fractions(binary_id, ierr) if (b% point_mass_i /= 2) then if (.not. b% do_wind_mass_transfer_2) then b% wind_xfer_fraction(2) = 0d0 - else if(.not. b% use_other_binary_wind_transfer) then + else if (.not. b% use_other_binary_wind_transfer) then call Bondi_Hoyle_wind_transfer(b% binary_id, 2, ierr) if (ierr /=0) then write(*,*) "Error in Bondi_Hoyle_wind_transfer(b% binary_id, 2, ierr)" @@ -71,7 +71,7 @@ subroutine eval_wind_xfer_fractions(binary_id, ierr) end if else call b% other_binary_wind_transfer(b% binary_id, 2, ierr) - if (ierr /=0) then + if (ierr /= 0) then write(*,*) "Error in other_binary_wind_transfer(b% binary_id, 2, ierr)" return end if diff --git a/binary/private/run_binary_support.f90 b/binary/private/run_binary_support.f90 index 3697b93ad..7dd87243b 100644 --- a/binary/private/run_binary_support.f90 +++ b/binary/private/run_binary_support.f90 @@ -425,11 +425,13 @@ end subroutine extras_binary_controls end if id = b% star_ids(j) + call star_ptr(id, s, ierr) + if (failed('star_ptr', ierr)) return ! Avoid repeating the accretor when using the implicit scheme plus ! rotation and implicit winds. When this happens the accretor won't ! usually care about the result of the evolution of the donor. - if (j == b% a_i .and. b% num_tries >0 .and. s% was_in_implicit_wind_limit) & + if (j == b% a_i .and. b% num_tries > 0 .and. s% was_in_implicit_wind_limit) & cycle ! fix sync timescales to zero. If synching stars these will be @@ -440,6 +442,7 @@ end subroutine extras_binary_controls b% t_sync_2 = 0 end if + ! evaluates mdot from wind for either star result = worst_result(result, & star_evolve_step_part1(id, first_try)) @@ -477,9 +480,12 @@ end subroutine extras_binary_controls j = 3 - k end if end if + id = b% star_ids(j) + call star_ptr(id, s, ierr) + if (failed('star_ptr', ierr)) return - if (j == b% a_i .and. b% num_tries >0 .and. s% was_in_implicit_wind_limit) & + if (j == b% a_i .and. b% num_tries > 0 .and. s% was_in_implicit_wind_limit) & cycle ! set accretion composition