COSMIC-PopSynth / COSMIC

COSMIC (Compact Object Synthesis and Monte Carlo Investigation Code)
GNU General Public License v3.0
45 stars 58 forks source link

[general issue] Divide by zero warning for `keep_singles=True` #604

Closed TomWagg closed 9 months ago

TomWagg commented 9 months ago

This warning appears

COSMIC/cosmic/utils.py:505 - RuntimeWarning: divide by zero encountered in divide
  q = M1 / M2

any time keep_singles=True for the sampler. This is likely because the "secondary" in the single has 0 mass when checking the Roche Lobe.

Unfortunately this was not fixed by making the single stars kstar=15

Here's a MWE

from cosmic.sample.initialbinarytable import InitialBinaryTable
from cosmic.evolve import Evolve

BSEDict = {'xi': 1.0, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'wdflag': 1, 'alpha1': 1.0,
           'pts1': 0.001, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 0.5,
           'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.0, 'mxns': 3.0, 'beta': -1.0, 'tflag': 1,
           'acc2': 1.5, 'grflag': 1, 'remnantflag': 4, 'ceflag': 0, 'eddfac': 1.0,
           'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -2.0, 'pisn': 45.0,
           'natal_kick_array': [[-100.0, -100.0, -100.0, -100.0, 0.0],
           [-100.0, -100.0, -100.0, -100.0, 0.0]], 'bhsigmafrac': 1.0,
           'polar_kick_angle': 90, '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],
           'cekickflag': 2, 'cehestarflag': 0, 'cemergeflag': 0, 'ecsn': 2.25,
           'ecsn_mlow': 1.6, 'aic': 1, 'ussn': 0, 'sigmadiv': -20.0, 'qcflag': 5,
           'eddlimflag': 0, 'fprimc_array': [2.0/21.0, 2.0/21.0, 2.0/21.0, 2.0/21.0,
                                             2.0/21.0, 2.0/21.0, 2.0/21.0, 2.0/21.0,
                                             2.0/21.0, 2.0/21.0, 2.0/21.0, 2.0/21.0,
                                             2.0/21.0, 2.0/21.0, 2.0/21.0, 2.0/21.0],
            'bhspinflag': 0, 'bhspinmag': 0.0, 'rejuv_fac': 1.0, 'rejuvflag': 0, 'htpmb': 1,
            'ST_cr': 1, 'ST_tide': 1, 'bdecayfac': 1, 'rembar_massloss': 0.5, 'kickflag': 0,
            'zsun': 0.014, 'bhms_coll_flag': 0, 'don_lim': -1, 'acc_lim': -1, 'binfrac': 0.5,
            'rtmsflag': 0}

initial_binaries = InitialBinaryTable.sampler('independent', np.linspace(0, 15, 16), np.linspace(0, 15, 16),
                                              binfrac_model=0.5, SF_start=10.0,
                                              SF_duration=0.0, met=0.02, size=10,
                                              primary_model='kroupa01', ecc_model='sana12', porb_model='sana12',
                                              keep_singles=True)[0]
Evolve.evolve(initialbinarytable=initial_binaries, BSEDict=BSEDict)