Converged develop with optimize: split alpha/acc_lim

docs/pages/evolve/interface.rst

@@ -236,7 +236,7 @@ new settings. For example, let's say we want to see how changing the common enve
 .. ipython:: python
 
     ce_alpha_10 = output.rerun_with_settings(
-        new_settings={'alpha1': 10}, inplace=False
+        new_settings={'alpha1_0': 10, 'alpha1_1': 10}, inplace=False
     )
     n_merger_original = len(output.final_bpp[output.final_bpp["sep"] == 0.0])
     n_merger_ce_alpha_10 = len(ce_alpha_10.final_bpp[ce_alpha_10.final_bpp["sep"] == 0.0])

examples/Params.ini

@@ -360,7 +360,7 @@ LBV_flag = 1
 ;    positive values - Sets the common-envelope efficiency parameter, \(\alpha\)
 ;    1.0 - Default value
 ; Default: 1.0
-alpha1 = 1.0
+alpha1 = [1.0, 1.0]
 
 ; lambdaf
 ; Binding energy factor for common envelope evolution. The initial binding energy of the stellar envelope is proportional to \(1 / \lambda\). See Hurley+2002, Eq. 69.
@@ -454,7 +454,7 @@ qcrit_array = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
 ;    -5 - As 5, but using the old Kiel & Hurley 2009 prescription for changing the orbital configuration of the binary, available for reproducibility purposes but not recommended for new work
 ;    -6 - As 6, but using the old Kiel & Hurley 2009 prescription for changing the orbital configuration of the binary, available for reproducibility purposes but not recommended for new work
 ; Default: 5
-kickflag = 5
+kickflag = 1
 
 ; sigma
 ; Sets the dispersion in the Maxwellian for the SN kick velocity in km/s
@@ -797,7 +797,7 @@ don_lim = -1
 ;    positive values - sets overall fraction of donor material that is accreted, with the rest being lost from the system (``acc_lim = 0.5`` assumes 50% accretion efficiency as in Belczynski+2008)
 ;    0.0 - assume all mass is lost from the system (fully nonconservative)
 ; Default: -1
-acc_lim = -1
+acc_lim = [-1, -1]
 
 ;;;;;;;;;;;;;
 ;;; Tides ;;;

src/cosmic/data/cosmic-settings.json

@@ -736,15 +736,15 @@
                 "settings-section-description": "Note: there are cases where a common envelope is forced regardless of the critical mass ratio for unstable mass transfer. In the following cases, a common envelope occurs regardless of the choices below:<ul><li><b>contact</b>: the stellar radii go into contact (common for similar ZAMS systems)</li><li><b>periapse contact</b>: the periapse distance is smaller than either of the stellar radii (common for highly eccentric systems)</li><li><b>core Roche overflow</b>: either of the stellar radii overflow their component's Roche radius (in this case, mass transfer from the convective core is always dynamically unstable)</li></ul>",
                 "name": "alpha1",
                 "description": "Common-envelope efficiency parameter which scales the efficiency of transferring orbital energy to the envelope. See <a href='https://ui.adsabs.harvard.edu/abs/2002MNRAS.329..897H/abstract'>Hurley+2002, Eq. 71</a>.",
-                "type": "number",
+                "type": "string",
                 "options-preface": "This value is equivalent to \\(\\alpha\\) in <a href='https://ui.adsabs.harvard.edu/abs/2002MNRAS.329..897H/abstract'>Hurley+2002, Eq. 71</a>.",
                 "options": [
                     {
                         "name": "positive values",
-                        "description": "Sets the common-envelope efficiency parameter, \\(\\alpha\\)"
+                        "description": "Sets the common-envelope efficiency parameter, \\(\\alpha\\), for each binary component"
                     },
                     {
-                        "name": 1.0,
+                        "name": "[1.0, 1.0]",
                         "description": "Default value",
                         "default": true
                     }
@@ -1627,14 +1627,13 @@
             },
             {
                 "name": "acc_lim",
-                "description": "Limits the amount of mass accreted during Roche-lobe overflow",
-                "type": "number",
+                "description": "Limits the amount of mass accreted during Roche-lobe overflow where each binary component must have its own limit specified in a list",
+                "type": "string",
                 "options-preface": "",
                 "options": [
                     {
                         "name": -1,
-                        "description": "limited to 10x the thermal rate of the accretor for MS/HG/CHeB and unlimited for GB/EAGB/AGB stars",
-                        "default": true
+                        "description": "limited to 10x the thermal rate of the accretor for MS/HG/CHeB and unlimited for GB/EAGB/AGB stars"
                     },
                     {
                         "name": -2,
@@ -1655,7 +1654,13 @@
                     {
                         "name": 0.0,
                         "description": "assume all mass is lost from the system (fully nonconservative)"
+                    },
+                    {
+                        "name": "[-1, -1]",
+                        "description": "Default choice",
+                        "default": true
                     }
+
                 ]
             },
             {

src/cosmic/evolve.py

@@ -124,7 +124,7 @@
 else:
     INITIAL_BINARY_TABLE_SAVE_COLUMNS = INITIAL_CONDITIONS_PASS_COLUMNS.copy()
 
-for col in ['natal_kick_array', 'qcrit_array', 'fprimc_array']:
+for col in ['natal_kick_array', 'qcrit_array', 'fprimc_array', 'alpha1', 'acc_lim']:
     INITIAL_BINARY_TABLE_SAVE_COLUMNS.remove(col)
 
 NATAL_KICK_COLUMNS = ['natal_kick',
@@ -140,10 +140,14 @@
 
 QCRIT_COLUMNS = ['qcrit_{0}'.format(kstar) for kstar in range(0, 16)]
 FPRIMC_COLUMNS = ['fprimc_{0}'.format(kstar) for kstar in range(0, 16)]
+ALPHA_COLUMNS = ['alpha1_{0}'.format(star) for star in range(0, 2)]
+ACCLIM_COLUMNS = ['acc_lim_{0}'.format(star) for star in range(0, 2)]
 
 INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(FLATTENED_NATAL_KICK_COLUMNS)
 INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(QCRIT_COLUMNS)
 INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(FPRIMC_COLUMNS)
+INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(ALPHA_COLUMNS)
+INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(ACCLIM_COLUMNS)
 
 # BSE doesn't need the binary fraction, so just add to columns for saving
 INITIAL_BINARY_TABLE_SAVE_COLUMNS.insert(7, 'binfrac')
@@ -344,6 +348,22 @@ def evolve(self, initialbinarytable, pool=None, bpp_columns=None, bcm_columns=No
                         initialbinarytable[col] = BSEDict['fprimc_array'][kstar]
                     else:
                         new_cols[col] = BSEDict['fprimc_array'][kstar]
+            elif k == 'alpha1':
+                columns_values = [BSEDict['alpha1']] * len(initialbinarytable)
+                initialbinarytable = initialbinarytable.assign(alpha1=columns_values)
+                for kstar in range(0,2):
+                    columns_values = pd.Series([BSEDict['alpha1'][kstar]] * len(initialbinarytable),
+                                               index=initialbinarytable.index,
+                                               name='alpha1_{0}'.format(kstar))
+                    initialbinarytable.loc[:, 'alpha1_{0}'.format(kstar)] = columns_values
+            elif k == 'acc_lim':
+                columns_values = [BSEDict['acc_lim']] * len(initialbinarytable)
+                initialbinarytable = initialbinarytable.assign(acc_lim=columns_values)
+                for kstar in range(0,2):
+                    columns_values = pd.Series([BSEDict['acc_lim'][kstar]] * len(initialbinarytable),
+                                               index=initialbinarytable.index,
+                                               name='acc_lim_{0}'.format(kstar))
+                    initialbinarytable.loc[:, 'acc_lim_{0}'.format(kstar)] = columns_values
             else:
                 # base case: if it's present, overwrite, if not, add to a list of new columns (see below)
                 if k in initialbinarytable.columns:
@@ -387,6 +407,12 @@ def evolve(self, initialbinarytable, pool=None, bpp_columns=None, bcm_columns=No
         if (pd.Series(FPRIMC_COLUMNS).isin(initialbinarytable.keys()).all()) and ('fprimc_array' not in BSEDict):
             initialbinarytable = initialbinarytable.assign(fprimc_array=initialbinarytable[FPRIMC_COLUMNS].values.tolist())
 
+        if (pd.Series(ALPHA_COLUMNS).isin(initialbinarytable.keys()).all()) and ('alpha1' not in BSEDict):
+            initialbinarytable = initialbinarytable.assign(alpha1=initialbinarytable[ALPHA_COLUMNS].values.tolist())
+
+        if (pd.Series(ACCLIM_COLUMNS).isin(initialbinarytable.keys()).all()) and ('acc_lim' not in BSEDict):
+            initialbinarytable = initialbinarytable.assign(acc_lim=initialbinarytable[ACCLIM_COLUMNS].values.tolist())
+
         # update timesteps based on mass modifier
         mass_modifier_cols = set(['pts1', 'pts2', 'pts3']).difference(columns_in_passed_initC)
         if dt_mass_modifiers and len(mass_modifier_cols) != 0:

src/cosmic/src/comenv.f

@@ -75,6 +75,7 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
       REAL*8 KW1_TEMP, KW2_TEMP
       REAL*8 rad(2),tms(2),lumin(2),B_0(2),bacc(2),tacc(2),epoch(2)
       REAL*8 menv_bpp(2),renv_bpp(2)
+      REAL*8 ALPHA_CE
 *
 * Initialize
 *
@@ -152,7 +153,12 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
 *
 * Calculate the final orbital energy without coalescence.
 *
-      EORBF = EORBI + EBINDI/ALPHA1
+      IF(switchedCE)THEN 
+         ALPHA_CE = ALPHA1(2)
+      ELSE
+         ALPHA_CE = ALPHA1(1)
+      ENDIF
+      EORBF = EORBI + EBINDI/ALPHA_CE
 *
 * If the secondary is on the main sequence see if it fills its Roche lobe.
 *
@@ -203,7 +209,12 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
 * Coalescence - calculate final binding energy.
 *
             EORBF = MAX(MC1*M2/(2.D0*SEPL),EORBI)
-            EBINDF = EBINDI - ALPHA1*(EORBF - EORBI)
+            IF(switchedCE)THEN 
+               ALPHA_CE = ALPHA1(2)
+            ELSE
+               ALPHA_CE = ALPHA1(1)
+            ENDIF
+            EBINDF = EBINDI - ALPHA_CE*(EORBF - EORBI)
             KW1_TEMP = KW
             KW2_TEMP = 15
          ELSE
@@ -441,7 +452,12 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
 * Calculate the final envelope binding energy.
 *
             EORBF = MAX(MC1*MC2/(2.D0*SEPL),EORBI)
-            EBINDF = EBINDI - ALPHA1*(EORBF - EORBI)
+            IF(switchedCE)THEN 
+               ALPHA_CE = ALPHA1(2)
+            ELSE
+               ALPHA_CE = ALPHA1(1)
+            ENDIF
+            EBINDF = EBINDI - ALPHA_CE*(EORBF - EORBI)
             if(output) write(*,*)'In dg or giant 1:',M01,M1,R1,M02,M2,
      & R2,MC1,MC2,MC3,KW1,KW2,KW,EORBF,EBINDF
 *

src/cosmic/src/const_bse.h

@@ -15,7 +15,7 @@
      &               qcflag,eddlimflag,bhspinflag,aic,rejuvflag,
      &               htpmb,ST_cr,ST_tide,bdecayfac,grflag,
      &               bhms_coll_flag,wd_mass_lim,rtmsflag,maltsev_mode
-      REAL*8 don_lim,acc_lim,Mbh_initial
+      REAL*8 don_lim,acc_lim(2),Mbh_initial
       COMMON /MTVARS/ don_lim,acc_lim,Mbh_initial
       INTEGER ceflag,cekickflag,cemergeflag,cehestarflag,ussn
       COMMON /CEFLAGS/ ceflag,cekickflag,cemergeflag,cehestarflag,ussn
@@ -29,7 +29,7 @@
       INTEGER LBV_flag
       COMMON /WINDVARS/ neta,bwind,hewind,beta,xi,acc2,epsnov,
      &                  eddfac,gamma,LBV_flag
-      REAL*8 alpha1,lambdaf
+      REAL*8 alpha1(2),lambdaf
       REAL*8 qcrit_array(16)
       COMMON /CEVARS/ qcrit_array,alpha1,lambdaf
       REAL*8 bconst,CK
@@ -61,7 +61,7 @@
       COMMON /POINTS/ pts1,pts2,pts3
       REAL*8 dmmax,drmax
       COMMON /TSTEPC/ dmmax,drmax
-      REAL*8 scm(50000,14),spp(20,3)
+      REAL*8 scm(50000,16),spp(25,20)
       COMMON /SINGLE/ scm,spp
       REAL*8 bcm(50000,52),bpp(1000,52)
       COMMON /BINARY/ bcm,bpp