[ci optional] Merge branch 'main' into abundance_timestep_controls

Author: evbauer

.gitattributes

@@ -10,9 +10,13 @@
 *.mp4 filter=lfs diff=lfs merge=lfs -text
 
 # Try to get github to better understand what our filetypes are
-*.inc linguist-langauge=Fortran
-*inlist* linguist-langauge=Fortran
+*.inc linguist-language=Fortran
+*inlist* linguist-language=Fortran
 *.mod linguist-generated
+*reaclib_results* -linguist-detectable
+*.ph -linguist-detectable
+*.[1-9] -linguist-detectable
+*.in -linguist-detectable
 
 # Things to ignore from the release
 star/dev_cases_compare_pulses export-ignore

.github/workflows/build_compile.yml

@@ -94,19 +94,22 @@ @ARTICLE{Paxton2019
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 
-@ARTICLE{Jermyn2022,
-       author = {{Jermyn}, Adam S. and {Bauer}, Evan B. and {Schwab}, Josiah and {Farmer}, R. and {Ball}, Warrick H. and {Bellinger}, Earl P. and {Dotter}, Aaron and {Joyce}, Meridith and {Marchant}, Pablo and {Mombarg}, Joey S.~G. and {Wolf}, William M. and {Wong}, Tin Long Sunny and {Cinquegrana}, Giulia C. and {Farrell}, Eoin and {Smolec}, R. and {Thoul}, Anne and {Cantiello}, Matteo and {Herwig}, Falk and {Toloza}, Odette and {Bildsten}, Lars and {Townsend}, Richard H.~D. and {Timmes}, F.~X.},
-        title = "{Modules for Experiments in Stellar Astrophysics (MESA): Time-Dependent Convection, Energy Conservation, Automatic Differentiation, and Infrastructure}",
-      journal = {arXiv e-prints},
-     keywords = {Astrophysics - Solar and Stellar Astrophysics, Astrophysics - Earth and Planetary Astrophysics, Astrophysics - High Energy Astrophysical Phenomena, Astrophysics - Instrumentation and Methods for Astrophysics},
-         year = 2022,
-        month = {dec},
-          eid = {arXiv:2208.03651},
-        pages = {arXiv:2208.03651},
+@ARTICLE{Jermyn2023,
+       author = {{Jermyn}, Adam S. and {Bauer}, Evan B. and {Schwab}, Josiah and {Farmer}, R. and {Ball}, Warrick H. and {Bellinger}, Earl P. and {Dotter}, Aaron and {Joyce}, Meridith and {Marchant}, Pablo and {Mombarg}, Joey S.~G. and {Wolf}, William M. and {Sunny Wong}, Tin Long and {Cinquegrana}, Giulia C. and {Farrell}, Eoin and {Smolec}, R. and {Thoul}, Anne and {Cantiello}, Matteo and {Herwig}, Falk and {Toloza}, Odette and {Bildsten}, Lars and {Townsend}, Richard H.~D. and {Timmes}, F.~X.},
+        title = "{Modules for Experiments in Stellar Astrophysics (MESA): Time-dependent Convection, Energy Conservation, Automatic Differentiation, and Infrastructure}",
+      journal = {\apjs},
+     keywords = {Stellar physics, Stellar evolution, Stellar evolutionary models, Computational methods, 1621, 1599, 2046, 1965, Astrophysics - Solar and Stellar Astrophysics, Astrophysics - Earth and Planetary Astrophysics, Astrophysics - High Energy Astrophysical Phenomena, Astrophysics - Instrumentation and Methods for Astrophysics},
+         year = 2023,
+        month = mar,
+       volume = {265},
+       number = {1},
+          eid = {15},
+        pages = {15},
+          doi = {10.3847/1538-4365/acae8d},
 archivePrefix = {arXiv},
        eprint = {2208.03651},
  primaryClass = {astro-ph.SR},
-       adsurl = {https://ui.adsabs.harvard.edu/abs/2022arXiv220803651J},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2023ApJS..265...15J},
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 

CITATIONS.bib

@@ -813,15 +813,19 @@
 ! wind_BH_alpha_*
 ! ~~~~~~~~~~~~~~~
 
-! Bondi-Hoyle accretion parameter for each star. The default is 3/2 taken
-! from Hurley et al. 2002, MNRAS, 329, 897, in agreement with
-! Boffin & Jorissen 1988, A&A, 205, 155
+! Bondi-Hoyle accretion parameter for each star. The default for alpha is 3/2
+! taken from Hurley et al. 2002, MNRAS, 329, 897, in agreement with
+! Boffin & Jorissen 1988, A&A, 205, 155.
+! The default for beta is 1/8=0.125 in accordance for results of cool
+! supergiants from Kucinskas A., 1999, Ap&SS, 262, 127
 ! "\_1" refers to first star, "\_2" to the second one.
 
 ! ::
 
     wind_BH_alpha_1 = 1.5d0
     wind_BH_alpha_2 = 1.5d0
+    wind_BH_beta_1 = 1.25d-1
+    wind_BH_beta_2 = 1.25d-1
 
 
 ! max_wind_transfer_fraction_*

binary/defaults/binary_controls.defaults

@@ -125,6 +125,8 @@ module binary_ctrls_io
          do_wind_mass_transfer_2, &
          wind_BH_alpha_1, &
          wind_BH_alpha_2, &
+         wind_BH_beta_1, &
+         wind_BH_beta_2, &
          max_wind_transfer_fraction_1, &
          max_wind_transfer_fraction_2, &
 
@@ -428,6 +430,8 @@ subroutine store_binary_controls(b, ierr)
          b% do_wind_mass_transfer_2 = do_wind_mass_transfer_2
          b% wind_BH_alpha_1 = wind_BH_alpha_1
          b% wind_BH_alpha_2 = wind_BH_alpha_2
+         b% wind_BH_beta_1 = wind_BH_beta_1
+         b% wind_BH_beta_2 = wind_BH_beta_2
          b% max_wind_transfer_fraction_1 = max_wind_transfer_fraction_1
          b% max_wind_transfer_fraction_2 = max_wind_transfer_fraction_2
 
@@ -618,6 +622,8 @@ subroutine set_binary_controls_for_writing(b, ierr)
          do_wind_mass_transfer_2 = b% do_wind_mass_transfer_2
          wind_BH_alpha_1 = b% wind_BH_alpha_1
          wind_BH_alpha_2 = b% wind_BH_alpha_2
+         wind_BH_beta_1 = b% wind_BH_beta_1
+         wind_BH_beta_2 = b% wind_BH_beta_2
          max_wind_transfer_fraction_1 = b% max_wind_transfer_fraction_1
          max_wind_transfer_fraction_2 = b% max_wind_transfer_fraction_2
 

binary/private/binary_ctrls_io.f90

@@ -188,6 +188,7 @@ subroutine binarydata_init(b, doing_restart)
             b% CE_lambda2 = 0d0
             b% CE_Ebind1 = 0d0
             b% CE_Ebind2 = 0d0
+            b% mtransfer_rate = 0
 
             b% num_tries = 0
 
@@ -337,6 +338,11 @@ integer function binary_evolve_step(b)
             b% eccentricity = b% eccentricity + get_edot(b) *b% time_step*secyer
             if (b% eccentricity < b% min_eccentricity) b% eccentricity = b% min_eccentricity
             if (b% eccentricity > b% max_eccentricity) b% eccentricity = b% max_eccentricity
+         else
+            b% edot_tidal = 0d0
+            b% edot_enhance = 0d0
+            b% extra_edot = 0d0
+            b% edot = 0d0
          end if
 
          !use new eccentricity to calculate new time coordinate
@@ -504,17 +510,19 @@ integer function binary_finish_step(b)
 
          binary_finish_step = keep_going
          ! update change factor in case mtransfer_rate has changed
-         if(b% mtransfer_rate_old /= b% mtransfer_rate .and. &
-             b% mtransfer_rate /= 0 .and. b% mtransfer_rate_old /= 0) then
-            if(b% mtransfer_rate < b% mtransfer_rate_old) then
-               b% change_factor = b% change_factor*(1d0-b% implicit_lambda) + b% implicit_lambda* &
-                  (1+b% change_factor_fraction*(b% mtransfer_rate/b% mtransfer_rate_old-1))
-            else
-               b% change_factor = b% change_factor*(1d0-b% implicit_lambda) + b% implicit_lambda* &
-                   (1+b% change_factor_fraction*(b% mtransfer_rate_old/b% mtransfer_rate-1))
+         if(.not. b% doing_first_model_of_run) then
+            if(b% mtransfer_rate_old /= b% mtransfer_rate .and. &
+               b% mtransfer_rate /= 0 .and. b% mtransfer_rate_old /= 0) then
+               if(b% mtransfer_rate < b% mtransfer_rate_old) then
+                  b% change_factor = b% change_factor*(1d0-b% implicit_lambda) + b% implicit_lambda* &
+                     (1+b% change_factor_fraction*(b% mtransfer_rate/b% mtransfer_rate_old-1))
+               else
+                  b% change_factor = b% change_factor*(1d0-b% implicit_lambda) + b% implicit_lambda* &
+                     (1+b% change_factor_fraction*(b% mtransfer_rate_old/b% mtransfer_rate-1))
+               end if
+               if(b% change_factor > b% max_change_factor) b% change_factor = b% max_change_factor
+               if(b% change_factor < b% min_change_factor) b% change_factor = b% min_change_factor
             end if
-            if(b% change_factor > b% max_change_factor) b% change_factor = b% max_change_factor
-            if(b% change_factor < b% min_change_factor) b% change_factor = b% min_change_factor
          end if
 
          ! store all variables into "old"

binary/private/binary_evolve.f90

@@ -126,39 +126,46 @@ integer function binary_pick_next_timestep(b)
 
          if (b% max_timestep < 0) b% max_timestep = b% s_donor% dt
 
-
          b% env(b% d_i) = s% star_mass - s% he_core_mass
          if (b% point_mass_i == 0) then
             b% env(b% a_i) = b% s_accretor% star_mass - b% s_accretor% he_core_mass
          end if
 
-         if (b% env_old(b% d_i) /= 0) then
-            env_change = b% env(b% d_i) - b% env_old(b% d_i)
+         if(.not. b% doing_first_model_of_run) then
+            if (b% env_old(b% d_i) /= 0) then
+               env_change = b% env(b% d_i) - b% env_old(b% d_i)
+            else
+               env_change = 0
+            end if
+            
+            if (b% rl_relative_gap_old(b% d_i) /= 0) then
+               rel_gap_change = b% rl_relative_gap_old(b% d_i) - b% rl_relative_gap(b% d_i)
+            else
+               rel_gap_change = 0
+            end if
+            
+            if (b% angular_momentum_j_old /= 0) then
+               j_change = b% angular_momentum_j - b% angular_momentum_j_old
+            else
+               j_change = 0
+            end if
+            
+            if (b% separation_old /= 0) then
+               sep_change = b% separation - b% separation_old
+            else
+               sep_change = 0
+            end if
+            if (b% eccentricity_old /= 0) then
+               e_change = b% eccentricity - b% eccentricity_old
+            else
+               e_change = 0
+            end if
          else
             env_change = 0
-         end if
-         
-         if (b% rl_relative_gap_old(b% d_i) /= 0) then
-            rel_gap_change = b% rl_relative_gap_old(b% d_i) - b% rl_relative_gap(b% d_i)
-         else
             rel_gap_change = 0
-         end if
-         
-         if (b% angular_momentum_j_old /= 0) then
-            j_change = b% angular_momentum_j - b% angular_momentum_j_old
-         else
             j_change = 0
-         end if
-         
-         if (b% separation_old /= 0) then
-            sep_change = b% separation - b% separation_old
-         else
             sep_change = 0
-         end if
-         if (b% eccentricity_old /= 0) then
-             e_change = b% eccentricity - b% eccentricity_old
-         else
-             e_change = 0
+            e_change = 0
          end if
 
          ! get limits for dt based on relative changes

binary/private/binary_timestep.f90

@@ -94,8 +94,8 @@ subroutine Bondi_Hoyle_wind_transfer(binary_id, s_i, ierr)
 
       type(binary_info), pointer :: b
       type (star_info), pointer :: s
-      real(dp) :: v_orb, v_wind, b_BH
-      real(dp) :: alpha  ! Bondi-Hoyle alpha for that star
+      real(dp) :: v_orb, v_wind
+      real(dp) :: alpha, beta  ! Bondi-Hoyle alpha, beta for that star
       real(dp) :: max_xfer  ! Maximum transfer fraction
 
       call binary_ptr(binary_id, b, ierr)
@@ -107,18 +107,20 @@ subroutine Bondi_Hoyle_wind_transfer(binary_id, s_i, ierr)
       if (s_i == 1) then
          s => b% s1
          alpha = b% wind_BH_alpha_1
+         beta = b% wind_BH_beta_1
          max_xfer = b% max_wind_transfer_fraction_1
       else
          s => b% s2
          alpha = b% wind_BH_alpha_2
+         beta = b% wind_BH_beta_2
          max_xfer = b% max_wind_transfer_fraction_2
       end if
 
       ! orbital speed Hurley et al 2002 eq. 8
-      v_orb = sqrt(standard_cgrav * b% m(s_i) / b% separation) !cm/s
+      v_orb = sqrt(standard_cgrav * (b% m(1) + b% m(2)) / b% separation) !cm/s
 
       ! windspeed from Hurley et al 2002 eq. 9
-      v_wind = sqrt( 2d0 / 8d0 *  standard_cgrav * b% m(s_i) / b% r(s_i) )
+      v_wind = sqrt(2d0 * beta *  standard_cgrav * b% m(s_i) / b% r(s_i))
 
       ! Bondi-Hoyle transfer fraction Hurley et al. 2002 eq. 6
       b% wind_xfer_fraction(s_i) = alpha / pow2(b% separation) /&

binary/private/binary_wind.f90

@@ -84,6 +84,7 @@ logical :: do_enhance_wind_1, do_enhance_wind_2
 real(dp) :: tout_B_wind_1, tout_B_wind_2
 logical :: do_wind_mass_transfer_1, do_wind_mass_transfer_2
 real(dp) :: wind_BH_alpha_1, wind_BH_alpha_2
+real(dp) :: wind_BH_beta_1, wind_BH_beta_2
 real(dp) :: max_wind_transfer_fraction_1, max_wind_transfer_fraction_2
 
 ! orbital jdot controls

binary/public/binary_controls.inc

@@ -22,6 +22,8 @@
 
    limit_retention_by_mdot_edd = .true.
    do_wind_mass_transfer_1 = .true.
+   wind_BH_beta_1 = 1.25d-1 ! wind_BH_beta_1 = 7d0  ! 7 since donor is O star, using this value would make more sense, but we keep the default of 1/8
+                           ! which results in wind accretion below and above the Eddington limit (which we want to test).
    do_jdot_mb = .false.
 
    report_rlo_solver_progress = .true.