Metisse integrate merge

.github/workflows/build_wheels_and_publish.yml

@@ -10,6 +10,7 @@ env:
   CIBW_ARCHS_LINUX: "x86_64"
   CIBW_SKIP: "*-win32 *musllinux*"
   CIBW_MANYLINUX_X86_64_IMAGE: manylinux2014
+  CIBW_ARCHS_MACOS: "x86_64 arm64"
   CIBW_BEFORE_BUILD: pip install numpy --config-settings=setup-args="-Dallow-noblas=true"
   CIBW_BUILD_VERBOSITY: "1"
   CIBW_BEFORE_ALL_LINUX: yum install -y gcc-gfortran
@@ -28,15 +29,16 @@ jobs:
     - uses: actions/checkout@v4
       with:
         fetch-depth: 0
-
-    # this step only uses a single python version since cibuildwheel
-    # builds multiple versions in one go on linux
+        submodules: recursive
+    - name: Switch METISSE submodule to track develop branch
+      run: |
+        git submodule set-branch --branch develop src/cosmic/src/METISSE
+        git submodule update --remote src/cosmic/src/METISSE
     - uses: actions/setup-python@v5
       name: Install Python
       with:
         python-version: "3.12"
-
-    - name: Install hdf5 (macOS)
+    - name: link gfortran and hdf5
       if: runner.os == 'macOS'
       run: |
         brew install hdf5

.github/workflows/python-package.yml

@@ -15,7 +15,15 @@ jobs:
         python-version: ['3.10', '3.11', '3.12', '3.13']
 
     steps:
-    - uses: actions/checkout@v4
+    - name: Checkout repository with submodules
+      uses: actions/checkout@v4
+      with:
+        submodules: recursive
+        fetch-depth: 0  # needed for full submodule history
+    - name: Switch METISSE submodule to track develop branch
+      run: |
+        git submodule set-branch --branch develop src/cosmic/src/METISSE
+        git submodule update --remote src/cosmic/src/METISSE
     - name: Set up Python ${{ matrix.python-version }}
       uses: actions/setup-python@v4
       with:

.gitmodules

@@ -0,0 +1,4 @@
+[submodule "src/cosmic/src/METISSE"]
+	path = src/cosmic/src/METISSE
+	url = https://github.com/TeamMETISSE/METISSE.git
+	branch = develop

bin/cosmic-pop

@@ -72,7 +72,7 @@ def parse_commandline():
 
     defaults = {}
     if not (args.inifile is None and (('-h' in remaining_argv) or ('--help' in remaining_argv))):
-        BSEDict, seed_int, filters, convergence, sampling = utils.parse_inifile(args.inifile)
+        BSEDict, SSEDict, seed_int, filters, convergence, sampling = utils.parse_inifile(args.inifile)
         defaults.update(sampling)
         defaults.update(filters)
         defaults.update(convergence)
@@ -169,8 +169,9 @@ if __name__ == '__main__':
     pool = schwimmbad.choose_pool(mpi=args.mpi, processes=args.nproc)
     if isinstance(pool, MPIPool):
         if not pool.is_master():
+            from sys import exit
             pool.wait()
-            sys.exit(0)
+            exit(0)
         nproc = len(pool.workers)
     else:
         nproc = args.nproc
@@ -181,7 +182,7 @@ if __name__ == '__main__':
 
     # READ AND PARSE INIFILE
     ###########################################################################
-    BSEDict, seed_int, filters, convergence, sampling = utils.parse_inifile(args.inifile)
+    BSEDict, SSEDict, seed_int, filters, convergence, sampling = utils.parse_inifile(args.inifile)
 
     # we now overwrite the inifile values with what was specified from the command line
     # (which could mean not overwriting anything at all because they are populated
@@ -217,8 +218,8 @@ if __name__ == '__main__':
                               "with {0}={2} from the commandline".format(argument, seed_int, getattr(args, argument)))
                 seed_int = getattr(args, argument)
 
-    # Check that the values in BSEDict, filters, and convergence are valid
-    utils.error_check(BSEDict, filters, convergence, sampling)
+    # Check that the values in BSEDict, SSEDict, filters, and convergence are valid
+    utils.error_check(BSEDict, SSEDict, filters, convergence, sampling)
 
     if seed_int != 0:
         np.random.seed(seed_int)
@@ -273,8 +274,8 @@ if __name__ == '__main__':
         # if config already exists, we need to overwrite it
         if "config" in f:
             del f["config"]
-        f["config"] = json.dumps({'BSEDict' : BSEDict, 'filters' : filters, 'convergence' : convergence,
-                                  'sampling' : sampling, 'rand_seed': seed_int})
+        f["config"] = json.dumps({'BSEDict' : BSEDict, 'SSEDict': SSEDict, 'filters' : filters, 
+                                  'convergence' : convergence, 'sampling' : sampling, 'rand_seed': seed_int})
         f.attrs['COSMIC_version'] = __version__
 
     # Initialize the step counter and convergence array/list
@@ -297,6 +298,9 @@ if __name__ == '__main__':
             log_file.write("You have specified both qmin and m2_min.\n")
             log_file.write("COSMIC will use qmin={} to determine the secondary masses in the initial sample.\n".format(args.qmin))
 
+    if sampling["sampling_method"] == "multidim" and args.Nstep < args.nproc ** 2:
+            warnings.warn("`Nstep` is less than the square of `nproc`. You will see repeats with multidim sampling. Consider raising `Nstep`.")
+
     while (Nstep < args.Niter) & (np.max(match) > convergence['match']) & ((time.time() - start_time) < args.max_wall_time):
         # Set random seed such that each iteration gets a unique, determinable seed
         rand_seed = seed_int + Nstep
@@ -322,8 +326,8 @@ if __name__ == '__main__':
                                                             met = sampling['metallicity'],
                                                             size = args.Nstep,
                                                             qmin = args.qmin,
-                                                            params = args.inifile)
-            # if m2_min is specified, use it to sample the initial binary table
+                                                            params = args.inifile,
+                                                            SSEDict = SSEDict)
             elif hasattr(args,'m2_min'):
                 init_samp_list = InitialBinaryTable.sampler(format_ = sampling['sampling_method'],
                                                             final_kstar1 = kstar1_range,
@@ -338,7 +342,8 @@ if __name__ == '__main__':
                                                             met = sampling['metallicity'],
                                                             size = args.Nstep,
                                                             m2_min = args.m2_min,
-                                                            params = args.inifile)
+                                                            params = args.inifile,
+                                                            SSEDict = SSEDict)
             else:
                 raise ValueError("You must specify either qmin or m2_min in the",
                                  " inifile if you are using the independent sampler")
@@ -389,6 +394,7 @@ if __name__ == '__main__':
         bpp, bcm, initCond, kick_info = Evolve.evolve(initialbinarytable=IBT,
                                                       pool=pool,
                                                       BSEDict=BSEDict,
+                                                      SSEDict=SSEDict,
                                                       idx=idx,
                                                       dtp=dtp,
                                                       timestep_conditions=filters['timestep_conditions'])

changelog.md

@@ -124,4 +124,4 @@ This release contains _several_ fixes to how CO core masses/remnant masses are h
 ## 3.7.6
 
 - Additions/changes:
-    - Added `pool` and `nproc` parameters to the `independent` sampler. Users can either specify an existing multiprocessing pool or a number of processors to use to create a pool that will be used for the sampling. This offers significant speed-ups for samples larger than 10,000 binaries.
+    - Added `pool` and `nproc` parameters to the `independent` sampler. Users can either specify an existing multiprocessing pool or a number of processors to use to create a pool that will be used for the sampling. This offers significant speed-ups for samples larger than 10,000 binaries.

ci/compile_benchmark.sh

@@ -1,2 +1,59 @@
-gfortran -coverage -fprofile-arcs -ftest-coverage -O0 src/cosmic/src/assign_remnant.f src/cosmic/src/benchmarkevolv2.f src/cosmic/src/comenv.f src/cosmic/src/corerd.f src/cosmic/src/deltat.f src/cosmic/src/dgcore.f src/cosmic/src/evolv2.f src/cosmic/src/gntage.f src/cosmic/src/hrdiag.f src/cosmic/src/hrdiag_remnant.f src/cosmic/src/instar.f src/cosmic/src/kick.f src/cosmic/src/mix.f src/cosmic/src/mlwind.f src/cosmic/src/mrenv.f src/cosmic/src/ran3.f src/cosmic/src/rl.f src/cosmic/src/star.f src/cosmic/src/zcnsts.f src/cosmic/src/zfuncs.f src/cosmic/src/concatkstars.f src/cosmic/src/bpp_array.f src/cosmic/src/checkstate.f  -o benchmarkevolv2.exe -I src/cosmic/src -Wl,-rpath,${CONDA_PREFIX}/lib
+#!/usr/bin/env bash
+set -euo pipefail
+
+# Directories
+SRC_DIR="src/cosmic/src"
+METISSE_DIR="$SRC_DIR/METISSE/src"
+
+# Compiler flags (removed rpath to CONDA_PREFIX)
+FFLAGS="-coverage -fprofile-arcs -ftest-coverage -O0 -J$SRC_DIR -I$SRC_DIR"
+
+# Phase 1: Compile METISSE modules in dependency order
+gfortran $FFLAGS -c \
+    $METISSE_DIR/track_support.f90 \
+    $METISSE_DIR/c_m_interface.f90 \
+    $METISSE_DIR/z_support.f90 \
+    $METISSE_DIR/sse_support.f90 \
+    $METISSE_DIR/remnant_support.f90 \
+    $METISSE_DIR/interp_support.f90 \
+    $METISSE_DIR/METISSE_gntage.f90 \
+    $METISSE_DIR/METISSE_deltat.f90 \
+    $METISSE_DIR/METISSE_mlwind.f90 \
+    $METISSE_DIR/METISSE_hrdiag.f90 \
+    $METISSE_DIR/METISSE_star.f90 \
+    $METISSE_DIR/METISSE_zcnsts.f90 \
+    $METISSE_DIR/comenv_lambda.f90 \
+    $METISSE_DIR/METISSE_miscellaneous.f90 \
+    $SRC_DIR/METISSE_utils.f90
+
+# Phase 2: Compile COSMIC and SSE sources + link everything
+gfortran $FFLAGS \
+    $SRC_DIR/hrdiag_remnant.f \
+    $SRC_DIR/assign_remnant.f \
+    $SRC_DIR/benchmarkevolv2.f \
+    $SRC_DIR/corerd.f \
+    $SRC_DIR/comenv.f \
+    $SRC_DIR/dgcore.f \
+    $SRC_DIR/evolv2.f \
+    $SRC_DIR/gntage.f \
+    $SRC_DIR/instar.f \
+    $SRC_DIR/kick.f \
+    $SRC_DIR/mix.f \
+    $SRC_DIR/mrenv.f \
+    $SRC_DIR/ran3.f \
+    $SRC_DIR/rl.f \
+    $SRC_DIR/concatkstars.f \
+    $SRC_DIR/comprad.f \
+    $SRC_DIR/bpp_array.f \
+    $SRC_DIR/checkstate.f \
+    $SRC_DIR/deltat.f \
+    $SRC_DIR/mlwind.f \
+    $SRC_DIR/hrdiag.f \
+    $SRC_DIR/star.f \
+    $SRC_DIR/zcnsts.f \
+    $SRC_DIR/SSE/SSE_*.f \
+    *.o \
+    -o benchmarkevolv2.exe
+
+# Run the benchmark
 ./benchmarkevolv2.exe

docs/generate_default_bsedict.py

@@ -28,7 +28,7 @@ def get_default_BSE_settings(to_python=False):
 
     if to_python:
         # ensure array settings are converted from strings to lists
-        for setting in ["qcrit_array", "natal_kick_array", "fprimc_array"]:
+        for setting in ["qcrit_array", "natal_kick_array", "alpha1", "acc_lim", "fprimc_array"]:
             # this one requires special handling because of the fractions
             if setting == "fprimc_array":
                 parts = defaults[setting].strip("[]").split(",")

docs/pages/about.rst

@@ -5,7 +5,13 @@ About
 #####
 COSMIC (Compact Object Synthesis and Monte Carlo Investigation Code) is a rapid binary population synthesis suite with a special focus of generating compact binary populations. 
 
-COSMIC currently implements stellar evolution using SSE (`Hurley, Pols, and Tout 2000 <https://ui.adsabs.harvard.edu/abs/2000MNRAS.315..543H/abstract>`_) and binary interactions using BSE (`Hurley, Tout, and Pols 2002 <https://ui.adsabs.harvard.edu/abs/2002MNRAS.329..897H/abstract>`_). Several modifications have been applied to BSE to account for recent updates to binary evolution especially important to compact binary formation (e.g. metallicity-dependent stellar winds or black hole natal kick strengths). For a detailed discussion of these modifications, see `Breivik et al. 2020 <https://ui.adsabs.harvard.edu/abs/2019arXiv191100903B/abstract>`_.
+COSMIC implements stellar evolution using the fitting formulae described in
+SSE (`Hurley, Pols, and Tout 2000 <https://ui.adsabs.harvard.edu/abs/2000MNRAS.315..543H/abstract>`_) 
+or METISSE which uses pre-computed detailed single star grids that follow the MIST-style equivalent 
+evolutionary phase (EEP) format. For more details on METISSE see the `METISSE documentation <https://cosmic-popsynth.github.io/Metisse/>`_. 
+
+Binary interactions are implemented based on the BSE algorithm (`Hurley, Tout, and Pols 2002 <https://ui.adsabs.harvard.edu/abs/2002MNRAS.329..897H/abstract>`_). 
+Several modifications have been applied to BSE to account for recent updates to binary evolution especially important to compact binary formation (e.g. metallicity-dependent stellar winds or black hole natal kick strengths). For a detailed discussion of these modifications, see `Breivik et al. 2020 <https://ui.adsabs.harvard.edu/abs/2019arXiv191100903B/abstract>`_.
 
 ************
 Using COSMIC

docs/pages/cite.rst

@@ -6,6 +6,9 @@ If you use COSMIC in your work, please cite both:
 
 - The COSMIC release paper: `Breivik et al. 2020 <https://ui.adsabs.harvard.edu/abs/2019arXiv191100903B/abstract>`_
 - The specific Zenodo DOI associated with version of COSMIC you used
+- The METISSE stellar evolution package: JOSS citation TBA
+
+If you used pre-supplied stellar tracks from METISSE, please also cite the relevant papers for those tracks as listed in the `METISSE documentation <https://cosmic-popsynth.github.io/Metisse/steller_models.html>`_.
 
 We recommend using the `Software Citation Station <https://tomwagg.com/software-citation-station>`_ to generate a citation
 for COSMIC (and any other software you use in your work!) - it's a great way to ensure you're citing software correctly and

docs/pages/evolve/evolve_sample.rst

@@ -45,4 +45,4 @@ guides:
 
     print(bcm.iloc[:10])
 
-    print(bpp)
+    print(bpp)

docs/pages/evolve/grid.rst

@@ -34,10 +34,12 @@ periods spaced evenly in log space.
         metallicity=np.ones(n_grid)*0.005
     )
 
+    SSEDict = {'stellar_engine': 'sse'}
+
     print(binary_grid)
 
-    bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=binary_grid, BSEDict=BSEDict)
+    bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=binary_grid, BSEDict=BSEDict, SSEDict=SSEDict)
 
     print(bpp)
 
-    print(bcm)
+    print(bcm)

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_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])

docs/pages/evolve/multiple.rst

@@ -12,7 +12,11 @@ Let's start by importing the necessary modules:
     from cosmic.evolve import Evolve
 
 
-And use the same BSE dict as before:
+And use the same SSE and BSE dictionaries as before:
+
+.. ipython:: python
+
+    SSEDict = {'stellar_engine': 'sse'}
 
 .. include:: ../../_generated/default_bsedict.rst
 
@@ -35,7 +39,7 @@ Below is an example for systems that could form GW150914 and GW170817 - like bin
 
     np.random.seed(5)
 
-    bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=binary_set, BSEDict=BSEDict)
+    bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=binary_set, BSEDict=BSEDict, SSEDict=SSEDict)
 
 As before, bpp, bcm, and initC are returned as pandas DataFrames which assign an
 index to each binary system we evolve. We can access each binary as follows:
@@ -59,7 +63,7 @@ progenitor, we expect most of the evolution to take place in the first ~60 Myr.
     :okexcept:
 
     from cosmic.plotting import evolve_and_plot
-    fig = evolve_and_plot(binary_set, t_min=None, t_max=[6.0, 60.0], BSEDict=BSEDict, sys_obs={})
+    fig = evolve_and_plot(binary_set, t_min=None, t_max=[6.0, 60.0], BSEDict=BSEDict, SSEDict=SSEDict, sys_obs={})
 
 
 .. plot::
@@ -69,5 +73,5 @@ progenitor, we expect most of the evolution to take place in the first ~60 Myr.
     import numpy as np
     np.random.seed(5)
     binary_set = InitialBinaryTable.InitialBinaries(m1=[85.543645, 11.171469], m2=[84.99784, 9.67305], porb=[446.795757, 370.758343], ecc=[0.448872, 0.370], tphysf=[13700.0, 13700.0], kstar1=[1, 1], kstar2=[1, 1], metallicity=[0.002, 0.02])
-    fig = evolve_and_plot(binary_set, t_min=None, t_max=[6.0, 60.0], BSEDict=default_BSEDict, sys_obs={})
+    fig = evolve_and_plot(binary_set, t_min=None, t_max=[6.0, 60.0], BSEDict=BSEDict, SSEDict=SSEDict, sys_obs={})
 

docs/pages/evolve/resolution.rst

@@ -40,9 +40,12 @@ First, print all time steps during mass transfer
         ecc=0.601408, tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.02
     )
 
+    SSEDict = {'stellar_engine': 'sse'}
 
 .. include:: ../../_generated/default_bsedict.rst
 
+    bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict,
+                                                SSEDict=SSEDict, timestep_conditions=timestep_conditions)
 
 .. ipython:: python
     :okwarning:
@@ -67,7 +70,8 @@ Second, pick a certain resolution for the bcm array until the system merges or i
     timestep_conditions = [['binstate=0', 'dtp=1.0']]
 
     bpp, bcm, initC, kick_info = Evolve.evolve(
-        initialbinarytable=single_binary, BSEDict=BSEDict,
+        initialbinarytable=single_binary, 
+        BSEDict=BSEDict, SSEDict=SSEDict, 
         timestep_conditions=timestep_conditions
     )
 
@@ -90,7 +94,8 @@ Finally, we show how to print a fine resolution only during the HMXB stage of th
     timestep_conditions = [['kstar_1=14', 'kstar_2<10','dtp=0.1'], ['kstar_2=14', 'kstar_1<10','dtp=0.1']]
 
     bpp, bcm, initC, kick_info = Evolve.evolve(
-        initialbinarytable=single_binary, BSEDict=BSEDict,
+        initialbinarytable=single_binary, 
+        BSEDict=BSEDict, SSEDict=SSEDict,
         timestep_conditions=timestep_conditions
     )
 

docs/pages/evolve/restart.rst

@@ -25,13 +25,15 @@ started from the beginning and three different points in the evolution:
         tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.002
     )
 
+    SSEDict = {'stellar_engine': 'sse'}
+
 .. include:: ../../_generated/default_bsedict.rst
 
 .. ipython:: python
     :okwarning:
 
     # evolve the binary
-    bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict)
+    bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict, SSEDict=SSEDict)
     print("From beginning")
     print(bpp)
 
@@ -82,7 +84,7 @@ restart the evolution after the mass transfer would occur. We can do this by usi
         tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.014*0.6
     )
     bpp, bcm, initC, kick_info = Evolve.evolve(
-        initialbinarytable=single_binary, BSEDict=BSEDict
+        initialbinarytable=single_binary, BSEDict=BSEDict, SSEDict=SSEDict
     )
 
     for column in bpp.columns:
@@ -98,7 +100,7 @@ restart the evolution after the mass transfer would occur. We can do this by usi
     initC['bin_num'] = np.linspace(0, 1000, 1000)
 
     bpp_restart, bcm_restart, initC_restart, kick_info_restart = Evolve.evolve(
-        initialbinarytable=initC, BSEDict={}
+        initialbinarytable=initC, BSEDict=BSEDict, SSEDict=SSEDict
     )
 
     bpp_BH = bpp_restart.loc[