diff --git a/compas_python_utils/preprocessing/compasConfigDefault.yaml b/compas_python_utils/preprocessing/compasConfigDefault.yaml index f01e21860..607817553 100644 --- a/compas_python_utils/preprocessing/compasConfigDefault.yaml +++ b/compas_python_utils/preprocessing/compasConfigDefault.yaml @@ -1,5 +1,5 @@ ##~!!~## COMPAS option values -##~!!~## File Created Mon Feb 2 10:03:24 2026 by COMPAS v03.29.00 +##~!!~## File Created Sun Nov 1 20:57:27 2026 by COMPAS v03.31.00 ##~!!~## ##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has ##~!!~## all COMPAS option entries commented so that the COMPAS default value for the @@ -53,6 +53,9 @@ booleanChoices: # --common-envelope-lambda-nanjing-use-rejuvenated-mass: False # Default: False # --revised-energy-formalism-nandez-ivanova: False # Default: False + ### CIRCUMBINARY DISK +# --circumbinary-disk: False # Default: False + ### SUPERNOVAE, KICKS AND REMNANTS # --allow-non-stripped-ECSN: False # Default: False # --pair-instability-supernovae: True # Default: True @@ -169,6 +172,17 @@ numericalChoices: # --common-envelope-second-stage-beta: 0.000000 # Default: 0.000000 # --common-envelope-slope-kruckow: -0.833333 # Default: -0.833333 # --maximum-mass-donor-nandez-ivanova: 2.000000 # Default: 2.000000 +# --post-common-envelope-eccentricity-cap: 0.100000 # Default: 0.100000 +# --post-common-envelope-eccentricity-fraction: 0.100000 # Default: 0.100000 + + ### CIRCUMBINARY DISK +# --circumbinary-disk-beta: 3.000000 # Default: 3.000000 +# --circumbinary-disk-eddington-cap-factor: 1.000000 # Default: 1.000000 +# --circumbinary-disk-eddington-ge23-factor: 0.666667 # Default: 0.666667 +# --circumbinary-disk-envelope-mass-fraction-supplied: 0.050000 # Default: 0.050000 +# --circumbinary-disk-inner-radius-over-separation: 2.500000 # Default: 2.500000 +# --circumbinary-disk-lifetime: 100000.00 # Default: 100000.00 +# --circumbinary-disk-structure-factor: 4.000000 # Default: 4.000000 ### SUPERNOVAE, KICKS AND REMNANTS # --eddington-accretion-factor: 1.000000 # Default: 1.000000 # multiplication Factor for eddington accretion onto NS&BH @@ -264,7 +278,7 @@ stringChoices: # --case-BB-stability-prescription: 'ALWAYS_STABLE' # Default: 'ALWAYS_STABLE' # Options: ['ALWAYS_UNSTABLE','TREAT_AS_OTHER_MT','ALWAYS_STABLE_ONTO_NSBH','ALWAYS_STABLE'] # --critical-mass-ratio-prescription: 'NONE' # Default: 'NONE' # Options: ['HURLEY_HJELLMING_WEBBINK','GE','GE_IC','CLAEYS','NONE'] # --stellar-zeta-prescription: 'SOBERMAN' # Default: 'SOBERMAN' # Options: ['ARBITRARY','HURLEY','SOBERMAN'] -# --mass-transfer-angular-momentum-loss-prescription: 'ISOTROPIC' # Default: 'ISOTROPIC' # Options: ['ARBITRARY','KLENCKI_LINEAR','MACLEOD_LINEAR','CIRCUMBINARY','ISOTROPIC','JEANS'] +# --mass-transfer-angular-momentum-loss-prescription: 'ISOTROPIC' # Default: 'ISOTROPIC' # Options: ['MIN','ARBITRARY','KLENCKI_LINEAR','MACLEOD_LINEAR','CIRCUMBINARY','ISOTROPIC','JEANS'] # --mass-transfer-accretion-efficiency-prescription: 'THERMAL' # Default: 'THERMAL' # Options: ['HAMSTARS','FIXED','THERMAL'] # --mass-transfer-rejuvenation-prescription: 'STARTRACK' # Default: 'STARTRACK' # Options: ['STARTRACK','HURLEY'] # --mass-transfer-thermal-limit-accretor-multiplier: 'CFACTOR' # Default: 'CFACTOR' # Options: ['ROCHELOBE','CFACTOR'] @@ -274,7 +288,11 @@ stringChoices: # --common-envelope-formalism: 'ENERGY' # Default: 'ENERGY' # Options: ['TWO_STAGE','ENERGY'] # --common-envelope-lambda-prescription: 'LAMBDA_NANJING' # Default: 'LAMBDA_NANJING' # Options: ['LAMBDA_DEWI','LAMBDA_KRUCKOW','LAMBDA_NANJING','LAMBDA_LOVERIDGE','LAMBDA_FIXED'] # Xu & Li 2010 # --common-envelope-mass-accretion-prescription: 'ZERO' # Default: 'ZERO' # Options: ['CHEVALIER','MACLEOD','UNIFORM','CONSTANT','ZERO'] -# --common-envelope-second-stage-gamma-prescription: 'ISOTROPIC' # Default: 'ISOTROPIC' # Options: ['ARBITRARY','KLENCKI_LINEAR','MACLEOD_LINEAR','CIRCUMBINARY','ISOTROPIC','JEANS'] +# --common-envelope-second-stage-gamma-prescription: 'ISOTROPIC' # Default: 'ISOTROPIC' # Options: ['MIN','ARBITRARY','KLENCKI_LINEAR','MACLEOD_LINEAR','CIRCUMBINARY','ISOTROPIC','JEANS'] +# --post-common-envelope-eccentricity-prescription: 'CIRCULAR' # Default: 'CIRCULAR' # Options: ['PRE_RLOF_FRACTION','PRE_RLOF_CAP','CIRCULAR'] + + ### CIRCUMBINARY DISK +# --circumbinary-disk-eddington-mode: 'GE23' # Default: 'GE23' # Options: ['TM23','GE23','CAP'] ### TIDES # --tides-prescription: 'NONE' # Default: 'NONE' # Options: ['ZAHN1977','KAPIL2026','PERFECT','NONE'] diff --git a/online-docs/pages/Getting started/getting-started.rst b/online-docs/pages/Getting started/getting-started.rst index f2cc6de8f..b59e9e203 100644 --- a/online-docs/pages/Getting started/getting-started.rst +++ b/online-docs/pages/Getting started/getting-started.rst @@ -11,7 +11,7 @@ To start using COMPAS, `git clone` the code, install dependencies and build COMP Install dependencies -------------------- -COMPAS requires a ``C++`` compiler, and the ``gsl``, ``boost``, and ``hdf5`` libraries. ``Python`` is required for the COMPAS post-processing tools. +COMPAS requires a ``C++`` compiler, and the ``gsl``, ``boost``, and ``hdf5`` libraries. ``Python`` is required for the COMPAS post-processing tools. Use of the circumbinary disk module requires installing `spindler-c `_. .. tabs:: diff --git a/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-BSE-circumbinary-disks.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-BSE-circumbinary-disks.rst new file mode 100644 index 000000000..31a2fed0a --- /dev/null +++ b/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-BSE-circumbinary-disks.rst @@ -0,0 +1,21 @@ +BSE circumbinary disks +====================== + +Default record definition for the BSE Circumbinary Disks log file:: + + const ANY_PROPERTY_VECTOR BSE_CIRCUMBINARY_DISKS_REC = { + BINARY_PROPERTY::RANDOM_SEED, + BINARY_PROPERTY::CIRCUMBINARY_DISK_FORMATION_TIME, + BINARY_PROPERTY::CIRCUMBINARY_DISK_INITIAL_MASS, + BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_SUPPLIED, + BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_RETAINED_1, + BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_RETAINED_2, + BINARY_PROPERTY::CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_PRE, + BINARY_PROPERTY::CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_POST, + BINARY_PROPERTY::CIRCUMBINARY_DISK_ECCENTRICITY_PRE, + BINARY_PROPERTY::CIRCUMBINARY_DISK_ECCENTRICITY_POST, + BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_1_AT_FORMATION, + BINARY_PROPERTY::CIRCUMBINARY_DISK_STELLAR_TYPE_1_AT_FORMATION, + BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_2_AT_FORMATION, + BINARY_PROPERTY::CIRCUMBINARY_DISK_STELLAR_TYPE_2_AT_FORMATION + }; diff --git a/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications.rst index 786733766..fb3d9566d 100644 --- a/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications.rst +++ b/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications.rst @@ -27,6 +27,7 @@ For Binary Star Evolution (BSE): standard-logfiles-default-record-specifications-BSE-supernovae standard-logfiles-default-record-specifications-BSE-dco standard-logfiles-default-record-specifications-BSE-ce + standard-logfiles-default-record-specifications-BSE-circumbinary-disks standard-logfiles-default-record-specifications-BSE-pulsars standard-logfiles-default-record-specifications-BSE-rlof standard-logfiles-default-record-specifications-BSE-switchlog diff --git a/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-binary.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-binary.rst index bd7cf724a..851c67241 100644 --- a/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-binary.rst +++ b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-binary.rst @@ -49,6 +49,240 @@ Binary Properties * - Header String: - Tau_Circ + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_FORMED** + - + * - Data type: + - BOOL + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.formed + * - Description: + - Flag to indicate if a post-common-envelope circumbinary disk formed and was evolved. + * - Header String: + - CBD + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_MASS_SUPPLIED** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.massSupplied + * - Description: + - Mass supplied to the binary through the circumbinary disk, potentially less than the initial reservoir if interrupted (\ :math:`M_\odot`). + * - Header String: + - CBD_Mass_Supplied + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_MASS_RETAINED_1** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.massRetained1 + * - Description: + - Mass retained by star 1 during circumbinary disk accretion (\ :math:`M_\odot`). + * - Header String: + - CBD_Mass_Retained(1) + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_MASS_RETAINED_2** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.massRetained2 + * - Description: + - Mass retained by star 2 during circumbinary disk accretion (\ :math:`M_\odot`). + * - Header String: + - CBD_Mass_Retained(2) + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_PRE** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.semiMajorAxisPre + * - Description: + - Semi-major axis before circumbinary disk evolution (\ :math:`R_\odot`). + * - Header String: + - SemiMajorAxisCBD + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_ECCENTRICITY_PRE** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.eccentricityPre + * - Description: + - Eccentricity before circumbinary disk evolution. + * - Header String: + - EccentricityCBD + + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_FORMATION_TIME** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.formationTime + * - Description: + - Time at which the disk formed (Myr). + * - Header String: + - Time + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_INITIAL_MASS** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.initialReservoirMass + * - Description: + - Initial mass assigned to the CBD reservoir (\ :math:`M_\odot`). + * - Header String: + - Initial_CBD_Mass + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_MASS_1_AT_FORMATION** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.mass1AtFormation + * - Description: + - Mass of star 1 when the CBD formed (\ :math:`M_\odot`). + * - Header String: + - Mass(1) + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_MASS_2_AT_FORMATION** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.mass2AtFormation + * - Description: + - Mass of star 2 when the CBD formed (\ :math:`M_\odot`). + * - Header String: + - Mass(2) + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_STELLAR_TYPE_1_AT_FORMATION** + - + * - Data type: + - STELLAR_TYPE + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.stellarType1AtFormation + * - Description: + - Stellar type of star 1 when the CBD formed. + * - Header String: + - Stellar_Type(1) + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_STELLAR_TYPE_2_AT_FORMATION** + - + * - Data type: + - STELLAR_TYPE + * - COMPAS variable: + - BaseBinaryStar::m_CircumbinaryDiskDetails.stellarType2AtFormation + * - Description: + - Stellar type of star 2 when the CBD formed. + * - Header String: + - Stellar_Type(2) + .. flat-table:: :widths: 25 75 1 1 :header-rows: 0 diff --git a/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-options.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-options.rst index cdb25945a..c11c29f7d 100644 --- a/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-options.rst +++ b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-options.rst @@ -194,6 +194,163 @@ DEPRECATION NOTICE: property ``BLACK_HOLE_KICKS`` has been deprecated and will s * - Header String: - PO_Circularise@\ MT +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK** + - + * - Data type: + - BOOL + * - COMPAS variable: + - Options::m_CircumbinaryDisk + * - Description: + - Value of program option ``--circumbinary-disk`` + * - Header String: + - PO_CBD + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_BETA** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - Options::m_CircumbinaryDiskBeta + * - Description: + - Value of program option ``--circumbinary-disk-beta`` + * - Header String: + - PO_CBD_Beta + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_ENVELOPE_MASS_FRACTION_SUPPLIED** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - Options::m_CircumbinaryDiskEnvelopeMassFractionSupplied + * - Description: + - Value of program option ``--circumbinary-disk-envelope-mass-fraction-supplied`` + * - Header String: + - PO_CBD_Envelope_Mass_Fraction_Supplied + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_INNER_RADIUS_OVER_SEPARATION** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - Options::m_CircumbinaryDiskInnerRadiusOverSeparation + * - Description: + - Value of program option ``--circumbinary-disk-inner-radius-over-separation`` + * - Header String: + - PO_CBD_Inner_Radius_Over_Separation + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_LIFETIME** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - Options::m_CircumbinaryDiskLifetime + * - Description: + - Value of program option ``--circumbinary-disk-lifetime`` + * - Header String: + - PO_CBD_Lifetime + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_STRUCTURE_FACTOR** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - Options::m_CircumbinaryDiskStructureFactor + * - Description: + - Value of program option ``--circumbinary-disk-structure-factor`` + * - Header String: + - PO_CBD_Structure_Factor + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_EDDINGTON_MODE** + - + * - Data type: + - INT + * - COMPAS variable: + - Options::m_CircumbinaryDiskEddingtonMode + * - Description: + - Value of program option ``--circumbinary-disk-eddington-mode`` + * - Header String: + - PO_CBD_Eddington_Mode + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_EDDINGTON_CAP_FACTOR** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - Options::m_CircumbinaryDiskEddingtonCapFactor + * - Description: + - Value of program option ``--circumbinary-disk-eddington-cap-factor`` + * - Header String: + - PO_CBD_Eddington_Cap_Factor + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **CIRCUMBINARY_DISK_EDDINGTON_GE23_FACTOR** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - Options::m_CircumbinaryDiskEddingtonGe23Factor + * - Description: + - Value of program option ``--circumbinary-disk-eddington-ge23-factor`` + * - Header String: + - PO_CBD_Eddington_GE23_Factor + + + + + + .. flat-table:: :widths: 25 75 1 1 :header-rows: 0 @@ -1620,6 +1777,57 @@ DEPRECATION NOTICE: property ``LBV_PRESCRIPTION`` has been deprecated and will s * - Header String: - PO_PISN_Upper_Limit +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION** + - + * - Data type: + - INT + * - COMPAS variable: + - Options::m_PostCommonEnvelopeEccentricityPrescription + * - Description: + - Value of program option ``--post-common-envelope-eccentricity-prescription`` + * - Header String: + - PO_Post_CE_Eccentricity_Prescription + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **POST_COMMON_ENVELOPE_ECCENTRICITY_CAP** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - Options::m_PostCommonEnvelopeEccentricityCap + * - Description: + - Value of program option ``--post-common-envelope-eccentricity-cap`` + * - Header String: + - PO_Post_CE_Eccentricity_Cap + + +.. flat-table:: + :widths: 25 75 1 1 + :header-rows: 0 + :class: aligned-text + + * - :cspan:`2` **POST_COMMON_ENVELOPE_ECCENTRICITY_FRACTION** + - + * - Data type: + - DOUBLE + * - COMPAS variable: + - Options::m_PostCommonEnvelopeEccentricityFraction + * - Description: + - Value of program option ``--post-common-envelope-eccentricity-fraction`` + * - Header String: + - PO_Post_CE_Eccentricity_Fraction + + .. flat-table:: :widths: 25 75 1 1 :header-rows: 0 diff --git a/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification.rst index d5742f7b0..dfb59d686 100644 --- a/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification.rst +++ b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification.rst @@ -50,6 +50,7 @@ Log file definitions file specification "BSE_DCO_REC" | # BSE only "BSE_SNE_REC" | # BSE only "BSE_CEE_REC" | # BSE only + "BSE_CBD_REC" | # BSE only "BSE_PULSARS_REC" | # BSE only "BSE_RLOF_REC" | # BSE only "BSE_DETAILED_REC" | # BSE only diff --git a/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-types-bse-circumbinary-disks.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-types-bse-circumbinary-disks.rst new file mode 100644 index 000000000..db1538bed --- /dev/null +++ b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-types-bse-circumbinary-disks.rst @@ -0,0 +1,7 @@ +BSE circumbinary disks file record types +======================================== + +Following is a list of the BSE Circumbinary Disks file record type numbers and corresponding symbolic names, and their meaning: + +1. DEFAULT |BR| + Default BSE_CIRCUMBINARY_DISKS file record type. \ No newline at end of file diff --git a/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-types.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-types.rst index 249e2d7b9..39658693a 100644 --- a/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-types.rst +++ b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-types.rst @@ -35,6 +35,7 @@ Each standard log file has its own set of record types - select a file below to standard-logfiles-record-types-bse-pulsar-evolution standard-logfiles-record-types-bse-double-compact-objects standard-logfiles-record-types-bse-common-envelopes + standard-logfiles-record-types-bse-circumbinary-disks standard-logfiles-record-types-bse-detailed-output standard-logfiles-record-types-bse-system-snapshot diff --git a/online-docs/pages/User guide/COMPAS output/standard-logfiles.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles.rst index e04efdbda..f3e616675 100644 --- a/online-docs/pages/User guide/COMPAS output/standard-logfiles.rst +++ b/online-docs/pages/User guide/COMPAS output/standard-logfiles.rst @@ -45,6 +45,12 @@ The standard log files are: - Records summary information for all binary systems that experience CEEs during BSE. * - - + * - Circumbinary Disks + - Records summary information for each post-common-envelope circumbinary disk that forms during BSE. + * - + - Enable with program option ``--circumbinary-disk``. + * - + - * - Pulsar Evolution - Records detailed Pulsar evolution information during SSE/BSE. * - diff --git a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst index c86d69f76..a542d63f4 100644 --- a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst +++ b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst @@ -124,6 +124,46 @@ Default = PESSIMISTIC |br| Circularise binary when it enters a Mass Transfer episode. |br| Default = TRUE +**--circumbinary-disk** |br| +Enable post-common-envelope circumbinary disk formation and evolution. +|br| Default = FALSE + +**--circumbinary-disk-beta** |br| +Ratio between the fraction of the binary angular momentum and the fraction of the common-envelope mass available for circumbinary disk formation. +|br| Default = 3.0 + +**--circumbinary-disk-envelope-mass-fraction-supplied** |br| +Fraction of the common-envelope mass supplied to the post-common-envelope binary through a circumbinary disk. +|br| Default = 0.05 + +**--circumbinary-disk-inner-radius-over-separation** |br| +Inner radius of the circumbinary disk in units of the post-common-envelope binary separation. +|br| Default = 2.5 + +**--circumbinary-disk-lifetime** |br| +Duration, in yr, over which the supplied circumbinary disk mass is delivered to the binary. +|br| Default = :math:`10^{5}` + +**--circumbinary-disk-structure-factor** |br| +Disk structure factor multiplying the circumbinary disk angular-momentum threshold. +|br| Default = 4.0 + +**--circumbinary-disk-eddington-mode** |br| +Eddington-limited accretion prescription used during circumbinary disk evolution. +|br| Options: { CAP, GE23, TM23 } +|br| ``CAP`` limits the accretion to ``--circumbinary-disk-eddington-cap-factor`` times the Eddington accretion rate. +|br| ``GE23`` reduces the super-Eddington retained accretion rate by a factor of ``--circumbinary-disk-eddington-ge23-factor``, following Ghodla & Eldridge 2023. +|br| ``TM23`` applies the Tuna & Metzger 2023 compact-object accretion prescription. +|br| Default = GE23 + +**--circumbinary-disk-eddington-cap-factor** |br| +Multiplier applied to the Eddington accretion rate in ``CAP`` mode. +|br| Default = 1.0 + +**--circumbinary-disk-eddington-ge23-factor** |br| +Retained fraction of the requested super-Eddington accretion rate in ``GE23`` mode. +|br| Default = :math:`2 / 3` + **--common-envelope-allow-immediate-RLOF-post-CE-survive** |br| Allow binaries that experience Roche lobe overflow immediately at the end of the CE phase to survive. |br| Default = FALSE @@ -207,7 +247,7 @@ Default = ZERO Recombination energy density (erg gāˆ’1). |br| Default = :math:`1.5 \times 10^{13}` -**common-envelope-second-stage-beta** |br| +**--common-envelope-second-stage-beta** |br| The fraction of the mass lost during the second stage of a 2-stage common envelope that is accreted. |br| Default = 0 @@ -1136,6 +1176,24 @@ Default = 135.0 Print details of population. |br| Default = FALSE +**--post-common-envelope-eccentricity-prescription** |br| +Prescription for setting the eccentricity immediately after common-envelope evolution. +|br| Options: { CIRCULAR, PRE_RLOF_CAP, PRE_RLOF_FRACTION } +|br| ``CIRCULAR`` assigns :math:`e_{\rm postCE} = 0`. +|br| ``PRE_RLOF_CAP`` assigns :math:`e_{\rm postCE} = \min(e_{\rm preRLOF}, e_{\rm CE,max})`, where :math:`e_{\rm CE,max}` is set by ``--post-common-envelope-eccentricity-cap``. +|br| ``PRE_RLOF_FRACTION`` assigns :math:`e_{\rm postCE} = f_e e_{\rm preRLOF}`, where :math:`f_e` is set by ``--post-common-envelope-eccentricity-fraction``. +|br| Default = CIRCULAR + +**--post-common-envelope-eccentricity-cap** |br| +Upper cap on the post-common-envelope eccentricity when ``--post-common-envelope-eccentricity-prescription`` is ``PRE_RLOF_CAP``. +|br| In this mode, :math:`e_{\rm postCE} = \min(e_{\rm preRLOF}, e_{\rm CE,max})`. +|br| Default = 0.1 + +**--post-common-envelope-eccentricity-fraction** |br| +Fraction of the pre-RLOF eccentricity retained after common-envelope evolution when ``--post-common-envelope-eccentricity-prescription`` is ``PRE_RLOF_FRACTION``. +|br| In this mode, :math:`e_{\rm postCE} = f_e e_{\rm preRLOF}`. +|br| Default = 0.1 + **--PPI-CO-Core-Shift-Hendriks** |br| Shift in CO core mass for PPI (in Msol) for the Hendriks+23 PPI prescription Default = 0.0 @@ -1580,7 +1638,11 @@ Go to :ref:`the top of this page ` for the full alphabetical --common-envelope-lambda-nanjing-use_rejuvenated-mass, --common-envelope-allow-main-sequence-survive, --common-envelope-allow-radiative-envelope-survive, --common-envelope-allow-immediate-RLOF-post-CE-survive, --common-envelope-mass-accretion-prescription, --common-envelope-mass-accretion-constant, --common-envelope-mass-accretion-min, --common-envelope-mass-accretion-max, --common-envelope-recombination-energy-density, --maximum-mass-donor-nandez-ivanova, ---revised-energy-formalism-nandez-ivanova, common-envelope-second-stage-beta, common-envelope-second-stage-gamma-prescription +--revised-energy-formalism-nandez-ivanova, --common-envelope-second-stage-beta, --common-envelope-second-stage-gamma-prescription, --post-common-envelope-eccentricity-prescription, --post-common-envelope-eccentricity-cap, --post-common-envelope-eccentricity-fraction + +--circumbinary-disk, --circumbinary-disk-beta, --circumbinary-disk-envelope-mass-fraction-supplied, --circumbinary-disk-inner-radius-over-separation, --circumbinary-disk-lifetime, +--circumbinary-disk-structure-factor, --circumbinary-disk-eddington-mode, --circumbinary-disk-eddington-cap-factor, --circumbinary-disk-eddington-ge23-factor + :ref:`Back to Top ` diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 04887c1c5..30ca591bd 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -3,6 +3,10 @@ What's new Following is a brief list of important updates to the COMPAS code. A complete record of changes can be found in the file ``changelog.h``. +**03.31.00 November 2, 2026** + +* Added post-common-envelope circumbinary disk formation and evolution following Grichener et al. (2026), including new CBD program options and the ``BSE_Circumbinary_Disks`` logfile. + **03.29.00 January 30, 2026** * Added new tidal prescription based on Zahn (1977) and Hurley et. al (2002), called ``ZAHN1977``. * Renamed the ``KAPIL2025`` tides prescription to ``KAPIL2026`` to match paper date. diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 19b575eb6..307fb7b30 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1,11 +1,19 @@ #include "BaseBinaryStar.h" +#include "CircumbinaryDisk.h" #include +#include #include "vector3d.h" + // gsl includes #include +// CBD reservoir bookkeeping tolerances used only by BaseBinaryStar. +namespace { +constexpr double CBD_RESERVOIR_COMPLETION_ABSOLUTE_TOLERANCE_MSOL = 1.0E-12; +constexpr double CBD_RESERVOIR_COMPLETION_RELATIVE_TOLERANCE = 1.0E-10; +} /* Constructor * @@ -268,6 +276,8 @@ void BaseBinaryStar::SetRemainingValues() { // Initialise other parameters m_SemiMajorAxisPrev = m_SemiMajorAxis; m_EccentricityPrev = m_Eccentricity; + m_SemiMajorAxisPreRLOF = m_SemiMajorAxis; + m_EccentricityPreRLOF = m_Eccentricity; // initial binary parameters - kept constant as a record of the initial parameters of the binary m_SemiMajorAxisInitial = m_SemiMajorAxis; @@ -282,12 +292,15 @@ void BaseBinaryStar::SetRemainingValues() { m_SemiMajorAxisAtDCOFormation = DEFAULT_INITIAL_DOUBLE_VALUE; m_EccentricityAtDCOFormation = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RestartTimestepAfterImmediateEvent = false; + double momentOfInertia1 = m_Star1->CalculateMomentOfInertiaAU(); double momentOfInertia2 = m_Star2->CalculateMomentOfInertiaAU(); m_TotalEnergy = CalculateTotalEnergy(m_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Omega(), m_Star2->Omega(), momentOfInertia1, momentOfInertia2); m_TotalAngularMomentum = CalculateAngularMomentum(m_SemiMajorAxis, m_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Omega(), m_Star2->Omega(), momentOfInertia1, momentOfInertia2); + m_TotalAngularMomentumPreRLOF = m_TotalAngularMomentum; m_TotalAngularMomentumPrev = m_TotalAngularMomentum; if (OPTIONS->CHEMode() != CHE_MODE::NONE) { // CHE enabled? @@ -524,6 +537,20 @@ COMPAS_VARIABLE BaseBinaryStar::BinaryPropertyValue(const T_ANY_PROPERTY p_Prope switch (property) { // which property? case BINARY_PROPERTY::CIRCULARIZATION_TIMESCALE: value = CircularizationTimescale(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_FORMED: value = CircumbinaryDiskFormed(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_FORMATION_TIME: value = CircumbinaryDiskFormationTime(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_1_AT_FORMATION: value = CircumbinaryDiskMass1AtFormation(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_2_AT_FORMATION: value = CircumbinaryDiskMass2AtFormation(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_STELLAR_TYPE_1_AT_FORMATION: value = CircumbinaryDiskStellarType1AtFormation(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_STELLAR_TYPE_2_AT_FORMATION: value = CircumbinaryDiskStellarType2AtFormation(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_INITIAL_MASS: value = CircumbinaryDiskInitialMass(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_RETAINED_1: value = CircumbinaryDiskMassRetained1(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_RETAINED_2: value = CircumbinaryDiskMassRetained2(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_SUPPLIED: value = CircumbinaryDiskMassSupplied(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_PRE: value = CircumbinaryDiskSemiMajorAxisPre(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_POST: value = CircumbinaryDiskSemiMajorAxisPost(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_ECCENTRICITY_PRE: value = CircumbinaryDiskEccentricityPre(); break; + case BINARY_PROPERTY::CIRCUMBINARY_DISK_ECCENTRICITY_POST: value = CircumbinaryDiskEccentricityPost(); break; case BINARY_PROPERTY::COMMON_ENVELOPE_AT_LEAST_ONCE: value = CEAtLeastOnce(); break; case BINARY_PROPERTY::COMMON_ENVELOPE_EVENT_COUNT: value = CommonEnvelopeEventCount(); break; case BINARY_PROPERTY::UNBOUND: value = Unbound(); break; @@ -1527,6 +1554,321 @@ void BaseBinaryStar::EvaluateSupernovae() { } +/* + * Apply optional post-CE circumbinary disk evolution, following the steps below: + * + * 1. Disk formation. The CircumbinaryDisk module assigns a fraction of the total + * binary mass lost during CE and a fraction of the binary's pre-RLOF + * total angular momentum to the disk, then checks whether that angular + * momentum is large enough for the material to orbit outside the adopted + * inner disk radius. If not, COMPAS continues with the ordinary post-CE state. + * + * 2. Disk evolution. The CBD orbital evolution in response to accreted mass + * is calculated using spindler-c (based on Valli+2024), modified to account for Eddington-limited accretion. + * Mass that is not accreted due to the Eddington limit is assumed to be ejected + * from the vicinity of the accretor, carrying its component-specific orbital + * angular momentum. The Spindler formalism assumes that the initial CBD + * reservoir is much smaller than the post-CE binary mass. + * + * A formed CBD is stored as an active reservoir. Later COMPAS timesteps + * process it before ordinary binary event checks, and the CBD event is logged + * when the reservoir is exhausted or the disk is truncated by another event. + * + * void ApplyCircumbinaryDiskAfterCommonEnvelope(const STELLAR_TYPE p_StellarType1PreCE, + * const STELLAR_TYPE p_StellarType2PreCE, + * const double p_Mass1PreCE, + * const double p_Mass2PreCE, + * const bool p_EnvelopeFlag1, + * const bool p_EnvelopeFlag2) + * + * @param [IN] p_StellarType1PreCE Stellar type of star 1 immediately before CE + * @param [IN] p_StellarType2PreCE Stellar type of star 2 immediately before CE + * @param [IN] p_Mass1PreCE Mass of star 1 immediately before CE (Msol) + * @param [IN] p_Mass2PreCE Mass of star 2 immediately before CE (Msol) + * @param [IN] p_EnvelopeFlag1 Whether star 1 had an envelope immediately before CE + * @param [IN] p_EnvelopeFlag2 Whether star 2 had an envelope immediately before CE + */ +void BaseBinaryStar::ApplyCircumbinaryDiskAfterCommonEnvelope(const STELLAR_TYPE p_StellarType1PreCE, + const STELLAR_TYPE p_StellarType2PreCE, + const double p_Mass1PreCE, + const double p_Mass2PreCE, + const bool p_EnvelopeFlag1, + const bool p_EnvelopeFlag2) { + + const bool supportedChannel = CircumbinaryDisk::IsSupportedCommonEnvelopeChannel( + p_StellarType1PreCE, + p_StellarType2PreCE, + p_Mass1PreCE, + p_Mass2PreCE, + p_EnvelopeFlag1, + p_EnvelopeFlag2 + ); + + if (!supportedChannel) { + return; + } + + const double diskLifetimeYears = OPTIONS->CircumbinaryDiskLifetime(); + const CircumbinaryDisk::FormationOptions formationOptions{CircumbinaryDisk::MakeFormationOptionsFromCompasOptions()}; + + const double mass1PostCEMsol = m_Star1->Mass(); + const double mass2PostCEMsol = m_Star2->Mass(); + const double semiMajorAxisPostCERsol = m_SemiMajorAxis * AU_TO_RSOL; + // The formation test uses the pre-RLOF orbit as the basis for the angular-momentum + // reservoir estimate. The post-CE separation is used only to define + // the inner CBD radius. + const CircumbinaryDisk::FormationInput circumbinaryDiskFormationInput{ + p_Mass1PreCE + p_Mass2PreCE, + mass1PostCEMsol + mass2PostCEMsol, + m_SemiMajorAxisPreRLOF * AU_TO_RSOL, + semiMajorAxisPostCERsol, + m_TotalAngularMomentumPreRLOF, + formationOptions + }; + + const CircumbinaryDisk::FormationResult formationResult{CircumbinaryDisk::EvaluateFormation(circumbinaryDiskFormationInput)}; + if (!formationResult.forms) { + return; + } + + // The formation step defines the intended CBD reservoir. The mass actually + // supplied can be smaller if another binary event truncates the disk. + const double initialReservoirMassMsol = formationResult.initialReservoirMassMsol; + + m_CircumbinaryDiskDetails.formed = true; + m_CircumbinaryDiskDetails.initialReservoirMass = initialReservoirMassMsol; + m_CircumbinaryDiskDetails.massSupplied = 0.0; + m_CircumbinaryDiskDetails.massRetained1 = 0.0; + m_CircumbinaryDiskDetails.massRetained2 = 0.0; + m_CircumbinaryDiskDetails.semiMajorAxisPre = semiMajorAxisPostCERsol; + m_CircumbinaryDiskDetails.semiMajorAxisPost = semiMajorAxisPostCERsol; + m_CircumbinaryDiskDetails.eccentricityPre = m_Eccentricity; + m_CircumbinaryDiskDetails.eccentricityPost = m_Eccentricity; + m_CircumbinaryDiskDetails.formationTime = m_Time; + m_CircumbinaryDiskDetails.mass1AtFormation = m_Star1->Mass(); + m_CircumbinaryDiskDetails.mass2AtFormation = m_Star2->Mass(); + m_CircumbinaryDiskDetails.stellarType1AtFormation = m_Star1->StellarType(); + m_CircumbinaryDiskDetails.stellarType2AtFormation = m_Star2->StellarType(); + + // Store the reservoir for evolution over subsequent COMPAS timesteps. + // This allows stellar radii, types, and Eddington limits to change while + // the disk is active. + m_ActiveCircumbinaryDisk.active = true; + m_ActiveCircumbinaryDisk.initialReservoirMass = initialReservoirMassMsol; + m_ActiveCircumbinaryDisk.remainingReservoirMass = initialReservoirMassMsol; + m_ActiveCircumbinaryDisk.supplyRateMsolPerYear = initialReservoirMassMsol / diskLifetimeYears; + m_ActiveCircumbinaryDisk.remainingLifetimeYears = diskLifetimeYears; + m_ActiveCircumbinaryDisk.massSupplied = 0.0; + m_ActiveCircumbinaryDisk.massRetained1 = 0.0; + m_ActiveCircumbinaryDisk.massRetained2 = 0.0; + m_RestartTimestepAfterImmediateEvent = true; + +} + + +/* + * Apply a CBD-driven stellar mass change through the normal COMPAS update path + * + * Non-compact accretors also have their effective initial mass, age, and + * rejuvenation state remapped after gaining mass. + * + * void ApplyCircumbinaryDiskMassChange(BinaryConstituentStar *p_Star, + * const double p_FinalMass) + * + * @param [IN] p_Star Star whose mass is to be updated + * @param [IN] p_FinalMass Final stellar mass after CBD accretion (Msol) + */ +void BaseBinaryStar::ApplyCircumbinaryDiskMassChange(BinaryConstituentStar *p_Star, const double p_FinalMass) { + + const double deltaMass = p_FinalMass - p_Star->Mass(); + if (utils::Compare(deltaMass, 0.0) == 0) return; + + const bool remapEffectiveInitialMass = utils::Compare(deltaMass, 0.0) > 0 && !p_Star->IsOneOf(COMPACT_OBJECTS); + + (void)p_Star->UpdateAttributes(deltaMass, 0.0, false); + + if (remapEffectiveInitialMass) { + p_Star->UpdateInitialMass(); + p_Star->UpdateAgeAfterMassLoss(); + p_Star->ApplyMassTransferRejuvenationFactor(); + p_Star->UpdateAttributes(0.0, 0.0, true); + } +} + + +/* + * FinaliseActiveCircumbinaryDisk() + * + * Deactivate the active finite-duration CBD and print the event record. + * + * The post-CBD orbital fields are deliberately not sampled here from the + * current binary state. They are cached immediately after each CBD increment + * in EvolveActiveCircumbinaryDiskOneTimestep(). This distinction matters when + * a later immediate event, such as a disruptive SN, terminates the disk: the CBD + * row should report the last state produced by CBD evolution, not the post-SN + * orbit. + * + * void FinaliseActiveCircumbinaryDisk() + */ +void BaseBinaryStar::FinaliseActiveCircumbinaryDisk() { + + const double initialReservoirMassMsol = std::max(m_ActiveCircumbinaryDisk.initialReservoirMass, 0.0); + const double completionToleranceMsol = std::max(CBD_RESERVOIR_COMPLETION_ABSOLUTE_TOLERANCE_MSOL, + CBD_RESERVOIR_COMPLETION_RELATIVE_TOLERANCE * std::max(initialReservoirMassMsol, 1.0)); + + const bool completedByMass = + utils::Compare(initialReservoirMassMsol, 0.0) > 0 && + utils::Compare(m_ActiveCircumbinaryDisk.remainingReservoirMass, 0.0, completionToleranceMsol) <= 0; + + const bool completedByDuration = + utils::Compare(initialReservoirMassMsol, 0.0) > 0 && + utils::Compare(m_ActiveCircumbinaryDisk.supplyRateMsolPerYear, 0.0) > 0 && + utils::Compare(m_ActiveCircumbinaryDisk.remainingLifetimeYears, 0.0) <= 0; + + const bool completed = completedByMass || completedByDuration; + + if (completed) { + // Under the constant-supply-rate prescription, exhausting the adopted + // lifetime and exhausting the intended reservoir are the same physical + // completion condition. Clamp away numerical residuals so the logged + // supplied mass exactly matches the initial reservoir. + m_ActiveCircumbinaryDisk.remainingReservoirMass = 0.0; + m_ActiveCircumbinaryDisk.remainingLifetimeYears = 0.0; + m_ActiveCircumbinaryDisk.massSupplied = initialReservoirMassMsol; + } + + // The live stellar masses may already have been + // changed by the event ending the CBD. Use the tracked CBD increments. + m_CircumbinaryDiskDetails.massSupplied = std::max(m_ActiveCircumbinaryDisk.massSupplied, 0.0); + m_CircumbinaryDiskDetails.massRetained1 = m_ActiveCircumbinaryDisk.massRetained1; + m_CircumbinaryDiskDetails.massRetained2 = m_ActiveCircumbinaryDisk.massRetained2; + m_ActiveCircumbinaryDisk.active = false; + + (void)PrintCircumbinaryDisk(); +} + + +/* + * EvolveActiveCircumbinaryDiskOneTimestep() + * + * Apply one finite-duration CBD increment during the current COMPAS timestep. + * This allows the stellar clock, wind mass loss, stellar radii/types, and + * Eddington limits to update between CBD increments instead of freezing those + * quantities for the full CBD lifetime at the CE event. The binary's orbital parameters + * can also change due to winds and tides. + * + * The disk mass supplied in this call is the smaller of supplyRate * dt and the + * remaining reservoir. If CBD accretion triggers an AIC, the ordinary + * post-increment supernova check in EvaluateBinary() handles it. + * + * void EvolveActiveCircumbinaryDiskOneTimestep(const double p_Dt) + * + * @param [IN] p_Dt Current COMPAS timestep (Myr) + */ +void BaseBinaryStar::EvolveActiveCircumbinaryDiskOneTimestep(const double p_Dt) { + + if (!m_ActiveCircumbinaryDisk.active) return; + + const double dtYears = std::min(p_Dt * MYR_TO_YEAR, m_ActiveCircumbinaryDisk.remainingLifetimeYears); + const double initialReservoirMassMsol = std::max(m_ActiveCircumbinaryDisk.initialReservoirMass, 0.0); + const double completionToleranceMsol = std::max(CBD_RESERVOIR_COMPLETION_ABSOLUTE_TOLERANCE_MSOL, + CBD_RESERVOIR_COMPLETION_RELATIVE_TOLERANCE * std::max(initialReservoirMassMsol, 1.0)); + + const bool reservoirExhausted = utils::Compare(m_ActiveCircumbinaryDisk.remainingReservoirMass, 0.0, completionToleranceMsol) <= 0; + if (utils::Compare(dtYears, 0.0) <= 0 || reservoirExhausted) { + FinaliseActiveCircumbinaryDisk(); + return; + } + + const double suppliedMassThisStepMsol = std::min( + m_ActiveCircumbinaryDisk.supplyRateMsolPerYear * dtYears, + m_ActiveCircumbinaryDisk.remainingReservoirMass + ); + + if (utils::Compare(suppliedMassThisStepMsol, 0.0) <= 0) { + FinaliseActiveCircumbinaryDisk(); + return; + } + + const double mass1BeforeCBDMsol = m_Star1->Mass(); + const double mass2BeforeCBDMsol = m_Star2->Mass(); + const CircumbinaryDisk::EvolutionOptions evolutionOptions{CircumbinaryDisk::MakeEvolutionOptionsFromCompasOptions()}; + + const CircumbinaryDisk::EvolutionInput evolutionInput{ + mass1BeforeCBDMsol, + mass2BeforeCBDMsol, + m_Star1->Radius(), + m_Star2->Radius(), + m_Star1->StellarType(), + m_Star2->StellarType(), + m_Eccentricity, + suppliedMassThisStepMsol, + dtYears, + evolutionOptions + }; + + const CircumbinaryDisk::EvolutionTrack track{CircumbinaryDisk::EvolveTrack(evolutionInput)}; + + const double finalMass1Msol = track.mass1Msol; + const double finalMass2Msol = track.mass2Msol; + + m_ActiveCircumbinaryDisk.massRetained1 += finalMass1Msol - mass1BeforeCBDMsol; + m_ActiveCircumbinaryDisk.massRetained2 += finalMass2Msol - mass2BeforeCBDMsol; + + m_SemiMajorAxis = m_SemiMajorAxis * track.semiMajorAxisFactor; + m_Eccentricity = track.eccentricity; + + ApplyCircumbinaryDiskMassChange(m_Star1, finalMass1Msol); + ApplyCircumbinaryDiskMassChange(m_Star2, finalMass2Msol); + + m_ActiveCircumbinaryDisk.remainingLifetimeYears = std::max(m_ActiveCircumbinaryDisk.remainingLifetimeYears - dtYears, 0.0); + m_ActiveCircumbinaryDisk.remainingReservoirMass = std::max(m_ActiveCircumbinaryDisk.remainingReservoirMass - suppliedMassThisStepMsol, 0.0); + m_ActiveCircumbinaryDisk.massSupplied += suppliedMassThisStepMsol; + + const bool completedByMass = utils::Compare(m_ActiveCircumbinaryDisk.remainingReservoirMass, 0.0, completionToleranceMsol) <= 0; + const bool completedByDuration = utils::Compare(m_ActiveCircumbinaryDisk.remainingLifetimeYears, 0.0) <= 0; + + // Mirror the cumulative disk increments before an interrupting event. + m_CircumbinaryDiskDetails.massSupplied = std::max(m_ActiveCircumbinaryDisk.massSupplied, 0.0); + m_CircumbinaryDiskDetails.massRetained1 = m_ActiveCircumbinaryDisk.massRetained1; + m_CircumbinaryDiskDetails.massRetained2 = m_ActiveCircumbinaryDisk.massRetained2; + m_CircumbinaryDiskDetails.semiMajorAxisPost = m_SemiMajorAxis * AU_TO_RSOL; + m_CircumbinaryDiskDetails.eccentricityPost = m_Eccentricity; + + if (completedByMass || completedByDuration) { + FinaliseActiveCircumbinaryDisk(); + } +} + + + +/* + * Calculate post-common-envelope eccentricity from user prescription + * + * double CalculatePostCommonEnvelopeEccentricity() + * + * @return Eccentricity to assign immediately after common envelope evolution + */ +double BaseBinaryStar::CalculatePostCommonEnvelopeEccentricity() const { + + // The PRE_RLOF_CAP and PRE_RLOF_FRACTION prescriptions use the eccentricity + // stored at first RLOF contact, before any RLOF circularisation. See InitialiseMassTransfer(). + switch (OPTIONS->PostCommonEnvelopeEccentricityPrescription()) { + case POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION::PRE_RLOF_CAP: { + return std::min(m_EccentricityPreRLOF, OPTIONS->PostCommonEnvelopeEccentricityCap()); + } + case POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION::PRE_RLOF_FRACTION: { + return OPTIONS->PostCommonEnvelopeEccentricityFraction() * m_EccentricityPreRLOF; + } + case POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION::CIRCULAR: { + return 0.0; + } + } + + return 0.0; +} + + /* * Resolve a Common Envelope Event * @@ -1543,6 +1885,11 @@ void BaseBinaryStar::EvaluateSupernovae() { */ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { + if (m_ActiveCircumbinaryDisk.active) { + // A new common-envelope event interrupts the finite-duration CBD. + FinaliseActiveCircumbinaryDisk(); + } + double alphaCE = OPTIONS->CommonEnvelopeAlpha(); // CE efficiency parameter double eccentricity = Eccentricity(); // current eccentricity (before CEE) @@ -1554,6 +1901,16 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { double omegaSpin2_pre_CE = m_Star2->Omega(); // star2 spin (before CEE) double semiMajorAxisAfterStage1 = 0.0; // semi-major axis after stage 1 (to remain zero unless using 2-stage CE formalism) + const STELLAR_TYPE stellarType1PreCE = m_Star1->StellarType(); // state before the CE/CBD update + const STELLAR_TYPE stellarType2PreCE = m_Star2->StellarType(); + const double mass1PreCE = m_Star1->Mass(); + const double mass2PreCE = m_Star2->Mass(); + + // Clear stale CBD event state at the start of a fresh CE. Any active + // disk has already been finalised at the start of this CE event. + if (!m_ActiveCircumbinaryDisk.active) { + m_CircumbinaryDiskDetails.formed = false; + } bool isDonorMS = false; // check for main sequence donor if (OPTIONS->AllowMainSequenceStarToSurviveCommonEnvelope()) { // allow main sequence stars to survive CEE? @@ -1694,11 +2051,12 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { THROW_ERROR(ERROR::UNKNOWN_CE_FORMALISM); // throw error } - double rRLdfin1 = m_SemiMajorAxis * CalculateRocheLobeRadius_Static(m_Mass1Final, m_Mass2Final); // Roche-lobe radius in AU after CEE, seen by star1 - double rRLdfin2 = m_SemiMajorAxis * CalculateRocheLobeRadius_Static(m_Mass2Final, m_Mass1Final); // Roche-lobe radius in AU after CEE, seen by star2 + m_Eccentricity = CalculatePostCommonEnvelopeEccentricity(); // assign post-CE eccentricity according to user prescription + double postCEPeriastronAU = m_SemiMajorAxis * (1.0 - m_Eccentricity); // Roche lobes should be checked at periastron + double rRLdfin1 = postCEPeriastronAU * CalculateRocheLobeRadius_Static(m_Mass1Final, m_Mass2Final); // Roche-lobe radius in AU after CEE, seen by star1 + double rRLdfin2 = postCEPeriastronAU * CalculateRocheLobeRadius_Static(m_Mass2Final, m_Mass1Final); // Roche-lobe radius in AU after CEE, seen by star2 double rRLdfin1Rsol = rRLdfin1 * AU_TO_RSOL; // Roche-lobe radius in Rsol after CEE, seen by star1 double rRLdfin2Rsol = rRLdfin2 * AU_TO_RSOL; // Roche-lobe radius in Rsol after CEE, seen by star2 - m_Eccentricity = 0.0; // we assume that a common envelope event (CEE) circularises the binary m_Star1->ResolveCommonEnvelopeAccretion(m_Mass1Final); // update star1's mass after CE accretion m_Star2->ResolveCommonEnvelopeAccretion(m_Mass2Final); // update star2's mass after CE accretion @@ -1754,12 +2112,24 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { m_Star2->SetPostCEEValues(); // squirrel away post CEE stellar values for star 2 SetPostCEEValues(m_SemiMajorAxis * AU_TO_RSOL, semiMajorAxisAfterStage1 * AU_TO_RSOL, m_Eccentricity, rRLdfin1Rsol, rRLdfin2Rsol); // squirrel away post CEE binary values (checks for post-CE RLOF, so should be done at end) - if (m_RLOFDetails.immediateRLOFPostCEE == true && !OPTIONS->AllowImmediateRLOFpostCEToSurviveCommonEnvelope()) {// is there immediate post-CE RLOF which is not allowed? + const bool allowImmediateRLOFPostCE = OPTIONS->AllowImmediateRLOFpostCEToSurviveCommonEnvelope(); + if ((!m_RLOFDetails.immediateRLOFPostCEE || allowImmediateRLOFPostCE) && OPTIONS->CircumbinaryDisk()) { // only evolve CBDs for CE survivors + // The CBD framework is a post-CE process. It is applied only after + // the ordinary CE calculation has produced a detached survivor; if + // the disk evolves during later timesteps, so the usual event + // checks handle any subsequent RLOF. + ApplyCircumbinaryDiskAfterCommonEnvelope(stellarType1PreCE, stellarType2PreCE, mass1PreCE, mass2PreCE, envelopeFlag1, envelopeFlag2); + } + + if (m_RLOFDetails.immediateRLOFPostCEE && !allowImmediateRLOFPostCE) { m_MassTransferTrackerHistory = MT_TRACKING::MERGER; m_Flags.stellarMerger = true; } + } + if (m_Flags.stellarMerger && m_ActiveCircumbinaryDisk.active) FinaliseActiveCircumbinaryDisk(); + (void)PrintCommonEnvelope(); // print (log) common envelope details } @@ -2210,6 +2580,14 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { // no m_MassTransferTrackerHistory = m_Donor == m_Star1 ? MT_TRACKING::STABLE_1_TO_2_SURV : MT_TRACKING::STABLE_2_TO_1_SURV; // record what happened - for later printing + if (m_ActiveCircumbinaryDisk.active && m_MassTransferTimescale == MT_TIMESCALE::THERMAL) { + // The CBD prescription is not coupled to the large discrete + // binary transformation represented by thermal-timescale stable + // MT. Nuclear-timescale stable MT may continue alongside the + // active CBD. + FinaliseActiveCircumbinaryDisk(); + } + double massGainAccretor = -massDiffDonor * m_FractionAccreted; // set accretor mass gain to mass loss * conservativeness double omegaDonor_pre_MT = m_Donor->Omega(); // used if full donor envelope is removed @@ -2370,10 +2748,27 @@ void BaseBinaryStar::InitialiseMassTransfer() { m_MassTransferTimescale = MT_TIMESCALE::NONE; m_MassLossRateInRLOF = 0.0; + // Store the binary state before this call can circularise the orbit, so + // post-CE eccentricity prescriptions can reference the pre-RLOF state. + const bool wasRLOF = m_RLOFDetails.isRLOF; + const double eccentricityAtRLOFCheck = m_Eccentricity; + const double semiMajorAxisAtRLOFCheck = m_SemiMajorAxis; + const double angularMomentumAtRLOFCheck = CalculateAngularMomentum(); + m_Star1->InitialiseMassTransfer(m_CEDetails.CEEnow, m_SemiMajorAxis, m_Eccentricity); // initialise mass transfer for star1 m_Star2->InitialiseMassTransfer(m_CEDetails.CEEnow, m_SemiMajorAxis, m_Eccentricity); // initialise mass transfer for star2 - if (m_Star1->IsRLOF() || m_Star2->IsRLOF()) { // either star overflowing its Roche Lobe? + const bool isRLOF = m_Star1->IsRLOF() || m_Star2->IsRLOF(); + + if (isRLOF && !wasRLOF) { + m_EccentricityPreRLOF = eccentricityAtRLOFCheck; // store eccentricity before any RLOF circularisation + m_SemiMajorAxisPreRLOF = semiMajorAxisAtRLOFCheck; // store semi-major axis before any RLOF circularisation + m_TotalAngularMomentumPreRLOF = angularMomentumAtRLOFCheck; + } + + m_RLOFDetails.isRLOF = isRLOF; + + if (isRLOF) { // either star overflowing its Roche Lobe? // yes - mass transfer if not both CH if (OPTIONS->CHEMode() != CHE_MODE::NONE && HasTwoOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS})) { // CHE enabled and both stars CH? // yes @@ -3262,8 +3657,18 @@ double BaseBinaryStar::ChooseTimestep(const double p_Factor) { } dt *= OPTIONS->TimestepMultiplier() * p_Factor; + dt = std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, TIDES_MINIMUM_FRACTIONAL_NUCLEAR_TIME * NUCLEAR_MINIMUM_TIMESTEP); // quantised and not less than minimum + + if (m_ActiveCircumbinaryDisk.active) { + dt = std::min(dt, m_ActiveCircumbinaryDisk.remainingLifetimeYears * YEAR_TO_MYR); + if (utils::Compare(m_ActiveCircumbinaryDisk.supplyRateMsolPerYear, 0.0) > 0) { + const double binaryMassMsol = m_Star1->Mass() + m_Star2->Mass(); + dt = std::min(dt, CBD_MAXIMUM_SUPPLIED_MASS_FRACTION_PER_TIMESTEP * binaryMassMsol / + m_ActiveCircumbinaryDisk.supplyRateMsolPerYear * YEAR_TO_MYR); + } + } - return std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, TIDES_MINIMUM_FRACTIONAL_NUCLEAR_TIME * NUCLEAR_MINIMUM_TIMESTEP); // quantised and not less than minimum + return dt; } @@ -3286,6 +3691,19 @@ double BaseBinaryStar::ChooseTimestep(const double p_Factor) { */ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { + // Finite-duration CBD evolution is applied before ordinary mass-transfer + // checks in each timestep. This lets a disk-driven orbital change affect + // whether ordinary COMPAS event checks, such as RLOF/SN handling, occur in + // that timestep. + if (m_ActiveCircumbinaryDisk.active && !m_RestartTimestepAfterImmediateEvent) { + EvolveActiveCircumbinaryDiskOneTimestep(p_Dt); + + if (m_RestartTimestepAfterImmediateEvent) return; + + m_SemiMajorAxisPrev = m_SemiMajorAxis; + m_EccentricityPrev = m_Eccentricity; + } + CalculateMassTransfer(p_Dt); // calculate mass transfer if necessary (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_MT); // print (log) detailed output @@ -3299,14 +3717,21 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { && !HasOneOf({STELLAR_TYPE::MASSLESS_REMNANT}) ) { // yes - avoid CEE if CH+CH or one star is a massless remnant ResolveCommonEnvelopeEvent(); // resolve CEE - immediate event - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_CEE); // print (log) detailed output + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_CEE); // print (log) detailed output } else if (m_Star1->IsSNevent() || m_Star2->IsSNevent()) { + const SN_EVENT supernovaEvents = // stash current SN types before EvaluateSupernovae() clears them + (m_Star1->IsSNevent() ? m_Star1->SN_Type() : SN_EVENT::NONE) | + (m_Star2->IsSNevent() ? m_Star2->SN_Type() : SN_EVENT::NONE); EvaluateSupernovae(); // evaluate supernovae (both stars) - immediate event (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_SN); // print (log) detailed output if (HasOneOf({ STELLAR_TYPE::NEUTRON_STAR })) { (void)PrintPulsarEvolutionParameters(BSE_PULSAR_RECORD_TYPE::POST_SN); // print (log) pulsar evolution parameters } + if (m_ActiveCircumbinaryDisk.active && // a live CBD survives non-disrupting bound events such as AIC + (CircumbinaryDisk::SupernovaEventsDisruptDisk(supernovaEvents) || IsUnbound())) { + FinaliseActiveCircumbinaryDisk(); // finalise and log the truncated CBD + } } else { ResolveMassChanges(); // apply mass loss and mass transfer as necessary @@ -3325,21 +3750,30 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { if (!StellarMerger() || (HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) && OPTIONS->EvolveMainSequenceMergers())) { // check stellar merger or evolving MS mergers // continue evolution if ((m_Star1->IsSNevent() || m_Star2->IsSNevent())) { + const SN_EVENT supernovaEvents = // stash current SN types before EvaluateSupernovae() clears them + (m_Star1->IsSNevent() ? m_Star1->SN_Type() : SN_EVENT::NONE) | + (m_Star2->IsSNevent() ? m_Star2->SN_Type() : SN_EVENT::NONE); EvaluateSupernovae(); // evaluate supernovae (both stars) if mass changes are responsible for a supernova (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_SN); // print (log) detailed output if (HasOneOf({ STELLAR_TYPE::NEUTRON_STAR })) { (void)PrintPulsarEvolutionParameters(BSE_PULSAR_RECORD_TYPE::POST_SN); // print (log) pulsar evolution parameters } + if (m_ActiveCircumbinaryDisk.active && // a live CBD survives non-disrupting bound events such as AIC + (CircumbinaryDisk::SupernovaEventsDisruptDisk(supernovaEvents) || IsUnbound())) { + FinaliseActiveCircumbinaryDisk(); // finalise and log the truncated CBD + } } - CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations - ProcessTides(p_Dt); // process tides if required - // assign new values to "previous" values, for following timestep - m_EccentricityPrev = m_Eccentricity; - m_SemiMajorAxisPrev = m_SemiMajorAxis; + if (!m_RestartTimestepAfterImmediateEvent) { // if an immediate event ended the timestep, do not process tides/spin with the old timestep + CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations + ProcessTides(p_Dt); // process tides if required + // assign new values to "previous" values, for following timestep + m_EccentricityPrev = m_Eccentricity; + m_SemiMajorAxisPrev = m_SemiMajorAxis; - m_Star1->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star1 - m_Star2->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star2 + m_Star1->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star1 + m_Star2->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star2 + } } } @@ -3518,15 +3952,21 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { // eventually, will want to implement a more careful prescription for the merger product, // perhaps allowing further evolution of the merger product (void)LogMergerToSwitchLog(); // log merger + if (m_ActiveCircumbinaryDisk.active) FinaliseActiveCircumbinaryDisk(); evolutionStatus = EVOLUTION_STATUS::STELLAR_MERGER; // stop evolution } } else if (HasStarsTouching()) { // binary components touching? (should usually be avoided as MT or CE or merger should happen prior to this) + if (m_ActiveCircumbinaryDisk.active) FinaliseActiveCircumbinaryDisk(); evolutionStatus = EVOLUTION_STATUS::STARS_TOUCHING; // yes - stop evolution } else if (IsUnbound()) { // binary is unbound? m_Flags.mergesInHubbleTime = false; // yes - won't merge in a Hubble time + if (m_ActiveCircumbinaryDisk.active) { // unbound binaries cannot retain an active CBD + FinaliseActiveCircumbinaryDisk(); // finalise and log the truncated CBD + } + if (IsDCO()) { // DCO (has two COs)? if (m_DCOFormationTime == DEFAULT_INITIAL_DOUBLE_VALUE) { // DCO not yet evaluated m_DCOFormationTime = m_Time; // set the DCO formation time @@ -3544,7 +3984,10 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { (void)PrintPulsarEvolutionParameters(BSE_PULSAR_RECORD_TYPE::POST_BINARY_TIMESTEP); // print (log) pulsar evolution parameters } - if (IsDCO() && !IsUnbound()) { // bound double compact object? + // DCO formation is not, by itself, a CBD-disrupting event. A bound DCO produced by a + // non-disrupting event (currently AIC) keeps the same disk, so defer DCO coalescence + // bookkeeping and the usual DCO stopping condition until the CBD has drained. + if (IsDCO() && !IsUnbound() && !m_ActiveCircumbinaryDisk.active) { // bound DCO with no active CBD? if (m_DCOFormationTime == DEFAULT_INITIAL_DOUBLE_VALUE) { // DCO not yet evaluated -- to ensure that the coalescence is only resolved once ResolveCoalescence(); // yes - resolve coalescence m_DCOFormationTime = m_Time; // set the DCO formation time @@ -3560,13 +4003,18 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? // yes // check for other reasons to stop evolution - if (IsDCO() && m_Time > (m_DCOFormationTime + m_TimeToCoalescence) && !IsUnbound()) { // evolution time exceeds DCO merger time? + if (IsDCO() && + !m_ActiveCircumbinaryDisk.active && + m_DCOFormationTime != DEFAULT_INITIAL_DOUBLE_VALUE && + m_Time > (m_DCOFormationTime + m_TimeToCoalescence) && + !IsUnbound()) { // evolution time exceeds DCO merger time? evolutionStatus = EVOLUTION_STATUS::DCO_MERGER_TIME; // yes - stop evolution } else if (m_Time > OPTIONS->MaxEvolutionTime()) { // evolution time exceeds maximum? evolutionStatus = EVOLUTION_STATUS::TIMES_UP; // yes - stop evolution } - else if (!OPTIONS->EvolveDoubleWhiteDwarfs() && IsWDandWD()) { // double WD and their evolution is not enabled? + else if (!m_ActiveCircumbinaryDisk.active && // allow a live finite-duration CBD to drain before stopping + !OPTIONS->EvolveDoubleWhiteDwarfs() && IsWDandWD()) { // double WD and their evolution is not enabled? evolutionStatus = EVOLUTION_STATUS::WD_WD; // yes - do not evolve double WD systems } else if ((HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) && !OPTIONS->EvolveMainSequenceMergers()) || @@ -3581,6 +4029,32 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { } } + if (evolutionStatus != EVOLUTION_STATUS::CONTINUE && m_ActiveCircumbinaryDisk.active) { // stop evolution with a truncated active CBD? + FinaliseActiveCircumbinaryDisk(); // yes - finalise and log the CBD before stopping + } + + if (evolutionStatus == EVOLUTION_STATUS::CONTINUE && m_RestartTimestepAfterImmediateEvent) { // immediate event changed the binary: choose a fresh timestep before aging stars + m_RestartTimestepAfterImmediateEvent = false; + + m_Star2->UpdatePreviousTimestepDuration(); + m_Star1->UpdatePreviousTimestepDuration(); + + if (usingProvidedTimesteps) { + if (stepNum >= timesteps.size()) { + evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED; + SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED); + } + else { + dt = timesteps[stepNum]; + stepNum++; + } + } + else { + if (OPTIONS->EmitGravitationalRadiation()) CalculateGravitationalRadiation(); + dt = ChooseTimestep(); + } + } + if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? // yes (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::PRE_STELLAR_TIMESTEP); // print (log) detailed output diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h index d6a416351..41c3c54dc 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -44,6 +44,9 @@ class BaseBinaryStar { m_CircularizationTimescale = p_Star.m_CircularizationTimescale; m_CEDetails = p_Star.m_CEDetails; + m_CircumbinaryDiskDetails = p_Star.m_CircumbinaryDiskDetails; + m_ActiveCircumbinaryDisk = p_Star.m_ActiveCircumbinaryDisk; + m_RestartTimestepAfterImmediateEvent = p_Star.m_RestartTimestepAfterImmediateEvent; m_Unbound = p_Star.m_Unbound; @@ -55,6 +58,7 @@ class BaseBinaryStar { m_EccentricityAtDCOFormation = p_Star.m_EccentricityAtDCOFormation; m_EccentricityInitial = p_Star.m_EccentricityInitial; m_EccentricityPreSN = p_Star.m_EccentricityPreSN; + m_EccentricityPreRLOF = p_Star.m_EccentricityPreRLOF; m_EccentricityPrev = p_Star.m_EccentricityPrev; m_Flags = p_Star.m_Flags; @@ -93,6 +97,7 @@ class BaseBinaryStar { m_SemiMajorAxisAtDCOFormation = p_Star.m_SemiMajorAxisAtDCOFormation; m_SemiMajorAxisInitial = p_Star.m_SemiMajorAxisInitial; m_SemiMajorAxisPreSN = p_Star.m_SemiMajorAxisPreSN; + m_SemiMajorAxisPreRLOF = p_Star.m_SemiMajorAxisPreRLOF; m_SemiMajorAxisPrev = p_Star.m_SemiMajorAxisPrev; m_SupernovaState = p_Star.m_SupernovaState; @@ -112,6 +117,7 @@ class BaseBinaryStar { m_TotalAngularMomentumPrev = p_Star.m_TotalAngularMomentumPrev; m_TotalAngularMomentum = p_Star.m_TotalAngularMomentum; + m_TotalAngularMomentumPreRLOF = p_Star.m_TotalAngularMomentumPreRLOF; m_TotalEnergy = p_Star.m_TotalEnergy; @@ -162,6 +168,20 @@ class BaseBinaryStar { bool CEAtLeastOnce() const { return m_CEDetails.CEEcount > 0; } unsigned int CEEventCount() const { return m_CEDetails.CEEcount; } double CircularizationTimescale() const { return m_CircularizationTimescale; } + bool CircumbinaryDiskFormed() const { return m_CircumbinaryDiskDetails.formed; } + double CircumbinaryDiskInitialMass() const { return m_CircumbinaryDiskDetails.initialReservoirMass; } + double CircumbinaryDiskMassSupplied() const { return m_CircumbinaryDiskDetails.massSupplied; } + double CircumbinaryDiskMassRetained1() const { return m_CircumbinaryDiskDetails.massRetained1; } + double CircumbinaryDiskMassRetained2() const { return m_CircumbinaryDiskDetails.massRetained2; } + double CircumbinaryDiskSemiMajorAxisPre() const { return m_CircumbinaryDiskDetails.semiMajorAxisPre; } + double CircumbinaryDiskSemiMajorAxisPost() const { return m_CircumbinaryDiskDetails.semiMajorAxisPost; } + double CircumbinaryDiskEccentricityPre() const { return m_CircumbinaryDiskDetails.eccentricityPre; } + double CircumbinaryDiskEccentricityPost() const { return m_CircumbinaryDiskDetails.eccentricityPost; } + double CircumbinaryDiskFormationTime() const { return m_CircumbinaryDiskDetails.formationTime; } + double CircumbinaryDiskMass1AtFormation() const { return m_CircumbinaryDiskDetails.mass1AtFormation; } + double CircumbinaryDiskMass2AtFormation() const { return m_CircumbinaryDiskDetails.mass2AtFormation; } + STELLAR_TYPE CircumbinaryDiskStellarType1AtFormation() const { return m_CircumbinaryDiskDetails.stellarType1AtFormation; } + STELLAR_TYPE CircumbinaryDiskStellarType2AtFormation() const { return m_CircumbinaryDiskDetails.stellarType2AtFormation; } unsigned int CommonEnvelopeEventCount() const { return m_CEDetails.CEEcount; } bool Unbound() const { return m_Unbound; } bool DoubleCoreCE() const { return m_CEDetails.doubleCoreCE; } @@ -309,6 +329,39 @@ class BaseBinaryStar { BinaryCEDetailsT m_CEDetails; // Common Event details + // Cached initial and final quantities written to one CBD event record. + struct CircumbinaryDiskDetailsT { + bool formed = false; // A post-CE circumbinary disk formed and was evolved + double initialReservoirMass = DEFAULT_INITIAL_DOUBLE_VALUE; // Initial CBD reservoir intended to be supplied (Msol) + double massSupplied = DEFAULT_INITIAL_DOUBLE_VALUE; // Mass actually supplied to the binary through the CBD (Msol) + double massRetained1 = DEFAULT_INITIAL_DOUBLE_VALUE; // Mass retained by star 1 during CBD accretion (Msol) + double massRetained2 = DEFAULT_INITIAL_DOUBLE_VALUE; // Mass retained by star 2 during CBD accretion (Msol) + double semiMajorAxisPre = DEFAULT_INITIAL_DOUBLE_VALUE; // Semi-major axis before CBD evolution (Rsol) + double semiMajorAxisPost = DEFAULT_INITIAL_DOUBLE_VALUE; // Semi-major axis after CBD evolution (Rsol) + double eccentricityPre = DEFAULT_INITIAL_DOUBLE_VALUE; // Eccentricity before CBD evolution + double eccentricityPost = DEFAULT_INITIAL_DOUBLE_VALUE; // Eccentricity after CBD evolution + double formationTime = DEFAULT_INITIAL_DOUBLE_VALUE; // COMPAS time at CBD formation (Myr) + double mass1AtFormation = DEFAULT_INITIAL_DOUBLE_VALUE; // Star 1 mass at CBD formation (Msol) + double mass2AtFormation = DEFAULT_INITIAL_DOUBLE_VALUE; // Star 2 mass at CBD formation (Msol) + STELLAR_TYPE stellarType1AtFormation = STELLAR_TYPE::NONE; // Star 1 type at CBD formation + STELLAR_TYPE stellarType2AtFormation = STELLAR_TYPE::NONE; // Star 2 type at CBD formation + } m_CircumbinaryDiskDetails{}; + + // Persistent reservoir and delivery state for an active CBD that is + // still being evolved alongside ordinary COMPAS timesteps. + struct ActiveCircumbinaryDiskT { + bool active = false; // CBD is currently draining over ordinary COMPAS timesteps + double initialReservoirMass = 0.0; // Intended mass of CBD reservoir at formation (Msol) + double remainingReservoirMass = 0.0; // CBD reservoir mass not yet supplied (Msol) + double supplyRateMsolPerYear = 0.0; // Uniform CBD supply rate (Msol yr^-1) + double remainingLifetimeYears = 0.0; // CBD lifetime still to be evolved (yr) + double massSupplied = 0.0; // Cumulative mass actually supplied through the CBD (Msol) + double massRetained1 = 0.0; // Cumulative CBD-supplied mass retained by star 1 (Msol) + double massRetained2 = 0.0; // Cumulative CBD-supplied mass retained by star 2 (Msol) + } m_ActiveCircumbinaryDisk{}; + + bool m_RestartTimestepAfterImmediateEvent; // Restart the timestep before aging stars after an instantaneous event + double m_CircularizationTimescale; bool m_Unbound; // Binary unbound? @@ -319,6 +372,7 @@ class BaseBinaryStar { double m_EccentricityAtDCOFormation; // Eccentricity at DCO formation double m_EccentricityInitial; // Record initial eccentricity double m_EccentricityPreSN; // Eccentricity prior to supernova + double m_EccentricityPreRLOF; // Eccentricity immediately before first RLOF in the current episode double m_EccentricityPrev; // Eccentricity at previous timestep struct FLAGS { // Miscellaneous flags @@ -365,6 +419,7 @@ class BaseBinaryStar { double m_SemiMajorAxisAtDCOFormation; // Semi-major axis at DCO formation double m_SemiMajorAxisInitial; // Record initial semi-major axis double m_SemiMajorAxisPreSN; // Semi-major axis prior to supernova + double m_SemiMajorAxisPreRLOF; // Semi-major axis immediately before first RLOF in the current episode double m_SemiMajorAxisPrev; // Semi-major axis at previous timestep SN_STATE m_SupernovaState; // Indicates which star (or stars) are undergoing / have undergone a supernova event @@ -384,6 +439,7 @@ class BaseBinaryStar { double m_DCOFormationTime; // Time of DCO formation double m_TotalAngularMomentum; + double m_TotalAngularMomentumPreRLOF; double m_TotalAngularMomentumPrev; double m_TotalEnergy; @@ -493,6 +549,22 @@ class BaseBinaryStar { void EvaluateSupernovae(); + // CE/CBD helper functions. The CBD physics lives in the CircumbinaryDisk* implementation files; + // these helpers keep ResolveCommonEnvelopeEvent() focused on the CE flow + // while handling COMPAS-specific state mutation, immediate-RLOF checks, + // AIC state changes, and CBD event logging. + double CalculatePostCommonEnvelopeEccentricity() const; + void ApplyCircumbinaryDiskAfterCommonEnvelope(const STELLAR_TYPE p_StellarType1PreCE, + const STELLAR_TYPE p_StellarType2PreCE, + const double p_Mass1PreCE, + const double p_Mass2PreCE, + const bool p_EnvelopeFlag1, + const bool p_EnvelopeFlag2); + void EvolveActiveCircumbinaryDiskOneTimestep(const double p_Dt); + void FinaliseActiveCircumbinaryDisk(); + void ApplyCircumbinaryDiskMassChange(BinaryConstituentStar *p_Star, + const double p_FinalMass); + ERROR EvolveOneTimestep(const double p_Dt); void EvolveOneTimestepPreamble(const double p_Dt); @@ -563,6 +635,11 @@ class BaseBinaryStar { bool PrintCommonEnvelope(const CE_RECORD_TYPE p_RecordType = CE_RECORD_TYPE::DEFAULT) const { return LOGGING->LogCommonEnvelope(this, p_RecordType); } + + // Write the current CBD event through the standard BSE CBD logfile. + bool PrintCircumbinaryDisk() const { + return LOGGING->LogCircumbinaryDisk(this); + } bool PrintPulsarEvolutionParameters(const BSE_PULSAR_RECORD_TYPE p_RecordType) const { return OPTIONS->EvolvePulsars() ? LOGGING->LogBSEPulsarEvolutionParameters(this, p_RecordType) : true; diff --git a/src/CircumbinaryDisk.h b/src/CircumbinaryDisk.h new file mode 100644 index 000000000..d77fefcda --- /dev/null +++ b/src/CircumbinaryDisk.h @@ -0,0 +1,124 @@ +#ifndef __CircumbinaryDisk_h__ +#define __CircumbinaryDisk_h__ + +#include "constants.h" +#include "typedefs.h" + +namespace CircumbinaryDisk { + +// Public API for the post-CE circumbinary-disk model. +// +// BaseBinaryStar owns COMPAS state, event ordering, logging, and mutation of +// stars/binaries. This module receives explicit physical inputs, evaluates the +// CBD formation/reservoir model, and returns CBD-driven binary changes without +// mutating COMPAS state. +// +// Validity domain: the Spindler-based evolution treats the CBD as a +// perturbative reservoir whose mass is dynamically negligible compared with +// the central binary. The formalism therefore assumes m_CBD << M_binary. +// Here m_CBD is the initial intended reservoir (FormationResult::initialReservoirMassMsol +// and the logged Initial_CBD_Mass), M_binary is the post-CE central-binary mass. The code +// does not impose a hard numerical cutoff on this asymptotic assumption, so +// parameter choices outside this regime are extrapolations. +// +// Public structs use COMPAS star labels: star 1 and star 2 are never reordered +// at the module boundary. Internally, the implementation handles the ordered +// q <= 1 convention required by the Siwek/Spindler tables and converts results +// back to COMPAS labels before returning them. +// +// Unit convention at the public boundary: +// masses : Msol +// radii : Rsol +// separations : Rsol, except semiMajorAxisFactor which is dimensionless +// durations : yr +// eccentricity : dimensionless + +// User/model parameters for the formation criterion. These are populated +// from COMPAS options by MakeFormationOptionsFromCompasOptions(); this struct +// intentionally carries no physical defaults. +struct FormationOptions { + double circumbinaryDiskBeta; + double envelopeMassFractionSupplied; // Assumed to yield initial m_CBD << post-CE M_binary + double diskStructureFactor; + double diskInnerRadiusOverSeparation; +}; + +// User/model parameters for the CBD evolution step. These are populated +// from COMPAS options by MakeEvolutionOptionsFromCompasOptions(). +struct EvolutionOptions { + CIRCUMBINARY_DISK_EDDINGTON_MODE eddingtonMode; + double eddingtonCapFactor; + double eddingtonGe23Factor; +}; + +// Complete physical state needed for the CBD formation criterion. +// The caller identifies the compact-object companion and the envelope donor +// from the pre-CE COMPAS star labels and supplies their physical properties. +// beta is fJ/fM, where fJ is the fraction of the pre-CE binary angular momentum +// assigned to the CBD reservoir and fM is the supplied fraction of the total +// binary mass lost during CE. Thus J_CBD = beta * fM * J_available. The +// formation criterion is evaluated as a total angular-momentum comparison in +// COMPAS physical units. The available mass is calculated from the decrease +// in total binary mass across CE. The resulting supplied mass is assumed to +// be much smaller than the post-CE binary mass in the subsequent Spindler evolution. +struct FormationInput { + double binaryMassPreCEMsol; + double binaryMassPostCEMsol; + double semiMajorAxisPreRLOFRsol; + double semiMajorAxisPostCERsol; + double totalAngularMomentumPreRLOF; + FormationOptions options{}; +}; + +// Outcome of the local CBD formation calculation, including the initial +// intended reservoir when the angular-momentum criterion is satisfied. +struct FormationResult { + bool forms = false; + double initialReservoirMassMsol = 0.0; // Intended initial CBD reservoir; assumed to be << post-CE binary mass +}; + +// Complete physical state needed for one Spindler/CBD evolution step after a +// CBD has formed. The supplied mass is the "bare" mass passed through the +// disk model; Eddington limiting determines how much of it is retained by each +// star. +struct EvolutionInput { + double mass1Msol; + double mass2Msol; + double radius1Rsol; + double radius2Rsol; + STELLAR_TYPE stellarType1; + STELLAR_TYPE stellarType2; + double eccentricity; + double suppliedMassMsol; + double durationYears; + EvolutionOptions options{}; +}; + +// COMPAS-labelled binary state returned after evolving one complete CBD +// increment. semiMajorAxisFactor multiplies the input separation. +struct EvolutionTrack { + double semiMajorAxisFactor = 1.0; + double eccentricity = 0.0; + double mass1Msol = 0.0; + double mass2Msol = 0.0; +}; + +bool SupernovaEventsDisruptDisk(const SN_EVENT p_Events); + +bool IsSupportedCommonEnvelopeChannel(const STELLAR_TYPE p_StellarType1PreCE, + const STELLAR_TYPE p_StellarType2PreCE, + const double p_Mass1PreCE, + const double p_Mass2PreCE, + const bool p_EnvelopeFlag1, + const bool p_EnvelopeFlag2); + +FormationOptions MakeFormationOptionsFromCompasOptions(); +EvolutionOptions MakeEvolutionOptionsFromCompasOptions(); + +FormationResult EvaluateFormation(const FormationInput& p_Input); + +EvolutionTrack EvolveTrack(const EvolutionInput& p_Input); + +} // namespace CircumbinaryDisk + +#endif // __CircumbinaryDisk_h__ diff --git a/src/CircumbinaryDiskEvolutionSpindler.cpp b/src/CircumbinaryDiskEvolutionSpindler.cpp new file mode 100644 index 000000000..69a5aeb9f --- /dev/null +++ b/src/CircumbinaryDiskEvolutionSpindler.cpp @@ -0,0 +1,666 @@ +#include "CircumbinaryDisk.h" + +#include "constants.h" +#include "utils.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +// spindler-c is a C library. Its public header currently does not wrap the +// declarations in extern "C", so we do it here for C++/COMPAS builds. +extern "C" { +#include "spindler.h" +} + +namespace odeint = boost::numeric::odeint; + +// Physics scope of this module: +// * Disk formation is handled in CircumbinaryDiskFormation.cpp. This file +// assumes a CBD reservoir has already formed and evolves the binary using +// the Siwek/Spindler tabulated derivatives for circumbinary-disk interaction. +// * The supplied reservoir is delivered at a constant rate over the chosen +// disk lifetime. Each component accretes the table-implied fraction unless +// limited by the selected Eddington prescription; unretained mass removes +// angular momentum from the binary using the accretor's specific angular momentum. + +namespace { + +// Backend details and ODE solver settings are local to the Spindler implementation. +constexpr char CBD_SPINDLER_MODEL_NAME[] = "Siwek23"; +constexpr double CBD_SPINDLER_ODE_RELATIVE_TOLERANCE = 1.0E-4; +constexpr double CBD_SPINDLER_ODE_ABSOLUTE_TOLERANCE = 1.0E-5; + +// Physical and dimensionless inputs used by the internal Spindler integration. +// suppliedMassFraction is the total mass supplied by the CBD in units of the +// initial post-CE binary mass and is assumed to be much smaller than unity; a +// small integration increment does not by itself make an initially large +// reservoir consistent with that assumption. initialMassRatio follows the ordered +// Spindler convention q = M2/M1 <= 1, not necessarily the original COMPAS +// star-2/star-1 label ratio. initialBinaryMassMsol converts between physical Msol +// units and the dimensionless ODE variables. +struct SpindlerEvolutionParams { + double suppliedMassFraction; + double initialMassRatio; + double initialEccentricity; + double initialBinaryMassMsol; + double durationYears; + + STELLAR_TYPE stellarType1 = STELLAR_TYPE::NONE; + STELLAR_TYPE stellarType2 = STELLAR_TYPE::NONE; + double radius1Rsol = 0.0; + double radius2Rsol = 0.0; + + CircumbinaryDisk::EvolutionOptions options{}; +}; + +// Final dimensionless state produced by the internal Spindler solver before +// component ordering is converted back to COMPAS labels. +struct SpindlerEvolutionResult { + double semiMajorAxisFactor = 1.0; + double mass1Fraction = 0.5; + double mass2Fraction = 0.5; + double eccentricity = 0.0; +}; + +// Owns one initialized spindler-c table set and evaluates CBD-driven binary +// evolution with COMPAS Eddington-limit and angular-momentum-loss extensions. +class SpindlerBackend { +public: + explicit SpindlerBackend(const std::string& p_ModelName = CBD_SPINDLER_MODEL_NAME); + ~SpindlerBackend(); + + // The backend uniquely owns spindler-c table data, so copying and moving are disabled. + SpindlerBackend(const SpindlerBackend&) = delete; + SpindlerBackend& operator=(const SpindlerBackend&) = delete; + SpindlerBackend(SpindlerBackend&&) = delete; + SpindlerBackend& operator=(SpindlerBackend&&) = delete; + // Evolve with Eddington-limited, non-conservative accretion. + SpindlerEvolutionResult Evolve(const SpindlerEvolutionParams& p_Params); + +private: + // Spindler logarithmic orbital and mass-ratio derivatives returned for one + // labelled binary state. + struct Derivatives { + double da; + double dq; + double de; + }; + + Derivatives DerivativesForLabeledQ(double p_QLabeled, double p_Eccentricity) const; + std::pair ComponentMassFractions(double p_QLabeled, double p_Eccentricity) const; + double ComputeEddingtonRate( + const STELLAR_TYPE p_StellarType, + double p_MassFraction, + double p_InitialBinaryMassMsol, + double p_RadiusRsol + ) const; + + static double ApplyEddingtonLimit( + double p_RequestedRate, + double p_EddingtonRate, + CIRCUMBINARY_DISK_EDDINGTON_MODE p_Mode, + double p_CapFactor, + double p_Ge23Factor + ); + + // Dimensionless Spindler-internal orbital angular momentum: G = 1, a0 = 1, and Mtot0 = 1. + static double SpindlerOrbitalAngularMomentum(double p_Mass1Fraction, double p_Mass2Fraction, double p_Eccentricity, double p_SemiMajorAxisFactor); + + struct spindler_data_t m_SpindlerData{}; +}; + +#ifndef SPINDLER_C_ROOT +#define SPINDLER_C_ROOT "" +#endif + +/* + * Validate the ordered initial mass ratio and eccentricity before constructing the ODE state + * + * void ValidateCommonInputs(const double p_InitialMassRatio, const double p_InitialEccentricity) + * + * @param [IN] p_InitialMassRatio Initial ordered mass ratio, M2/M1 <= 1 + * @param [IN] p_InitialEccentricity Initial eccentricity + */ +void ValidateCommonInputs(const double p_InitialMassRatio, const double p_InitialEccentricity) { + if (utils::Compare(p_InitialMassRatio, 0.0) <= 0 || !std::isfinite(p_InitialMassRatio)) { + throw std::invalid_argument("q0 must be positive and finite"); + } + if (!std::isfinite(p_InitialEccentricity) || utils::Compare(p_InitialEccentricity, 0.0) < 0 || utils::Compare(p_InitialEccentricity, 1.0) >= 0) { + throw std::invalid_argument("e0 must satisfy 0 <= e0 < 1"); + } +} + +/* + * Convert the dimensionless ODE state into the internal Spindler result + * + * SpindlerEvolutionResult MakeResultFromState(const std::vector& p_State) + * + * @param [IN] p_State ODE state [a/a0, M1/M0, M2/M0, e] + * @return Internal evolution result for this state + */ +SpindlerEvolutionResult MakeResultFromState(const std::vector& p_State) { + SpindlerEvolutionResult result{}; + result.semiMajorAxisFactor = p_State[0]; + result.mass1Fraction = p_State[1]; + result.mass2Fraction = p_State[2]; + result.eccentricity = p_State[3]; + return result; +} + +std::mutex spindlerInitMutex; + +// spindler_init() resolves its table paths relative to the process working +// directory. This RAII guard temporarily switches to the spindler-c root and +// guarantees that COMPAS's original working directory is restored on exit. +class ScopedWorkingDirectory { +public: + /* + * Temporarily change the process working directory + * + * explicit ScopedWorkingDirectory(const std::string& p_NewDirectory) + * + * @param [IN] p_NewDirectory Directory required while spindler_init() resolves table paths + */ + explicit ScopedWorkingDirectory(const std::string& p_NewDirectory) { + if (p_NewDirectory.empty()) return; + + char* cwd = getcwd(nullptr, 0); + if (cwd == nullptr) { + std::ostringstream msg; + msg << "getcwd failed before spindler_init: " << std::strerror(errno); + throw std::runtime_error(msg.str()); + } + m_OldDirectory = cwd; + std::free(cwd); + + if (chdir(p_NewDirectory.c_str()) != 0) { + std::ostringstream msg; + msg << "Could not chdir to SPINDLER_C_ROOT='" << p_NewDirectory + << "' before spindler_init: " << std::strerror(errno); + throw std::runtime_error(msg.str()); + } + m_Changed = true; + } + + /* + * Restore the process working directory saved by the constructor + * + * ~ScopedWorkingDirectory() + */ + ~ScopedWorkingDirectory() { + if (m_Changed) { + (void)chdir(m_OldDirectory.c_str()); + } + } + + // Keep directory-restoration ownership unique to a single scope guard. + ScopedWorkingDirectory(const ScopedWorkingDirectory&) = delete; + ScopedWorkingDirectory& operator=(const ScopedWorkingDirectory&) = delete; + +private: + std::string m_OldDirectory; + bool m_Changed = false; +}; + +/* + * Initialise one spindler-c table backend + * + * SpindlerBackend(const std::string& p_ModelName) + * + * @param [IN] p_ModelName Name of the spindler-c table model to load + */ +SpindlerBackend::SpindlerBackend(const std::string& p_ModelName) { + std::vector modelNameCString(p_ModelName.begin(), p_ModelName.end()); + modelNameCString.push_back('\0'); + + // If enabled, the Makefile defines SPINDLER_C_ROOT to the user-supplied + // spindler-c root. spindler-c currently looks for tables//... + // relative to the current working directory during spindler_init(), so we + // temporarily initialize from that directory and then restore COMPAS's + // original working directory. + + std::lock_guard initLock(spindlerInitMutex); + ScopedWorkingDirectory cwdGuard{SPINDLER_C_ROOT}; + + const int err = spindler_init(modelNameCString.data(), &m_SpindlerData); + if (err != SPINDLER_NO_ERROR) { + std::ostringstream msg; + msg << "spindler_init failed for model '" << p_ModelName + << "' with error code " << err + << ". Check SPINDLER_C_ROOT='" << SPINDLER_C_ROOT + << "' and verify that tables/" << p_ModelName << " exists there."; + throw std::runtime_error(msg.str()); + } +} + +/* + * Release the table data allocated by spindler_init() + * + * ~SpindlerBackend() + */ +SpindlerBackend::~SpindlerBackend() { + spindler_free_data(&m_SpindlerData); +} + +/* + * Query the tabulated logarithmic orbital derivatives + * + * Derivatives DerivativesForLabeledQ(double p_QLabeled, double p_Eccentricity) const + * + * @param [IN] p_QLabeled Mass ratio in the caller's current component labels + * @param [IN] p_Eccentricity Current eccentricity + * @return Da, Dq, and De mapped back to the caller's labels + */ +SpindlerBackend::Derivatives SpindlerBackend::DerivativesForLabeledQ(double p_QLabeled, double p_Eccentricity) const { + if (utils::Compare(p_QLabeled, 0.0) <= 0 || !std::isfinite(p_QLabeled)) { + throw std::runtime_error("Encountered non-positive or non-finite mass ratio during Spindler integration"); + } + + const bool swapped = utils::Compare(p_QLabeled, 1.0) > 0; + const double qOrderedRaw = swapped ? 1.0 / p_QLabeled : p_QLabeled; + const double eccentricityRaw = std::max(p_Eccentricity, 0.0); + + // The table coefficients are logarithmic responses of the binary to CBD + // accretion/torques. Da multiplies a, De multiplies e, and Dq describes the + // evolution of the ordered mass ratio. If COMPAS's labelled mass ratio is + // inverted relative to q <= 1, the sign of the q derivative must be inverted + // before returning to the caller's labels. + const double qLookup = std::clamp(qOrderedRaw, CBD_SIWEK23_Q_MIN, CBD_SIWEK23_Q_MAX); + const double eccentricityLookup = std::clamp(eccentricityRaw, CBD_SIWEK23_E_MIN, CBD_SIWEK23_E_MAX); + + Derivatives derivatives{}; + derivatives.da = spindler_get_Da(qLookup, eccentricityLookup, const_cast(&m_SpindlerData)); + + const double deLookup = spindler_get_De(qLookup, eccentricityLookup, const_cast(&m_SpindlerData)); + derivatives.de = utils::Compare(eccentricityRaw, 0.0) > 0 ? deLookup * eccentricityLookup / eccentricityRaw : 0.0; + + const double dqLookup = spindler_get_Dq(qLookup, eccentricityLookup, const_cast(&m_SpindlerData)); + const double dqOrdered = dqLookup * qLookup / qOrderedRaw; + derivatives.dq = swapped ? -dqOrdered : dqOrdered; + return derivatives; +} + +/* + * Infer the component shares of the supplied CBD mass + * + * std::pair ComponentMassFractions(double p_QLabeled, double p_Eccentricity) const + * + * @param [IN] p_QLabeled Mass ratio in the caller's current component labels + * @param [IN] p_Eccentricity Current eccentricity + * @return Fractions of supplied mass assigned to components 1 and 2 + */ +std::pair SpindlerBackend::ComponentMassFractions(double p_QLabeled, double p_Eccentricity) const { + if (utils::Compare(p_QLabeled, 0.0) <= 0 || !std::isfinite(p_QLabeled)) { + throw std::runtime_error("Encountered non-positive or non-finite mass ratio while computing accretion split"); + } + + const bool swapped = utils::Compare(p_QLabeled, 1.0) > 0; + const double qOrderedRaw = swapped ? 1.0 / p_QLabeled : p_QLabeled; + const double eccentricityRaw = std::max(p_Eccentricity, 0.0); + const double qLookup = std::clamp(qOrderedRaw, CBD_SIWEK23_Q_MIN, CBD_SIWEK23_Q_MAX); + const double eccentricityLookup = std::clamp(eccentricityRaw, CBD_SIWEK23_E_MIN, CBD_SIWEK23_E_MAX); + + const double dqLookup = spindler_get_Dq(qLookup, eccentricityLookup, const_cast(&m_SpindlerData)); + const double dqOrdered = dqLookup * qLookup / qOrderedRaw; + + // Infer how the supplied CBD mass is split between the two components from + // the tabulated q evolution. This keeps the mass growth consistent with the + // same table derivative that drives the orbital response instead of imposing + // an independent accretion-ratio prescription. + double secondaryAccretionFraction = qOrderedRaw * (dqOrdered + 1.0 + qOrderedRaw) / std::pow(1.0 + qOrderedRaw, 2.0); + secondaryAccretionFraction = std::clamp(secondaryAccretionFraction, 0.0, 1.0); + const double primaryAccretionFraction = 1.0 - secondaryAccretionFraction; + + if (swapped) { + return {secondaryAccretionFraction, primaryAccretionFraction}; + } + return {primaryAccretionFraction, secondaryAccretionFraction}; +} + +/* + * Calculate the Eddington accretion rate for one component + * + * double ComputeEddingtonRate(const STELLAR_TYPE p_StellarType, + * double p_MassFraction, + * double p_InitialBinaryMassMsol, + * double p_RadiusRsol) const + * + * @param [IN] p_StellarType Current stellar type of the component + * @param [IN] p_MassFraction Component mass divided by the initial binary mass + * @param [IN] p_InitialBinaryMassMsol Initial binary mass used for ODE normalisation (Msol) + * @param [IN] p_RadiusRsol Current component radius (Rsol) + * @return Eddington rate in dimensionless Spindler mass units per year + */ +double SpindlerBackend::ComputeEddingtonRate( + const STELLAR_TYPE p_StellarType, + double p_MassFraction, + double p_InitialBinaryMassMsol, + double p_RadiusRsol +) const { + if (p_StellarType == STELLAR_TYPE::BLACK_HOLE) { + // Black-hole Eddington mass accretion rate = 4.42e-7 * (M_phys / 10) Msol/yr. Because the ODE masses + // are normalised by the initial binary mass, this is returned in + // dimensionless Spindler mass units per year. + return CBD_BH_EDDINGTON_RATE_AT_REFERENCE_MASS * (p_MassFraction / CBD_BH_EDDINGTON_REFERENCE_MASS_MSOL); + } + + const bool isWhiteDwarf = p_StellarType == STELLAR_TYPE::HELIUM_WHITE_DWARF || + p_StellarType == STELLAR_TYPE::CARBON_OXYGEN_WHITE_DWARF || + p_StellarType == STELLAR_TYPE::OXYGEN_NEON_WHITE_DWARF; + + if (p_StellarType == STELLAR_TYPE::NEUTRON_STAR || isWhiteDwarf) { + if (utils::Compare(p_RadiusRsol, 0.0) <= 0 || !std::isfinite(p_RadiusRsol)) { + throw std::invalid_argument("radius must be positive and finite for radius-scaled Eddington limits"); + } + if (utils::Compare(p_InitialBinaryMassMsol, 0.0) <= 0 || !std::isfinite(p_InitialBinaryMassMsol)) { + throw std::invalid_argument("initialBinaryMassMsol must be positive for radius-scaled Eddington limits"); + } + + // COMPAS uses a radius-scaled Eddington rate for white dwarfs and other + // non-BH accretors. Preserve the existing CBD neutron-star normalisation + // while applying the standard COMPAS/Hurley coefficient to white dwarfs. + const double eddingtonRateMsolPerYear = p_StellarType == STELLAR_TYPE::NEUTRON_STAR + ? CBD_NS_EDDINGTON_RATE_AT_REFERENCE_RADIUS * (p_RadiusRsol / CBD_NS_EDDINGTON_REFERENCE_RADIUS_RSOL) + : CBD_WD_EDDINGTON_RATE_PER_RADIUS * p_RadiusRsol; + return eddingtonRateMsolPerYear / p_InitialBinaryMassMsol; + } + + // The Eddington limit is assumed to be irrelevant for non-compact accretors. + return std::numeric_limits::infinity(); +} + +/* + * Apply the selected Eddington-limited retention prescription + * + * CAP imposes a hard multiple of the Eddington rate, TM23 reduces the + * super-Eddington rate following Tuna & Metzger (2023), and GE23 retains a + * fixed fraction of the requested super-Eddington rate. + * + * double ApplyEddingtonLimit(double p_RequestedRate, + * double p_EddingtonRate, + * CIRCUMBINARY_DISK_EDDINGTON_MODE p_Mode, + * double p_CapFactor, + * double p_Ge23Factor) + * + * @param [IN] p_RequestedRate Accretion rate requested by the Spindler mass split + * @param [IN] p_EddingtonRate Eddington accretion rate in matching units + * @param [IN] p_Mode Eddington-limiting prescription + * @param [IN] p_CapFactor Multiple of the Eddington rate allowed in CAP mode + * @param [IN] p_Ge23Factor Retained fraction above Eddington in GE23 mode + * @return Retained accretion rate + */ +double SpindlerBackend::ApplyEddingtonLimit( + double p_RequestedRate, + double p_EddingtonRate, + CIRCUMBINARY_DISK_EDDINGTON_MODE p_Mode, + double p_CapFactor, + double p_Ge23Factor +) { + switch (p_Mode) { + case CIRCUMBINARY_DISK_EDDINGTON_MODE::CAP: { + return std::min(p_RequestedRate, p_CapFactor * p_EddingtonRate); + } + + case CIRCUMBINARY_DISK_EDDINGTON_MODE::TM23: { + if (std::isfinite(p_EddingtonRate) && utils::Compare(p_EddingtonRate, 0.0) > 0 && utils::Compare(p_RequestedRate, p_EddingtonRate) > 0) { + const double ratio = p_RequestedRate / p_EddingtonRate; + const double tm23RetentionFraction = std::min(1.0, std::pow(CBD_TM23_RETENTION_COEFFICIENT / ratio, CBD_TM23_RETENTION_EXPONENT)); + return p_RequestedRate * tm23RetentionFraction; + } + return p_RequestedRate; + } + + case CIRCUMBINARY_DISK_EDDINGTON_MODE::GE23: { + return utils::Compare(p_RequestedRate, p_EddingtonRate) > 0 ? p_RequestedRate * p_Ge23Factor : p_RequestedRate; + } + } + + throw std::invalid_argument("eddingtonMode must be 'cap', 'ge23', or 'tm23'"); +} + +/* + * Calculate dimensionless orbital angular momentum in the Spindler normalisation + * + * double SpindlerOrbitalAngularMomentum(double p_Mass1Fraction, + * double p_Mass2Fraction, + * double p_Eccentricity, + * double p_SemiMajorAxisFactor) + * + * @param [IN] p_Mass1Fraction Ordered primary mass in Spindler units + * @param [IN] p_Mass2Fraction Ordered secondary mass in Spindler units + * @param [IN] p_Eccentricity Current eccentricity + * @param [IN] p_SemiMajorAxisFactor Current semi-major axis divided by its initial value + * @return Orbital angular momentum for G = Mtot,0 = a0 = 1 + */ +double SpindlerBackend::SpindlerOrbitalAngularMomentum(double p_Mass1Fraction, double p_Mass2Fraction, double p_Eccentricity, double p_SemiMajorAxisFactor) { + const double totalMassFraction = p_Mass1Fraction + p_Mass2Fraction; + const double reducedMassFraction = p_Mass1Fraction * p_Mass2Fraction / totalMassFraction; + const double oneMinusE2 = std::max(1.0 - p_Eccentricity * p_Eccentricity, std::numeric_limits::min()); + return reducedMassFraction * std::sqrt(totalMassFraction * p_SemiMajorAxisFactor * oneMinusE2); +} + +/* + * Evolve the Spindler system with Eddington-limited accretion + * + * State vector: [a/a0, M1/M0, M2/M0, e] in the ordered Spindler labels. The + * independent variable is time in years. The supplied CBD mass is spread + * uniformly over p_Params.durationYears, then split between the two stars using + * the Spindler qdot-derived accretion fractions. The encompassing CBD model assumes the initial + * reservoir, Initial_CBD_Mass, satisfies m_CBD << M_binary. + * + * The orbital terms have two pieces: + * - CBD torque/accretion response from the Spindler Da and De tables. + * - Extra non-conservative widening/shrinking from angular momentum carried + * away by material that was supplied but ejected due to Eddington limits. + * + * The table torques scale with the supplied mass flux. Eccentricity is evolved + * only through the tabulated CBD response; the simple wind-loss term affects a + * through Jdot but does not add an independent de/dt prescription. + * + * SpindlerEvolutionResult Evolve(const SpindlerEvolutionParams& p_Params) + * + * @param [IN] p_Params Ordered, dimensionless initial state and model parameters + * @return State reached at the requested duration + */ +SpindlerEvolutionResult SpindlerBackend::Evolve(const SpindlerEvolutionParams& p_Params) { + ValidateCommonInputs(p_Params.initialMassRatio, p_Params.initialEccentricity); + + if (utils::Compare(p_Params.suppliedMassFraction, 0.0) < 0 || !std::isfinite(p_Params.suppliedMassFraction)) { + throw std::invalid_argument("suppliedMassFraction must be non-negative and finite"); + } + if (utils::Compare(p_Params.initialBinaryMassMsol, 0.0) <= 0 || !std::isfinite(p_Params.initialBinaryMassMsol)) { + throw std::invalid_argument("initialBinaryMassMsol must be positive and finite"); + } + if (utils::Compare(p_Params.durationYears, 0.0) <= 0 || !std::isfinite(p_Params.durationYears)) { + throw std::invalid_argument("duration must be positive and finite"); + } + if (utils::Compare(p_Params.options.eddingtonCapFactor, 0.0) <= 0 || !std::isfinite(p_Params.options.eddingtonCapFactor)) { + throw std::invalid_argument("eddingtonCapFactor must be positive and finite"); + } + if (!(utils::Compare(p_Params.options.eddingtonGe23Factor, 0.0) > 0 && utils::Compare(p_Params.options.eddingtonGe23Factor, 1.0) <= 0) || + !std::isfinite(p_Params.options.eddingtonGe23Factor)) { + throw std::invalid_argument("eddingtonGe23Factor must satisfy 0 < f <= 1"); + } + + using StateType = std::vector; + StateType state = {1.0, 1.0 / (1.0 + p_Params.initialMassRatio), p_Params.initialMassRatio / (1.0 + p_Params.initialMassRatio), p_Params.initialEccentricity}; + + if (utils::Compare(p_Params.suppliedMassFraction, 0.0) == 0) { + return MakeResultFromState(state); + } + + const double supplyRateFractionPerYear = p_Params.suppliedMassFraction / p_Params.durationYears; + + // Define the ODE derivatives for orbital size, component masses, and + // eccentricity under CBD torques plus non-conservative mass loss. + auto system = [this, &p_Params, supplyRateFractionPerYear](const StateType& p_XState, StateType& p_DxDt, const double /*p_Time*/) { + p_DxDt.resize(4); + const double semiMajorAxisFactor = p_XState[0]; + const double mass1Fraction = p_XState[1]; + const double mass2Fraction = p_XState[2]; + const double eccentricity = p_XState[3]; + + if (utils::Compare(semiMajorAxisFactor, 0.0) <= 0 || utils::Compare(mass1Fraction, 0.0) <= 0 || utils::Compare(mass2Fraction, 0.0) <= 0) { + throw std::runtime_error("Non-positive a, M1, or M2 encountered during Spindler integration"); + } + + const double totalMassFraction = mass1Fraction + mass2Fraction; + const double massRatio = mass2Fraction / mass1Fraction; + const double orbitalAngularMomentum = SpindlerOrbitalAngularMomentum(mass1Fraction, mass2Fraction, eccentricity, semiMajorAxisFactor); + + const double mdotEdd1 = ComputeEddingtonRate( + p_Params.stellarType1, mass1Fraction, p_Params.initialBinaryMassMsol, p_Params.radius1Rsol + ); + const double mdotEdd2 = ComputeEddingtonRate( + p_Params.stellarType2, mass2Fraction, p_Params.initialBinaryMassMsol, p_Params.radius2Rsol + ); + + const auto componentMassFractions = ComponentMassFractions(massRatio, eccentricity); + const double mdotReq1 = componentMassFractions.first * supplyRateFractionPerYear; + const double mdotReq2 = componentMassFractions.second * supplyRateFractionPerYear; + + const double mdotAcc1 = ApplyEddingtonLimit( + mdotReq1, + mdotEdd1, + p_Params.options.eddingtonMode, + p_Params.options.eddingtonCapFactor, + p_Params.options.eddingtonGe23Factor + ); + const double mdotAcc2 = ApplyEddingtonLimit( + mdotReq2, + mdotEdd2, + p_Params.options.eddingtonMode, + p_Params.options.eddingtonCapFactor, + p_Params.options.eddingtonGe23Factor + ); + const double mdotLost1 = std::max(mdotReq1 - mdotAcc1, 0.0); + const double mdotLost2 = std::max(mdotReq2 - mdotAcc2, 0.0); + + // The Spindler tables give the CBD-driven response per unit supplied mass + // through the disk/binary interaction. + const auto derivatives = DerivativesForLabeledQ(massRatio, eccentricity); + const double dadtDisk = derivatives.da * semiMajorAxisFactor * supplyRateFractionPerYear / totalMassFraction; + const double dedt = derivatives.de * eccentricity * supplyRateFractionPerYear / totalMassFraction; + + const double jBinary = orbitalAngularMomentum / totalMassFraction; + const double j1 = massRatio * jBinary; + const double j2 = jBinary / massRatio; + const double dJdtWind = -(j1 * mdotLost1 + j2 * mdotLost2); + const double dadtWind = utils::Compare(orbitalAngularMomentum, 0.0) > 0 ? 2.0 * semiMajorAxisFactor * dJdtWind / orbitalAngularMomentum : 0.0; + + p_DxDt[0] = dadtDisk + dadtWind; + p_DxDt[1] = mdotAcc1; + p_DxDt[2] = mdotAcc2; + p_DxDt[3] = dedt; + }; + + auto stepper = odeint::make_controlled( + CBD_SPINDLER_ODE_ABSOLUTE_TOLERANCE, + CBD_SPINDLER_ODE_RELATIVE_TOLERANCE, + odeint::runge_kutta_cash_karp54() + ); + + double time = 0.0; + double dt = std::max(p_Params.durationYears / 1000.0, 1.0e-8); + + // Deliberately use strict comparisons here rather than utils::Compare(): + // the ODE controller must advance to, and clip exactly at, the requested + // endpoint instead of treating a nearby time as equivalent. + while (time < p_Params.durationYears) { + if (time + dt > p_Params.durationYears) { + dt = p_Params.durationYears - time; + } + + boost::numeric::odeint::controlled_step_result stepResult; + do { + stepResult = stepper.try_step(system, state, time, dt); + } while (stepResult == boost::numeric::odeint::fail); + } + + return MakeResultFromState(state); +} + +/* + * Translate COMPAS-labelled input to Spindler variables, evolve it, and translate the result back + * + * CircumbinaryDisk::EvolutionTrack EvolveTrackWithBackend(const CircumbinaryDisk::EvolutionInput& p_Input, + * SpindlerBackend& p_Spindler) + * + * @param [IN] p_Input COMPAS-labelled physical input for one binary-CBD interaction segment + * @param [IN] p_Spindler Initialised Spindler table backend + * @return Final COMPAS-labelled state after the requested increment + */ +CircumbinaryDisk::EvolutionTrack EvolveTrackWithBackend( + const CircumbinaryDisk::EvolutionInput& p_Input, + SpindlerBackend& p_Spindler +) { + const double binaryMassMsol = p_Input.mass1Msol + p_Input.mass2Msol; + const bool componentsSwapped = utils::Compare(p_Input.mass1Msol, p_Input.mass2Msol) < 0; + const double primaryMassMsol = componentsSwapped ? p_Input.mass2Msol : p_Input.mass1Msol; + const double secondaryMassMsol = componentsSwapped ? p_Input.mass1Msol : p_Input.mass2Msol; + + // Convert the public COMPAS-labelled, physical-unit input into the ordered, + // dimensionless variables expected by the Spindler ODE. componentsSwapped + // records whether this translation exchanged the COMPAS labels, so the + // output can be mapped back to the original star ordering. + SpindlerEvolutionParams params{}; + params.suppliedMassFraction = p_Input.suppliedMassMsol / binaryMassMsol; + params.initialMassRatio = secondaryMassMsol / primaryMassMsol; + params.initialEccentricity = p_Input.eccentricity; + params.initialBinaryMassMsol = binaryMassMsol; + params.durationYears = p_Input.durationYears; + params.stellarType1 = componentsSwapped ? p_Input.stellarType2 : p_Input.stellarType1; + params.stellarType2 = componentsSwapped ? p_Input.stellarType1 : p_Input.stellarType2; + params.radius1Rsol = componentsSwapped ? p_Input.radius2Rsol : p_Input.radius1Rsol; + params.radius2Rsol = componentsSwapped ? p_Input.radius1Rsol : p_Input.radius2Rsol; + params.options = p_Input.options; + + const SpindlerEvolutionResult spindlerResult{p_Spindler.Evolve(params)}; + + CircumbinaryDisk::EvolutionTrack track{}; + track.semiMajorAxisFactor = spindlerResult.semiMajorAxisFactor; + track.eccentricity = spindlerResult.eccentricity; + track.mass1Msol = (componentsSwapped ? spindlerResult.mass2Fraction : spindlerResult.mass1Fraction) * binaryMassMsol; + track.mass2Msol = (componentsSwapped ? spindlerResult.mass1Fraction : spindlerResult.mass2Fraction) * binaryMassMsol; + + return track; +} + +} // namespace + +namespace CircumbinaryDisk { + +/* + * Evolve a formed CBD and return COMPAS-labelled final state + * + * spindler-c is initialized lazily on the first CBD call in a COMPAS process + * and then reused for all subsequent binaries in that process. Spindler evolves + * a/a0, so the returned semi-major-axis value is a multiplicative factor that + * the caller applies to its current binary separation. + * + * EvolutionTrack EvolveTrack(const EvolutionInput& p_Input) + * + * @param [IN] p_Input COMPAS-labelled physical input; per-call supplied mass is assumed to be << binary mass + * @return Final COMPAS-labelled state after the requested increment + */ +EvolutionTrack EvolveTrack(const EvolutionInput& p_Input) { + static SpindlerBackend spindler{CBD_SPINDLER_MODEL_NAME}; + return EvolveTrackWithBackend(p_Input, spindler); +} + + +} // namespace CircumbinaryDisk diff --git a/src/CircumbinaryDiskEvolutionUnavailable.cpp b/src/CircumbinaryDiskEvolutionUnavailable.cpp new file mode 100644 index 000000000..657e73073 --- /dev/null +++ b/src/CircumbinaryDiskEvolutionUnavailable.cpp @@ -0,0 +1,24 @@ +#include "CircumbinaryDisk.h" + +// Stub implementation used when COMPAS is built without spindler-c support. + +#include + +namespace CircumbinaryDisk { + +/* + * Reject CBD evolution when COMPAS was built without the spindler-c backend + * + * EvolutionTrack EvolveTrack(const EvolutionInput&) + * + * @return This function does not return; it throws a runtime error + */ +EvolutionTrack EvolveTrack(const EvolutionInput&) { + throw std::runtime_error( + "Circumbinary-disk evolution was requested, but this COMPAS executable " + "was built without spindler-c support. Rebuild with SPINDLER_C_ROOT " + "set to the spindler-c build tree." + ); +} + +} // namespace CircumbinaryDisk diff --git a/src/CircumbinaryDiskFormation.cpp b/src/CircumbinaryDiskFormation.cpp new file mode 100644 index 000000000..a319991b5 --- /dev/null +++ b/src/CircumbinaryDiskFormation.cpp @@ -0,0 +1,158 @@ +#include "CircumbinaryDisk.h" + +#include "constants.h" +#include "Options.h" +#include "utils.h" + +#include + +namespace CircumbinaryDisk { + + +/* + * Copy the current global COMPAS options needed by the CBD formation model + * + * FormationOptions MakeFormationOptionsFromCompasOptions() + * + * @return CBD formation options + */ +FormationOptions MakeFormationOptionsFromCompasOptions() { + return { + OPTIONS->CircumbinaryDiskBeta(), + OPTIONS->CircumbinaryDiskEnvelopeMassFractionSupplied(), + OPTIONS->CircumbinaryDiskStructureFactor(), + OPTIONS->CircumbinaryDiskInnerRadiusOverSeparation() + }; +} + +/* + * Copy the current global COMPAS options needed by the CBD evolution model + * + * EvolutionOptions MakeEvolutionOptionsFromCompasOptions() + * + * @return CBD evolution options + */ +EvolutionOptions MakeEvolutionOptionsFromCompasOptions() { + return { + OPTIONS->CircumbinaryDiskEddingtonMode(), + OPTIONS->CircumbinaryDiskEddingtonCapFactor(), + OPTIONS->CircumbinaryDiskEddingtonGe23Factor() + }; +} + +/* + * Determine whether a set of supernova-like events destroys an active CBD + * + * AIC is treated as a change to the live COMPAS stellar state, so a bound + * system continues evolving through the same physical disk. + * + * bool SupernovaEventsDisruptDisk(const SN_EVENT p_Events) + * + * @param [IN] p_Events Supernova-event bitmask generated during the current event + * @return true for any non-AIC supernova event; false otherwise + */ +bool SupernovaEventsDisruptDisk(const SN_EVENT p_Events) { + return p_Events != SN_EVENT::NONE && p_Events != SN_EVENT::AIC; +} + +/* + * Determine whether a CE has the supported compact-object plus envelope-donor geometry + * + * bool IsSupportedCommonEnvelopeChannel(const STELLAR_TYPE p_StellarType1PreCE, + * const STELLAR_TYPE p_StellarType2PreCE, + * const double p_Mass1PreCE, + * const double p_Mass2PreCE, + * const bool p_EnvelopeFlag1, + * const bool p_EnvelopeFlag2) + * + * @param [IN] p_StellarType1PreCE Stellar type of star 1 immediately before CE + * @param [IN] p_StellarType2PreCE Stellar type of star 2 immediately before CE + * @param [IN] p_Mass1PreCE Mass of star 1 immediately before CE (Msol) + * @param [IN] p_Mass2PreCE Mass of star 2 immediately before CE (Msol) + * @param [IN] p_EnvelopeFlag1 Whether star 1 has an envelope immediately before CE + * @param [IN] p_EnvelopeFlag2 Whether star 2 has an envelope immediately before CE + * @return Whether the CE channel is supported + */ +bool IsSupportedCommonEnvelopeChannel(const STELLAR_TYPE p_StellarType1PreCE, + const STELLAR_TYPE p_StellarType2PreCE, + const double p_Mass1PreCE, + const double p_Mass2PreCE, + const bool p_EnvelopeFlag1, + const bool p_EnvelopeFlag2) { + + const STELLAR_TYPE_LIST supportedCompactObjects = { + STELLAR_TYPE::HELIUM_WHITE_DWARF, + STELLAR_TYPE::CARBON_OXYGEN_WHITE_DWARF, + STELLAR_TYPE::OXYGEN_NEON_WHITE_DWARF, + STELLAR_TYPE::NEUTRON_STAR, + STELLAR_TYPE::BLACK_HOLE + }; + const bool star1IsCompactObject = utils::IsOneOf(p_StellarType1PreCE, supportedCompactObjects); + const bool star2IsCompactObject = utils::IsOneOf(p_StellarType2PreCE, supportedCompactObjects); + + if (!((star1IsCompactObject && p_EnvelopeFlag2 && !p_EnvelopeFlag1) || + (star2IsCompactObject && p_EnvelopeFlag1 && !p_EnvelopeFlag2))) { + return false; + } + + const double binaryMassPreCEMsol = p_Mass1PreCE + p_Mass2PreCE; + if (utils::Compare(binaryMassPreCEMsol, 0.0) <= 0) { + return false; + } + + return true; +} + +/* + * Evaluate CBD formation criterion + * + * The mass supplied to the potential CBD is a fraction f_M of the decrease in total binary mass + * across CE: (M1 + M2)_pre-CE - (M1 + M2)_post-CE. The resulting initial CBD + * reservoir is assumed to satisfy m_CBD << M_binary during the evolution; formation does not + * check this assumption. + * + * The criterion compares the total angular momentum assigned to the + * CBD reservoir against the total angular momentum required for + * supporting the CBD, given by the Keplerian AM at the adopted CBD inner radius times a disk structure factor. + * + * FormationResult EvaluateFormation(const FormationInput& p_Input) + * + * @param [IN] p_Input Physical state and model options for the formation criterion + * @return Formation outcome and supplied mass + */ +FormationResult EvaluateFormation(const FormationInput& p_Input) { + FormationResult result{}; + + const double massLostDuringCEMsol = p_Input.binaryMassPreCEMsol - p_Input.binaryMassPostCEMsol; + if (utils::Compare(massLostDuringCEMsol, 0.0) <= 0 || utils::Compare(p_Input.semiMajorAxisPreRLOFRsol, 0.0) <= 0) { + return result; + } + + result.initialReservoirMassMsol = p_Input.options.envelopeMassFractionSupplied * massLostDuringCEMsol; + if (utils::Compare(result.initialReservoirMassMsol, 0.0) <= 0) { + return result; + } + + const double totalAngularMomentumPreRLOF = p_Input.totalAngularMomentumPreRLOF; + // beta = fJ/fM, so the angular momentum assigned to the supplied CBD + // reservoir is J_CBD = beta * fM * J_available. + const double suppliedDiskAngularMomentum = p_Input.options.circumbinaryDiskBeta * + p_Input.options.envelopeMassFractionSupplied * + totalAngularMomentumPreRLOF; + + // The adopted inner CBD radius scales with the post-CE separation because + // the disk is assumed to orbit the post-CE binary, not the + // larger pre-CE orbit. Formation is allowed when the assigned angular + // momentum meets or exceeds the threshold. + const double diskInnerRadiusAU = p_Input.options.diskInnerRadiusOverSeparation * p_Input.semiMajorAxisPostCERsol * RSOL_TO_AU; + const double minimumDiskAngularMomentum = result.initialReservoirMassMsol * + p_Input.options.diskStructureFactor * + std::sqrt(G_AU_Msol_yr * + p_Input.binaryMassPostCEMsol * + diskInnerRadiusAU); + + result.forms = utils::Compare(suppliedDiskAngularMomentum, minimumDiskAngularMomentum) >= 0; + return result; +} + +} // namespace CircumbinaryDisk diff --git a/src/Log.cpp b/src/Log.cpp index ca4e3824c..a751586bf 100644 --- a/src/Log.cpp +++ b/src/Log.cpp @@ -285,6 +285,7 @@ void Log::Start(const string p_LogBasePathString, // BSE if (NotesPropertyPresent(m_BSE_CEE_Rec )) m_BSE_CEE_Notes = BOOL_VECTOR(m_BSE_CEE_Notes.size(), true); + if (NotesPropertyPresent(m_BSE_CBD_Rec )) m_BSE_CBD_Notes = BOOL_VECTOR(m_BSE_CBD_Notes.size(), true); if (NotesPropertyPresent(m_BSE_DCO_Rec )) m_BSE_DCO_Notes = BOOL_VECTOR(m_BSE_DCO_Notes.size(), true); if (NotesPropertyPresent(m_BSE_Detailed_Rec )) m_BSE_Detailed_Notes = BOOL_VECTOR(m_BSE_Detailed_Notes.size(), true); if (NotesPropertyPresent(m_BSE_Pulsars_Rec )) m_BSE_Pulsars_Notes = BOOL_VECTOR(m_BSE_Pulsars_Notes.size(), true); @@ -2071,6 +2072,11 @@ std::tuple Log::GetStandardLogFile annotations = m_BSE_CEE_Notes; // logfile annotations break; + case LOGFILE::BSE_CIRCUMBINARY_DISKS: // BSE_CIRCUMBINARY_DISKS + recordProperties = m_BSE_CBD_Rec; // record properties + annotations = m_BSE_CBD_Notes; // logfile annotations + break; + case LOGFILE::BSE_DETAILED_OUTPUT: // BSE_DETAILED_OUTPUT recordProperties = m_BSE_Detailed_Rec; // record properties annotations = m_BSE_Detailed_Notes; // logfile annotations @@ -2440,6 +2446,13 @@ LogfileDetailsT Log::StandardLogFileDetails(const LOGFILE p_Logfile, const strin fileDetails.annotations = m_BSE_CEE_Notes; break; + case LOGFILE::BSE_CIRCUMBINARY_DISKS: // BSE_CIRCUMBINARY_DISKS + fileDetails.filename = std::get<0>(LOGFILE_DESCRIPTOR.at(LOGFILE::BSE_CIRCUMBINARY_DISKS)); + fileDetails.recordTypes = {1U}; + fileDetails.recordProperties = m_BSE_CBD_Rec; + fileDetails.annotations = m_BSE_CBD_Notes; + break; + case LOGFILE::BSE_DOUBLE_COMPACT_OBJECTS: // BSE_DOUBLE_COMPACT_OBJECTS fileDetails.filename = OPTIONS->LogfileDoubleCompactObjects(); fileDetails.recordTypes = OPTIONS->LogfileDoubleCompactObjectsRecordTypes(); @@ -3147,6 +3160,10 @@ void Log::UpdateLogfileRecordSpecs(const LOGFILE p_Logfile, if (p_UseDefaultProps) baseProps = m_BSE_CEE_Rec; baseNotes = m_BSE_CEE_Notes; break; + case LOGFILE::BSE_CIRCUMBINARY_DISKS: + if (p_UseDefaultProps) baseProps = m_BSE_CBD_Rec; + baseNotes = m_BSE_CBD_Notes; + break; case LOGFILE::BSE_DETAILED_OUTPUT: if (p_UseDefaultProps) baseProps = m_BSE_Detailed_Rec; baseNotes = m_BSE_Detailed_Notes; @@ -3294,6 +3311,7 @@ void Log::UpdateLogfileRecordSpecs(const LOGFILE p_Logfile, // replace existing props and annotations vector for given logfile switch (p_Logfile) { case LOGFILE::BSE_COMMON_ENVELOPES : m_BSE_CEE_Rec = newProps; m_BSE_CEE_Notes = newNotes; break; + case LOGFILE::BSE_CIRCUMBINARY_DISKS : m_BSE_CBD_Rec = newProps; m_BSE_CBD_Notes = newNotes; break; case LOGFILE::BSE_DETAILED_OUTPUT : m_BSE_Detailed_Rec = newProps; m_BSE_Detailed_Notes = newNotes; break; case LOGFILE::BSE_DOUBLE_COMPACT_OBJECTS: m_BSE_DCO_Rec = newProps; m_BSE_DCO_Notes = newNotes; break; case LOGFILE::BSE_PULSAR_EVOLUTION : m_BSE_Pulsars_Rec = newProps; m_BSE_Pulsars_Notes = newNotes; break; diff --git a/src/Log.h b/src/Log.h index 3ca691105..7525aa21a 100755 --- a/src/Log.h +++ b/src/Log.h @@ -570,6 +570,7 @@ class Log { // logfile record specifications // BSE ANY_PROPERTY_VECTOR m_BSE_CEE_Rec = BSE_COMMON_ENVELOPES_REC; // default specification + ANY_PROPERTY_VECTOR m_BSE_CBD_Rec = BSE_CIRCUMBINARY_DISKS_REC; // default specification ANY_PROPERTY_VECTOR m_BSE_DCO_Rec = BSE_DOUBLE_COMPACT_OBJECTS_REC; // default specification ANY_PROPERTY_VECTOR m_BSE_Detailed_Rec = BSE_DETAILED_OUTPUT_REC; // default specification ANY_PROPERTY_VECTOR m_BSE_Pulsars_Rec = BSE_PULSAR_EVOLUTION_REC; // default specification @@ -603,6 +604,7 @@ class Log { // BSE BOOL_VECTOR m_BSE_CEE_Notes = BOOL_VECTOR(OPTIONS->NotesHdrs().size(), false); + BOOL_VECTOR m_BSE_CBD_Notes = BOOL_VECTOR(OPTIONS->NotesHdrs().size(), false); BOOL_VECTOR m_BSE_DCO_Notes = BOOL_VECTOR(OPTIONS->NotesHdrs().size(), false); BOOL_VECTOR m_BSE_Detailed_Notes = BOOL_VECTOR(OPTIONS->NotesHdrs().size(), false); BOOL_VECTOR m_BSE_Pulsars_Notes = BOOL_VECTOR(OPTIONS->NotesHdrs().size(), false); @@ -1244,6 +1246,10 @@ class Log { bool LogCommonEnvelope(const T* const p_Binary, const CE_RECORD_TYPE p_RecordType) { return LogStandardRecord(std::get<2>(LOGFILE_DESCRIPTOR.at(LOGFILE::BSE_COMMON_ENVELOPES)), 0, LOGFILE::BSE_COMMON_ENVELOPES, static_cast(p_RecordType), p_Binary); } + // Write one circumbinary-disk event to the standard CBD logfile. + template + bool LogCircumbinaryDisk(const T* const p_Binary) { return LogStandardRecord(std::get<2>(LOGFILE_DESCRIPTOR.at(LOGFILE::BSE_CIRCUMBINARY_DISKS)), 0, LOGFILE::BSE_CIRCUMBINARY_DISKS, 1U, p_Binary); } + template bool LogDoubleCompactObject(const T* const p_Binary, const DCO_RECORD_TYPE p_RecordType) { return LogStandardRecord(std::get<2>(LOGFILE_DESCRIPTOR.at(LOGFILE::BSE_DOUBLE_COMPACT_OBJECTS)), 0, LOGFILE::BSE_DOUBLE_COMPACT_OBJECTS, static_cast(p_RecordType), p_Binary); } diff --git a/src/LogTypedefs.h b/src/LogTypedefs.h index cb99c198f..26e5f8d4b 100644 --- a/src/LogTypedefs.h +++ b/src/LogTypedefs.h @@ -477,6 +477,20 @@ enum class ANY_STAR_PROPERTY: int { STAR_PROPERTIES }; enum class BINARY_PROPERTY: int { NONE, CIRCULARIZATION_TIMESCALE, + CIRCUMBINARY_DISK_FORMED, + CIRCUMBINARY_DISK_FORMATION_TIME, + CIRCUMBINARY_DISK_MASS_1_AT_FORMATION, + CIRCUMBINARY_DISK_MASS_2_AT_FORMATION, + CIRCUMBINARY_DISK_STELLAR_TYPE_1_AT_FORMATION, + CIRCUMBINARY_DISK_STELLAR_TYPE_2_AT_FORMATION, + CIRCUMBINARY_DISK_INITIAL_MASS, + CIRCUMBINARY_DISK_MASS_RETAINED_1, + CIRCUMBINARY_DISK_MASS_RETAINED_2, + CIRCUMBINARY_DISK_MASS_SUPPLIED, + CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_PRE, + CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_POST, + CIRCUMBINARY_DISK_ECCENTRICITY_PRE, + CIRCUMBINARY_DISK_ECCENTRICITY_POST, COMMON_ENVELOPE_AT_LEAST_ONCE, COMMON_ENVELOPE_EVENT_COUNT, UNBOUND, @@ -635,6 +649,20 @@ enum class BINARY_PROPERTY: int { const COMPASUnorderedMap BINARY_PROPERTY_LABEL = { { BINARY_PROPERTY::NONE, "NONE" }, { BINARY_PROPERTY::CIRCULARIZATION_TIMESCALE, "CIRCULARIZATION_TIMESCALE" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_FORMED, "CIRCUMBINARY_DISK_FORMED" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_FORMATION_TIME, "CIRCUMBINARY_DISK_FORMATION_TIME" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_1_AT_FORMATION, "CIRCUMBINARY_DISK_MASS_1_AT_FORMATION" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_2_AT_FORMATION, "CIRCUMBINARY_DISK_MASS_2_AT_FORMATION" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_STELLAR_TYPE_1_AT_FORMATION, "CIRCUMBINARY_DISK_STELLAR_TYPE_1_AT_FORMATION" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_STELLAR_TYPE_2_AT_FORMATION, "CIRCUMBINARY_DISK_STELLAR_TYPE_2_AT_FORMATION" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_INITIAL_MASS, "CIRCUMBINARY_DISK_INITIAL_MASS" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_RETAINED_1, "CIRCUMBINARY_DISK_MASS_RETAINED_1" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_RETAINED_2, "CIRCUMBINARY_DISK_MASS_RETAINED_2" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_SUPPLIED, "CIRCUMBINARY_DISK_MASS_SUPPLIED" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_PRE, "CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_PRE" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_POST, "CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_POST" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_ECCENTRICITY_PRE, "CIRCUMBINARY_DISK_ECCENTRICITY_PRE" }, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_ECCENTRICITY_POST, "CIRCUMBINARY_DISK_ECCENTRICITY_POST" }, { BINARY_PROPERTY::COMMON_ENVELOPE_AT_LEAST_ONCE, "COMMON_ENVELOPE_AT_LEAST_ONCE" }, { BINARY_PROPERTY::COMMON_ENVELOPE_EVENT_COUNT, "COMMON_ENVELOPE_EVENT_COUNT" }, { BINARY_PROPERTY::UNBOUND, "UNBOUND" }, @@ -814,6 +842,19 @@ enum class PROGRAM_OPTION: int { CIRCULARISE_BINARY_DURING_MT, + CIRCUMBINARY_DISK, + CIRCUMBINARY_DISK_BETA, + CIRCUMBINARY_DISK_EDDINGTON_CAP_FACTOR, + CIRCUMBINARY_DISK_EDDINGTON_GE23_FACTOR, + CIRCUMBINARY_DISK_EDDINGTON_MODE, + CIRCUMBINARY_DISK_INNER_RADIUS_OVER_SEPARATION, + CIRCUMBINARY_DISK_LIFETIME, + CIRCUMBINARY_DISK_STRUCTURE_FACTOR, + CIRCUMBINARY_DISK_ENVELOPE_MASS_FRACTION_SUPPLIED, + POST_COMMON_ENVELOPE_ECCENTRICITY_CAP, + POST_COMMON_ENVELOPE_ECCENTRICITY_FRACTION, + POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION, + COMMON_ENVELOPE_ALPHA, COMMON_ENVELOPE_ALPHA_THERMAL, COMMON_ENVELOPE_FORMALISM, @@ -1051,6 +1092,19 @@ const COMPASUnorderedMap PROGRAM_OPTION_LABEL = { { PROGRAM_OPTION::CIRCULARISE_BINARY_DURING_MT, "CIRCULARISE_BINARY_DURING_MT" }, + { PROGRAM_OPTION::CIRCUMBINARY_DISK, "CIRCUMBINARY_DISK" }, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_BETA, "CIRCUMBINARY_DISK_BETA" }, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_EDDINGTON_CAP_FACTOR, "CIRCUMBINARY_DISK_EDDINGTON_CAP_FACTOR" }, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_EDDINGTON_GE23_FACTOR, "CIRCUMBINARY_DISK_EDDINGTON_GE23_FACTOR" }, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_EDDINGTON_MODE, "CIRCUMBINARY_DISK_EDDINGTON_MODE" }, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_INNER_RADIUS_OVER_SEPARATION, "CIRCUMBINARY_DISK_INNER_RADIUS_OVER_SEPARATION" }, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_LIFETIME, "CIRCUMBINARY_DISK_LIFETIME" }, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_STRUCTURE_FACTOR, "CIRCUMBINARY_DISK_STRUCTURE_FACTOR" }, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_ENVELOPE_MASS_FRACTION_SUPPLIED, "CIRCUMBINARY_DISK_ENVELOPE_MASS_FRACTION_SUPPLIED" }, + { PROGRAM_OPTION::POST_COMMON_ENVELOPE_ECCENTRICITY_CAP, "POST_COMMON_ENVELOPE_ECCENTRICITY_CAP" }, + { PROGRAM_OPTION::POST_COMMON_ENVELOPE_ECCENTRICITY_FRACTION, "POST_COMMON_ENVELOPE_ECCENTRICITY_FRACTION" }, + { PROGRAM_OPTION::POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION, "POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION" }, + { PROGRAM_OPTION::COMMON_ENVELOPE_ALPHA, "COMMON_ENVELOPE_ALPHA" }, { PROGRAM_OPTION::COMMON_ENVELOPE_ALPHA_THERMAL, "COMMON_ENVELOPE_ALPHA_THERMAL" }, { PROGRAM_OPTION::COMMON_ENVELOPE_FORMALISM, "COMMON_ENVELOPE_FORMALISM" }, @@ -1463,6 +1517,20 @@ const std::map ANY_STAR_PROPERTY_DETAIL = { // the logfiles - all keys present here should also be in BINARY_PROPERTY and BINARY_PROPERTY_LABEL const std::map BINARY_PROPERTY_DETAIL = { { BINARY_PROPERTY::CIRCULARIZATION_TIMESCALE, { TYPENAME::DOUBLE, "Tau_Circ", "Myr", 24, 15}}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_FORMED, { TYPENAME::BOOL, "CBD", "Event", 0, 0 }}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_FORMATION_TIME, { TYPENAME::DOUBLE, "Time", "Myr", 24, 15}}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_1_AT_FORMATION, { TYPENAME::DOUBLE, "Mass(1)", "Msol", 24, 15}}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_2_AT_FORMATION, { TYPENAME::DOUBLE, "Mass(2)", "Msol", 24, 15}}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_STELLAR_TYPE_1_AT_FORMATION, { TYPENAME::STELLAR_TYPE, "Stellar_Type(1)", "-", 4, 1 }}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_STELLAR_TYPE_2_AT_FORMATION, { TYPENAME::STELLAR_TYPE, "Stellar_Type(2)", "-", 4, 1 }}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_INITIAL_MASS, { TYPENAME::DOUBLE, "Initial_CBD_Mass", "Msol", 24, 15}}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_RETAINED_1, { TYPENAME::DOUBLE, "CBD_Mass_Retained(1)", "Msol", 24, 15}}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_RETAINED_2, { TYPENAME::DOUBLE, "CBD_Mass_Retained(2)", "Msol", 24, 15}}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_SUPPLIED, { TYPENAME::DOUBLE, "CBD_Mass_Supplied", "Msol", 24, 15}}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_PRE, { TYPENAME::DOUBLE, "SemiMajorAxisCBD", "Rsol", 24, 15}}, + { BINARY_PROPERTY::CIRCUMBINARY_DISK_ECCENTRICITY_PRE, { TYPENAME::DOUBLE, "EccentricityCBD", "-", 24, 15}}, { BINARY_PROPERTY::COMMON_ENVELOPE_AT_LEAST_ONCE, { TYPENAME::BOOL, "CEE", "Event", 0, 0 }}, { BINARY_PROPERTY::COMMON_ENVELOPE_EVENT_COUNT, { TYPENAME::UINT, "CE_Event_Counter", "Count", 6, 1 }}, { BINARY_PROPERTY::DOUBLE_CORE_COMMON_ENVELOPE, { TYPENAME::BOOL, "Double_Core_CE", "Event", 0, 0 }}, @@ -1642,6 +1710,19 @@ const std::map PROGRAM_OPTION_DETAIL = { { PROGRAM_OPTION::CIRCULARISE_BINARY_DURING_MT, { TYPENAME::BOOL, "PO_Circularise@MT", "Flag", 0, 0 }}, + { PROGRAM_OPTION::CIRCUMBINARY_DISK, { TYPENAME::BOOL, "PO_CBD", "Flag", 0, 0 }}, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_EDDINGTON_MODE, { TYPENAME::INT, "PO_CBD_Eddington_Mode", "-", 4, 1 }}, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_EDDINGTON_CAP_FACTOR, { TYPENAME::DOUBLE, "PO_CBD_Eddington_Cap_Factor", "-", 24, 15}}, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_EDDINGTON_GE23_FACTOR, { TYPENAME::DOUBLE, "PO_CBD_Eddington_GE23_Factor", "-", 24, 15}}, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_BETA, { TYPENAME::DOUBLE, "PO_CBD_Beta", "-", 24, 15}}, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_INNER_RADIUS_OVER_SEPARATION, { TYPENAME::DOUBLE, "PO_CBD_Inner_Radius_Over_Separation", "-", 24, 15}}, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_LIFETIME, { TYPENAME::DOUBLE, "PO_CBD_Lifetime", "yr", 24, 15}}, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_STRUCTURE_FACTOR, { TYPENAME::DOUBLE, "PO_CBD_Structure_Factor", "-", 24, 15}}, + { PROGRAM_OPTION::CIRCUMBINARY_DISK_ENVELOPE_MASS_FRACTION_SUPPLIED, { TYPENAME::DOUBLE, "PO_CBD_Envelope_Mass_Fraction_Supplied", "-", 24, 15}}, + { PROGRAM_OPTION::POST_COMMON_ENVELOPE_ECCENTRICITY_CAP, { TYPENAME::DOUBLE, "PO_Post_CE_Eccentricity_Cap", "-", 24, 15}}, + { PROGRAM_OPTION::POST_COMMON_ENVELOPE_ECCENTRICITY_FRACTION, { TYPENAME::DOUBLE, "PO_Post_CE_Eccentricity_Fraction", "-", 24, 15}}, + { PROGRAM_OPTION::POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION, { TYPENAME::INT, "PO_Post_CE_Eccentricity_Prescription", "-", 4, 1}}, + { PROGRAM_OPTION::COMMON_ENVELOPE_ALPHA, { TYPENAME::DOUBLE, "PO_CE_Alpha", "-", 24, 15}}, { PROGRAM_OPTION::COMMON_ENVELOPE_ALPHA_THERMAL, { TYPENAME::DOUBLE, "PO_CE_Alpha_Thermal", "-", 24, 15}}, { PROGRAM_OPTION::COMMON_ENVELOPE_FORMALISM, { TYPENAME::INT, "PO_CE_Formalism", "-", 4, 1 }}, @@ -1856,6 +1937,7 @@ enum class LOGFILE: int { ERROR_LOG, BSE_COMMON_ENVELOPES, + BSE_CIRCUMBINARY_DISKS, BSE_DETAILED_OUTPUT, BSE_DOUBLE_COMPACT_OBJECTS, BSE_PULSAR_EVOLUTION, @@ -2157,6 +2239,34 @@ const ANY_PROPERTY_VECTOR BSE_DETAILED_OUTPUT_REC = { }; +// BSE_CIRCUMBINARY_DISKS_REC +// +// Default record definition for the BSE Circumbinary Disks logfile. +// +// This file is event-based: one row is written for each CBD that forms. +// The row is written when the active disk is exhausted or truncated by an +// ordinary COMPAS termination/event boundary. Time, Mass(1/2), and +// Stellar_Type(1/2) are cached at CBD formation rather than sampled at write time. +// Initial_CBD_Mass is the disk +// reservoir defined at formation; CBD_Mass_Supplied is the amount actually +// processed before the row is written. +// +const ANY_PROPERTY_VECTOR BSE_CIRCUMBINARY_DISKS_REC = { + BINARY_PROPERTY::RANDOM_SEED, + BINARY_PROPERTY::CIRCUMBINARY_DISK_FORMATION_TIME, + BINARY_PROPERTY::CIRCUMBINARY_DISK_INITIAL_MASS, + BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_SUPPLIED, + BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_RETAINED_1, + BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_RETAINED_2, + BINARY_PROPERTY::CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_PRE, + BINARY_PROPERTY::CIRCUMBINARY_DISK_SEMI_MAJOR_AXIS_POST, + BINARY_PROPERTY::CIRCUMBINARY_DISK_ECCENTRICITY_PRE, + BINARY_PROPERTY::CIRCUMBINARY_DISK_ECCENTRICITY_POST, + BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_1_AT_FORMATION, + BINARY_PROPERTY::CIRCUMBINARY_DISK_STELLAR_TYPE_1_AT_FORMATION, + BINARY_PROPERTY::CIRCUMBINARY_DISK_MASS_2_AT_FORMATION, + BINARY_PROPERTY::CIRCUMBINARY_DISK_STELLAR_TYPE_2_AT_FORMATION +}; // BSE_DOUBLE_COMPACT_OBJECT_REC // // Default record definition for the Double Compact Objects logfile @@ -2601,6 +2711,7 @@ const std::map LOGFILE_DESCRIPTOR = { { LOGFILE::ERROR_LOG, { "Error_Log", {}, "", "", LOGFILE_CLASS::NONE }}, { LOGFILE::BSE_COMMON_ENVELOPES, { "BSE_Common_Envelopes", BSE_COMMON_ENVELOPES_REC, "BSE_CEE", "BSE_CEE_REC", LOGFILE_CLASS::BINARY }}, + { LOGFILE::BSE_CIRCUMBINARY_DISKS, { "BSE_Circumbinary_Disks", BSE_CIRCUMBINARY_DISKS_REC, "BSE_CBD", "BSE_CBD_REC", LOGFILE_CLASS::BINARY }}, { LOGFILE::BSE_DETAILED_OUTPUT, { "BSE_Detailed_Output", BSE_DETAILED_OUTPUT_REC, "BSE_DETAILED", "BSE_DETAILED_REC", LOGFILE_CLASS::BINARY }}, { LOGFILE::BSE_DOUBLE_COMPACT_OBJECTS, { "BSE_Double_Compact_Objects", BSE_DOUBLE_COMPACT_OBJECTS_REC, "BSE_DCO", "BSE_DCO_REC", LOGFILE_CLASS::BINARY }}, { LOGFILE::BSE_PULSAR_EVOLUTION, { "BSE_Pulsar_Evolution", BSE_PULSAR_EVOLUTION_REC, "BSE_PULSARS", "BSE_PULSARS_REC", LOGFILE_CLASS::BINARY }}, diff --git a/src/Makefile b/src/Makefile index c5b9ba5ed..cb5481f55 100644 --- a/src/Makefile +++ b/src/Makefile @@ -68,10 +68,23 @@ SOURCES := \ \ BinaryConstituentStar.cpp \ BaseBinaryStar.cpp \ + CircumbinaryDiskFormation.cpp \ BinaryStar.cpp \ \ main.cpp +# Optional spindler-c backend. Default OFF keeps ordinary COMPAS builds +# independent of the external spindler-c and librinterpolate libraries. +# To enable the backend, either set this here or provide it on the make command +# line, e.g. make SPINDLER_C_ROOT=/path/to/spindler-c +SPINDLER_C_ROOT ?= + +ifneq ($(strip $(SPINDLER_C_ROOT)),) +SOURCES += CircumbinaryDiskEvolutionSpindler.cpp +else +SOURCES += CircumbinaryDiskEvolutionUnavailable.cpp +endif + OBJS := $(patsubst %.cpp,$(ODIR)/%.o,$(SOURCES)) # GSL, Boost, and HDF5 directories @@ -86,6 +99,50 @@ BOOSTLIBDIR ?= /lib HDF5INCDIR ?= /usr/include/hdf5/serial HDF5LIBDIR ?= /usr/lib/$(shell $(CPP) -print-multiarch)/hdf5/serial +# spindler-c include/link flags are populated only when SPINDLER_C_ROOT is +# nonempty. +SPINDLER_C_CPPFLAGS := +SPINDLER_C_LDFLAGS := +SPINDLER_C_LIBS := + +ifneq ($(strip $(SPINDLER_C_ROOT)),) + +SPINDLER_C_INCDIR ?= $(SPINDLER_C_ROOT)/include +SPINDLER_C_LIBDIR ?= $(SPINDLER_C_ROOT)/lib + +# librinterpolate is a dependency of spindler-c, linked as -lrinterpolate. +# Leave RINTERPOLATE_ROOT empty if librinterpolate is installed in a standard +# compiler/linker search path. Otherwise, set it here or on the make command +# line, e.g. make SPINDLER_C_ROOT=/path/to/spindler-c RINTERPOLATE_ROOT=/path/to/rinterpolate +RINTERPOLATE_ROOT ?= +RINTERPOLATE_INCDIR ?= +RINTERPOLATE_LIBDIR ?= + +ifneq ($(strip $(RINTERPOLATE_ROOT)),) +ifeq ($(strip $(RINTERPOLATE_INCDIR)),) +RINTERPOLATE_INCDIR := $(RINTERPOLATE_ROOT)/include +endif +ifeq ($(strip $(RINTERPOLATE_LIBDIR)),) +RINTERPOLATE_LIBDIR := $(RINTERPOLATE_ROOT)/lib +endif +endif + +SPINDLER_C_CPPFLAGS += -DCOMPAS_ENABLE_SPINDLER_C +SPINDLER_C_CPPFLAGS += -I$(SPINDLER_C_INCDIR) +ifneq ($(strip $(RINTERPOLATE_INCDIR)),) +SPINDLER_C_CPPFLAGS += -I$(RINTERPOLATE_INCDIR) +endif +SPINDLER_C_CPPFLAGS += -DSPINDLER_C_ROOT=\"$(SPINDLER_C_ROOT)\" +SPINDLER_C_LDFLAGS += -L$(SPINDLER_C_LIBDIR) \ + -Wl,-rpath,$(SPINDLER_C_LIBDIR) +ifneq ($(strip $(RINTERPOLATE_LIBDIR)),) +SPINDLER_C_LDFLAGS += -L$(RINTERPOLATE_LIBDIR) \ + -Wl,-rpath,$(RINTERPOLATE_LIBDIR) +endif +SPINDLER_C_LIBS += -lspindler -lrinterpolate + +endif + # Homebrew defaults for macOS local builds ifeq ($(UNAME_S),Darwin) ifneq ($(wildcard $(HOMEBREW_PREFIX)/opt/boost/include/boost/version.hpp),) @@ -108,7 +165,7 @@ endif # might need -Wno-vla-cxx-extension for clang>18? EXTRA_WARN_SUPPRESS := -Wno-vla-extension -CXXFLAGS := -std=c++17 -Wall $(OPTFLAGS) $(CXXFLAGS_SPECIFIC) $(EXTRA_WARN_SUPPRESS) +CXXFLAGS := -std=c++17 -Wall $(OPTFLAGS) $(CXXFLAGS_SPECIFIC) $(EXTRA_WARN_SUPPRESS) $(SPINDLER_C_CPPFLAGS) ICFLAGS := -I$(GSLINCDIR) -I$(BOOSTINCDIR) -I$(HDF5INCDIR) -I. LIBS := -lm -lz -ldl -lpthread @@ -120,7 +177,7 @@ ifeq ($(wildcard $(BOOSTLIBDIR)/libboost_system*),) endif endif HDF5LIBS := -lhdf5_hl_cpp -lhdf5_cpp -lhdf5_hl -lhdf5 -LFLAGS := -L$(GSLLIBDIR) -L$(BOOSTLIBDIR) -L$(HDF5LIBDIR) -Xlinker -rpath -Xlinker $(BOOSTLIBDIR) $(HDF5LIBS) $(LIBS) $(GSLLIBS) $(BOOSTLIBS) $(LFLAGS_SPECIFIC) +LFLAGS := -L$(GSLLIBDIR) -L$(BOOSTLIBDIR) -L$(HDF5LIBDIR) $(SPINDLER_C_LDFLAGS) -Xlinker -rpath -Xlinker $(BOOSTLIBDIR) $(HDF5LIBS) $(LIBS) $(GSLLIBS) $(BOOSTLIBS) $(SPINDLER_C_LIBS) $(LFLAGS_SPECIFIC) .PHONY: all static fast staticfast clean link diff --git a/src/Options.cpp b/src/Options.cpp index ec62db46c..c61a5e72e 100644 --- a/src/Options.cpp +++ b/src/Options.cpp @@ -555,6 +555,21 @@ void Options::OptionValues::Initialise() { m_CommonEnvelopeMassAccretionConstant = 0.0; // Common envelope formalism + m_PostCommonEnvelopeEccentricityPrescription.type = POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION::CIRCULAR; + m_PostCommonEnvelopeEccentricityPrescription.typeString = POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION_LABEL.at(m_PostCommonEnvelopeEccentricityPrescription.type); + m_PostCommonEnvelopeEccentricityCap = 0.1; + m_PostCommonEnvelopeEccentricityFraction = 0.1; + m_CircumbinaryDisk = false; + m_CircumbinaryDiskBeta = 3.0; + m_CircumbinaryDiskEddingtonMode.type = CIRCUMBINARY_DISK_EDDINGTON_MODE::GE23; + m_CircumbinaryDiskEddingtonMode.typeString = CIRCUMBINARY_DISK_EDDINGTON_MODE_LABEL.at(m_CircumbinaryDiskEddingtonMode.type); + m_CircumbinaryDiskEddingtonCapFactor = 1.0; + m_CircumbinaryDiskEddingtonGe23Factor = 2.0 / 3.0; + m_CircumbinaryDiskInnerRadiusOverSeparation = 2.5; + m_CircumbinaryDiskLifetime = 1.0e5; + m_CircumbinaryDiskStructureFactor = 4.0; + m_CircumbinaryDiskEnvelopeMassFractionSupplied = 0.05; + m_CommonEnvelopeFormalism.type = CE_FORMALISM::ENERGY; m_CommonEnvelopeFormalism.typeString = CE_FORMALISM_LABEL.at(m_CommonEnvelopeFormalism.type); @@ -850,6 +865,11 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt po::value(&p_Options->m_CirculariseBinaryDuringMassTransfer)->default_value(p_Options->m_CirculariseBinaryDuringMassTransfer)->implicit_value(true), ("Circularise binary when it enters a Mass Transfer episode (default = " + std::string(p_Options->m_CirculariseBinaryDuringMassTransfer ? "TRUE" : "FALSE") + ")").c_str() ) + ( + "circumbinary-disk", + po::value(&p_Options->m_CircumbinaryDisk)->default_value(p_Options->m_CircumbinaryDisk)->implicit_value(true), + ("Apply post-common-envelope circumbinary-disk evolution where the CBD formation criterion is satisfied (default = " + std::string(p_Options->m_CircumbinaryDisk ? "TRUE" : "FALSE") + ")").c_str() + ) ( "common-envelope-allow-immediate-RLOF-post-CE-survive", po::value(&p_Options->m_AllowImmediateRLOFpostCEToSurviveCommonEnvelope)->default_value(p_Options->m_AllowImmediateRLOFpostCEToSurviveCommonEnvelope)->implicit_value(true), @@ -1137,6 +1157,62 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt // double + // Circumbinary-disk and post-CE eccentricity options. + ( + "circumbinary-disk-beta", + po::value(&p_Options->m_CircumbinaryDiskBeta)->default_value(p_Options->m_CircumbinaryDiskBeta), + ("Circumbinary disk angular-momentum retention parameter beta = fJ/fM; must be non-negative, and beta * fM must be <= 1 (default = " + std::to_string(p_Options->m_CircumbinaryDiskBeta) + ")").c_str() + ) + ( + "circumbinary-disk-eddington-cap-factor", + po::value(&p_Options->m_CircumbinaryDiskEddingtonCapFactor)->default_value(p_Options->m_CircumbinaryDiskEddingtonCapFactor), + ("Multiplier applied to the Eddington accretion rate in CBD Eddington mode CAP (default = " + std::to_string(p_Options->m_CircumbinaryDiskEddingtonCapFactor) + ")").c_str() + ) + ( + "circumbinary-disk-eddington-ge23-factor", + po::value(&p_Options->m_CircumbinaryDiskEddingtonGe23Factor)->default_value(p_Options->m_CircumbinaryDiskEddingtonGe23Factor), + ("Retained fraction of the requested super-Eddington CBD accretion rate in Eddington mode GE23; must satisfy 0 < f <= 1 (default = " + std::to_string(p_Options->m_CircumbinaryDiskEddingtonGe23Factor) + ")").c_str() + ) + ( + "circumbinary-disk-inner-radius-over-separation", + po::value(&p_Options->m_CircumbinaryDiskInnerRadiusOverSeparation)->default_value(p_Options->m_CircumbinaryDiskInnerRadiusOverSeparation), + ("Inner circumbinary-disk radius in units of post-CE binary separation; must be >= 1 (default = " + std::to_string(p_Options->m_CircumbinaryDiskInnerRadiusOverSeparation) + ")").c_str() + ) + ( + "circumbinary-disk-lifetime", + po::value(&p_Options->m_CircumbinaryDiskLifetime)->default_value(p_Options->m_CircumbinaryDiskLifetime), + ("Circumbinary-disk lifetime in years; must be > 0 (default = " + std::to_string(p_Options->m_CircumbinaryDiskLifetime) + ")").c_str() + ) + ( + "circumbinary-disk-structure-factor", + po::value(&p_Options->m_CircumbinaryDiskStructureFactor)->default_value(p_Options->m_CircumbinaryDiskStructureFactor), + ("Disk-structure factor multiplying the CBD angular-momentum threshold; must be > 1 (default = " + std::to_string(p_Options->m_CircumbinaryDiskStructureFactor) + ")").c_str() + ) + ( + "circumbinary-disk-envelope-mass-fraction-supplied", + po::value(&p_Options->m_CircumbinaryDiskEnvelopeMassFractionSupplied)->default_value(p_Options->m_CircumbinaryDiskEnvelopeMassFractionSupplied), + ("Fraction of the total binary mass lost during common-envelope evolution that is supplied through a circumbinary disk; the resulting initial CBD mass (logged as Initial_CBD_Mass) is assumed to be much smaller than the post-CE binary mass in the Spindler formalism; must be between 0 and 0.2 (default = " + std::to_string(p_Options->m_CircumbinaryDiskEnvelopeMassFractionSupplied) + ")").c_str() + ) + ( + "post-common-envelope-eccentricity-prescription", + po::value(&p_Options->m_PostCommonEnvelopeEccentricityPrescription.typeString)->default_value(p_Options->m_PostCommonEnvelopeEccentricityPrescription.typeString), + ("Post-common-envelope eccentricity prescription: CIRCULAR, PRE_RLOF_CAP, or PRE_RLOF_FRACTION (default = " + p_Options->m_PostCommonEnvelopeEccentricityPrescription.typeString + ")").c_str() + ) + ( + "circumbinary-disk-eddington-mode", + po::value(&p_Options->m_CircumbinaryDiskEddingtonMode.typeString)->default_value(p_Options->m_CircumbinaryDiskEddingtonMode.typeString), + ("CBD Eddington-limited accretion mode: CAP, TM23, or GE23 (default = " + p_Options->m_CircumbinaryDiskEddingtonMode.typeString + ")").c_str() + ) + ( + "post-common-envelope-eccentricity-cap", + po::value(&p_Options->m_PostCommonEnvelopeEccentricityCap)->default_value(p_Options->m_PostCommonEnvelopeEccentricityCap), + ("Upper cap on the post-common-envelope eccentricity when --post-common-envelope-eccentricity-prescription is PRE_RLOF_CAP; must be in [0, 1) (default = " + std::to_string(p_Options->m_PostCommonEnvelopeEccentricityCap) + ")").c_str() + ) + ( + "post-common-envelope-eccentricity-fraction", + po::value(&p_Options->m_PostCommonEnvelopeEccentricityFraction)->default_value(p_Options->m_PostCommonEnvelopeEccentricityFraction), + ("Fraction of the pre-RLOF eccentricity retained after common envelope evolution when --post-common-envelope-eccentricity-prescription is PRE_RLOF_FRACTION; must be in [0, 1] (default = " + std::to_string(p_Options->m_PostCommonEnvelopeEccentricityFraction) + ")").c_str() + ) ( "common-envelope-alpha", po::value(&p_Options->m_CommonEnvelopeAlpha)->default_value(p_Options->m_CommonEnvelopeAlpha), @@ -2376,6 +2452,16 @@ std::string Options::OptionValues::CheckAndSetOptions() { COMPLAIN_IF(!found, "Unknown Mass Transfer Angular Momentum Loss Prescription for 2-stage CE"); } + if (!DEFAULTED("circumbinary-disk-eddington-mode")) { + std::tie(found, m_CircumbinaryDiskEddingtonMode.type) = utils::GetMapKey(m_CircumbinaryDiskEddingtonMode.typeString, CIRCUMBINARY_DISK_EDDINGTON_MODE_LABEL, m_CircumbinaryDiskEddingtonMode.type); + COMPLAIN_IF(!found, "Unknown circumbinary disk Eddington mode: use CAP, GE23, or TM23"); + } + + if (!DEFAULTED("post-common-envelope-eccentricity-prescription")) { + std::tie(found, m_PostCommonEnvelopeEccentricityPrescription.type) = utils::GetMapKey(m_PostCommonEnvelopeEccentricityPrescription.typeString, POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION_LABEL, m_PostCommonEnvelopeEccentricityPrescription.type); + COMPLAIN_IF(!found, "Unknown post-CE eccentricity prescription (--post-common-envelope-eccentricity-prescription). Allowed: CIRCULAR, PRE_RLOF_CAP, PRE_RLOF_FRACTION"); + } + if (!DEFAULTED("critical-mass-ratio-prescription")) { // critical mass ratio prescription std::tie(found, m_QCritPrescription.type) = utils::GetMapKey(m_QCritPrescription.typeString, QCRIT_PRESCRIPTION_LABEL, m_QCritPrescription.type); COMPLAIN_IF(!found, "Unknown qCrit Prescription"); @@ -2566,6 +2652,32 @@ std::string Options::OptionValues::CheckAndSetOptions() { // constraint/value/range checks - alphabetically (where possible) + COMPLAIN_IF(utils::Compare(m_PostCommonEnvelopeEccentricityCap, 0.0) < 0 || utils::Compare(m_PostCommonEnvelopeEccentricityCap, 1.0) >= 0, + "Post-CE eccentricity cap (--post-common-envelope-eccentricity-cap) must be in [0, 1)"); + COMPLAIN_IF(utils::Compare(m_PostCommonEnvelopeEccentricityFraction, 0.0) < 0 || utils::Compare(m_PostCommonEnvelopeEccentricityFraction, 1.0) > 0, + "Post-CE eccentricity fraction (--post-common-envelope-eccentricity-fraction) must be in [0, 1]"); + + COMPLAIN_IF(utils::Compare(m_CircumbinaryDiskBeta, 0.0) < 0, + "CBD beta (--circumbinary-disk-beta) must be >= 0"); + COMPLAIN_IF(utils::Compare(m_CircumbinaryDiskEddingtonCapFactor, 0.0) <= 0, + "CBD Eddington CAP factor (--circumbinary-disk-eddington-cap-factor) must be > 0"); + COMPLAIN_IF(utils::Compare(m_CircumbinaryDiskEddingtonGe23Factor, 0.0) <= 0 || utils::Compare(m_CircumbinaryDiskEddingtonGe23Factor, 1.0) > 0, + "CBD Eddington GE23 factor (--circumbinary-disk-eddington-ge23-factor) must satisfy 0 < f <= 1"); + COMPLAIN_IF(utils::Compare(m_CircumbinaryDiskEnvelopeMassFractionSupplied, 0.0) < 0 || utils::Compare(m_CircumbinaryDiskEnvelopeMassFractionSupplied, CBD_MAXIMUM_ENVELOPE_MASS_FRACTION_SUPPLIED) > 0, + "CBD CE mass-loss fraction supplied (--circumbinary-disk-envelope-mass-fraction-supplied) must be between 0 and 0.2"); + COMPLAIN_IF(utils::Compare(m_CircumbinaryDiskBeta * m_CircumbinaryDiskEnvelopeMassFractionSupplied, 1.0) > 0, + "CBD beta * CE mass-loss fraction supplied must be <= 1"); + COMPLAIN_IF(utils::Compare(m_CircumbinaryDiskInnerRadiusOverSeparation, 1.0) < 0, + "CBD inner radius over separation (--circumbinary-disk-inner-radius-over-separation) must be >= 1"); + COMPLAIN_IF(utils::Compare(m_CircumbinaryDiskLifetime, 0.0) <= 0, + "CBD lifetime (--circumbinary-disk-lifetime) must be > 0 yr"); + COMPLAIN_IF(utils::Compare(m_CircumbinaryDiskStructureFactor, 1.0) <= 0, + "CBD structure factor (--circumbinary-disk-structure-factor) must be > 1"); +#ifndef COMPAS_ENABLE_SPINDLER_C + COMPLAIN_IF(m_CircumbinaryDisk, + "Circumbinary-disk evolution was requested, but this COMPAS executable was built without spindler-c support. Rebuild with SPINDLER_C_ROOT set to the spindler-c build tree."); +#endif + COMPLAIN_IF(m_CommonEnvelopeAlpha < 0.0, "CE alpha (--common-envelope-alpha) < 0"); COMPLAIN_IF(m_CommonEnvelopeAlphaThermal < 0.0 || m_CommonEnvelopeAlphaThermal > 1.0, "CE alpha thermal (--common-envelope-alpha-thermal) must be between 0 and 1"); COMPLAIN_IF(m_CommonEnvelopeLambdaMultiplier < 0.0, "CE lambda multiplie (--common-envelope-lambda-multiplier < 0"); @@ -2819,6 +2931,8 @@ STR_VECTOR Options::AllowedOptionValues(const std::string p_OptionString) { case _("chemically-homogeneous-evolution-mode") : POPULATE_RET(CHE_MODE_LABEL); break; case _("common-envelope-formalism") : POPULATE_RET(CE_FORMALISM_LABEL); break; case _("common-envelope-lambda-prescription") : POPULATE_RET(CE_LAMBDA_PRESCRIPTION_LABEL); break; + case _("circumbinary-disk-eddington-mode") : POPULATE_RET(CIRCUMBINARY_DISK_EDDINGTON_MODE_LABEL); break; + case _("post-common-envelope-eccentricity-prescription") : POPULATE_RET(POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION_LABEL); break; case _("common-envelope-mass-accretion-prescription") : POPULATE_RET(CE_ACCRETION_PRESCRIPTION_LABEL); break; case _("common-envelope-second-stage-gamma-prescription") : POPULATE_RET(MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION_LABEL); break; case _("critical-mass-ratio-prescription") : POPULATE_RET(QCRIT_PRESCRIPTION_LABEL); break; @@ -4929,6 +5043,19 @@ COMPAS_VARIABLE Options::OptionValue(const T_ANY_PROPERTY p_Property) const { case PROGRAM_OPTION::CIRCULARISE_BINARY_DURING_MT : value = CirculariseBinaryDuringMassTransfer(); break; + case PROGRAM_OPTION::CIRCUMBINARY_DISK : value = CircumbinaryDisk(); break; + case PROGRAM_OPTION::CIRCUMBINARY_DISK_BETA : value = CircumbinaryDiskBeta(); break; + case PROGRAM_OPTION::CIRCUMBINARY_DISK_EDDINGTON_CAP_FACTOR : value = CircumbinaryDiskEddingtonCapFactor(); break; + case PROGRAM_OPTION::CIRCUMBINARY_DISK_EDDINGTON_GE23_FACTOR : value = CircumbinaryDiskEddingtonGe23Factor(); break; + case PROGRAM_OPTION::CIRCUMBINARY_DISK_EDDINGTON_MODE : value = static_cast(CircumbinaryDiskEddingtonMode()); break; + case PROGRAM_OPTION::CIRCUMBINARY_DISK_INNER_RADIUS_OVER_SEPARATION : value = CircumbinaryDiskInnerRadiusOverSeparation(); break; + case PROGRAM_OPTION::CIRCUMBINARY_DISK_LIFETIME : value = CircumbinaryDiskLifetime(); break; + case PROGRAM_OPTION::CIRCUMBINARY_DISK_STRUCTURE_FACTOR : value = CircumbinaryDiskStructureFactor(); break; + case PROGRAM_OPTION::CIRCUMBINARY_DISK_ENVELOPE_MASS_FRACTION_SUPPLIED : value = CircumbinaryDiskEnvelopeMassFractionSupplied(); break; + case PROGRAM_OPTION::POST_COMMON_ENVELOPE_ECCENTRICITY_CAP : value = PostCommonEnvelopeEccentricityCap(); break; + case PROGRAM_OPTION::POST_COMMON_ENVELOPE_ECCENTRICITY_FRACTION : value = PostCommonEnvelopeEccentricityFraction(); break; + case PROGRAM_OPTION::POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION : value = static_cast(PostCommonEnvelopeEccentricityPrescription()); break; + case PROGRAM_OPTION::COMMON_ENVELOPE_ALPHA : value = CommonEnvelopeAlpha(); break; case PROGRAM_OPTION::COMMON_ENVELOPE_ALPHA_THERMAL : value = CommonEnvelopeAlphaThermal(); break; case PROGRAM_OPTION::COMMON_ENVELOPE_FORMALISM : value = static_cast(CommonEnvelopeFormalism()); break; diff --git a/src/Options.h b/src/Options.h index a9c1eb162..0e288a298 100755 --- a/src/Options.h +++ b/src/Options.h @@ -19,7 +19,6 @@ #include #include "constants.h" - #include // Boost string manipulation #include // Boost command line options tools #include // Boost filesystem tools for handling paths etc. @@ -469,6 +468,15 @@ class Options { "case-BB-stability-prescription", "circularise-binary-during-mass-transfer", + "circumbinary-disk", + "circumbinary-disk-beta", + "circumbinary-disk-eddington-cap-factor", + "circumbinary-disk-eddington-ge23-factor", + "circumbinary-disk-eddington-mode", + "circumbinary-disk-envelope-mass-fraction-supplied", + "circumbinary-disk-inner-radius-over-separation", + "circumbinary-disk-lifetime", + "circumbinary-disk-structure-factor", "common-envelope-allow-main-sequence-survive", "common-envelope-alpha", "common-envelope-alpha-thermal", @@ -544,6 +552,10 @@ class Options { "orbital-period-max", "orbital-period-min", + "post-common-envelope-eccentricity-cap", + "post-common-envelope-eccentricity-fraction", + "post-common-envelope-eccentricity-prescription", + "rlof-printing", "rotational-frequency-1", "rotational-frequency-2", @@ -589,6 +601,15 @@ class Options { "check-photon-tiring-limit", "chemically-homogeneous-evolution-mode", "circularise-binary-during-mass-transfer", + "circumbinary-disk", + "circumbinary-disk-beta", + "circumbinary-disk-eddington-cap-factor", + "circumbinary-disk-eddington-ge23-factor", + "circumbinary-disk-eddington-mode", + "circumbinary-disk-envelope-mass-fraction-supplied", + "circumbinary-disk-inner-radius-over-separation", + "circumbinary-disk-lifetime", + "circumbinary-disk-structure-factor", "common-envelope-allow-main-sequence-survive", "common-envelope-formalism", "common-envelope-lambda-prescription", @@ -1102,6 +1123,20 @@ class Options { double m_MassTransferCriticalMassRatioWhiteDwarfNonDegenerateAccretor; // Critical mass ratio for MT from a white dwarf double m_MassTransferCriticalMassRatioWhiteDwarfDegenerateAccretor; // Critical mass ratio for MT from a white dwarf on to a degenerate accretor + // Circumbinary disk and post-CE eccentricity options + ENUM_OPT m_PostCommonEnvelopeEccentricityPrescription; // Prescription for setting eccentricity after common envelope evolution + double m_PostCommonEnvelopeEccentricityCap; // Eccentricity cap assigned after common envelope evolution when requested + double m_PostCommonEnvelopeEccentricityFraction; // Fraction of the pre-RLOF eccentricity retained after common envelope evolution when requested + bool m_CircumbinaryDisk; // Whether to apply post-CE CBD evolution where the formation criterion is satisfied + double m_CircumbinaryDiskBeta; // CBD angular-momentum parameter beta = fJ/fM + ENUM_OPT m_CircumbinaryDiskEddingtonMode; // CBD Eddington-limited accretion mode: CAP, GE23, or TM23 + double m_CircumbinaryDiskEddingtonCapFactor; // CAP-mode multiplier applied to the Eddington rate + double m_CircumbinaryDiskEddingtonGe23Factor; // GE23-mode retained fraction for super-Eddington requested accretion + double m_CircumbinaryDiskInnerRadiusOverSeparation; // Inner CBD radius in units of post-CE binary separation + double m_CircumbinaryDiskLifetime; // CBD lifetime in yr + double m_CircumbinaryDiskStructureFactor; // Disk-structure factor multiplying the CBD angular-momentum threshold + double m_CircumbinaryDiskEnvelopeMassFractionSupplied; // Fraction of CE mass loss supplied through a CBD; assumed to yield initial m_CBD << post-CE binary mass + // Common Envelope options double m_CommonEnvelopeAlpha; // Common envelope efficiency alpha parameter double m_CommonEnvelopeLambda; // Common envelope Lambda parameter @@ -1432,6 +1467,18 @@ class Options { bool CommandLineGrid() const { return m_CmdLine.complexOptionValues.size() != 0; } + bool CircumbinaryDisk() const { return OPT_VALUE("circumbinary-disk", m_CircumbinaryDisk, true); } + double CircumbinaryDiskBeta() const { return OPT_VALUE("circumbinary-disk-beta", m_CircumbinaryDiskBeta, true); } + CIRCUMBINARY_DISK_EDDINGTON_MODE CircumbinaryDiskEddingtonMode() const { return OPT_VALUE("circumbinary-disk-eddington-mode", m_CircumbinaryDiskEddingtonMode.type, true); } + double CircumbinaryDiskEddingtonCapFactor() const { return OPT_VALUE("circumbinary-disk-eddington-cap-factor", m_CircumbinaryDiskEddingtonCapFactor, true); } + double CircumbinaryDiskEddingtonGe23Factor() const { return OPT_VALUE("circumbinary-disk-eddington-ge23-factor", m_CircumbinaryDiskEddingtonGe23Factor, true); } + double CircumbinaryDiskInnerRadiusOverSeparation() const { return OPT_VALUE("circumbinary-disk-inner-radius-over-separation", m_CircumbinaryDiskInnerRadiusOverSeparation, true); } + double CircumbinaryDiskLifetime() const { return OPT_VALUE("circumbinary-disk-lifetime", m_CircumbinaryDiskLifetime, true); } + double CircumbinaryDiskStructureFactor() const { return OPT_VALUE("circumbinary-disk-structure-factor", m_CircumbinaryDiskStructureFactor, true); } + double CircumbinaryDiskEnvelopeMassFractionSupplied() const { return OPT_VALUE("circumbinary-disk-envelope-mass-fraction-supplied", m_CircumbinaryDiskEnvelopeMassFractionSupplied, true); } + double PostCommonEnvelopeEccentricityCap() const { return OPT_VALUE("post-common-envelope-eccentricity-cap", m_PostCommonEnvelopeEccentricityCap, true); } + double PostCommonEnvelopeEccentricityFraction() const { return OPT_VALUE("post-common-envelope-eccentricity-fraction", m_PostCommonEnvelopeEccentricityFraction, true); } + POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION PostCommonEnvelopeEccentricityPrescription() const { return OPT_VALUE("post-common-envelope-eccentricity-prescription", m_PostCommonEnvelopeEccentricityPrescription.type, true); } double CommonEnvelopeAlpha() const { return OPT_VALUE("common-envelope-alpha", m_CommonEnvelopeAlpha, true); } double CommonEnvelopeAlphaThermal() const { return OPT_VALUE("common-envelope-alpha-thermal", m_CommonEnvelopeAlphaThermal, true); } CE_FORMALISM CommonEnvelopeFormalism() const { return OPT_VALUE("common-envelope-formalism", m_CommonEnvelopeFormalism.type, true); } diff --git a/src/changelog.h b/src/changelog.h index 1e22d11e0..975820a21 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1712,7 +1712,25 @@ // - Add option (MIN) to lose specific angular momentum during non-conservative mass transfer that is the smaller of isotropic re-emission and Jeans mode mass loss // - Defect repair: // - Corrected orbital separation re-computation following pulsational (L/M threshold) mass loss -// +// 03.31.00 AG - November 01, 2026 - Enhancement: +// Added a post-common-envelope circumbinary disk (CBD) framework ā€“ disabled by default - based on the formalism +// described in Grichener et al. (2026). The framework assumes that after common-envelope +// evolution (CEE), a fraction of the pre-RLOF envelope mass and the binary's total angular +// momentum (AM) might remain available for CBD formation. If the available AM suffices to support a CBD of the retained mass, one is formed. +// The CBD is then accreted back by the binary, with the evolution based on spindler-c (*), modified to account for the Eddington accretion limit. +// Non-accreted material is assumed to be lost in winds carrying the accretor's specific orbital angular momentum. +// CBD related output is written in a new file called BSE_Circumbinary_Disks. +// Main user options: +// --circumbinary-disk enables the CBD framework +// --circumbinary-disk-beta controls the likelihood of CBD formation +// --circumbinary-disk-envelope-mass-fraction-supplied controls the CBD mass +// --circumbinary-disk-eddington-mode controls the Eddington-limited accretion prescription +// --post-common-envelope-eccentricity-prescription controls the post-CEE eccentricity (non-zero values tend to be significantly increased by the CBD-binary interaction) +// See documentation for more information on the available options. +// (*) https://github.com/ruggero-valli/spindler-c, installation required to use the CBD module + + + // Version string format is MM.mm.rr, where // // MM is the MAJOR release number: this should be incremented whenever major new functionality is introduced @@ -1722,7 +1740,7 @@ // if MM is incremented, set mm and rr to 00, even if defect repairs and minor enhancements were also made // if mm is incremented, set rr to 00, even if defect repairs were also made -const std::string VERSION_STRING = "03.30.00"; +const std::string VERSION_STRING = "03.31.00"; # endif // __changelog_h__ diff --git a/src/constants.h b/src/constants.h index ddfb8ec3e..ddb789d34 100755 --- a/src/constants.h +++ b/src/constants.h @@ -255,6 +255,26 @@ constexpr unsigned int ABSOLUTE_MAXIMUM_TIMESTEPS = 1000000; constexpr int MAX_BSE_INITIAL_CONDITIONS_ITERATIONS = 100; // Maximum loop iterations looking for initial conditions for binary systems constexpr int MAX_TIMESTEP_RETRIES = 30; // Maximum retries to find a good timestep for stellar evolution +// Circumbinary disk (CBD) model limits and physical prescription coefficients. +// General timestep and option bounds are kept here, as for tides and mass loss. +constexpr double CBD_MAXIMUM_SUPPLIED_MASS_FRACTION_PER_TIMESTEP = 1.0E-3; // Maximum CBD mass supplied per timestep as a fraction of current binary mass +constexpr double CBD_MAXIMUM_ENVELOPE_MASS_FRACTION_SUPPLIED = 0.2; // Maximum fraction of CE mass lost allowed for the CBD reservoir + +// Siwek23 lookup-table domain (values outside the domain are clamped). +constexpr double CBD_SIWEK23_Q_MIN = 0.1; +constexpr double CBD_SIWEK23_Q_MAX = 1.0; +constexpr double CBD_SIWEK23_E_MIN = 0.0; +constexpr double CBD_SIWEK23_E_MAX = 0.8; + +// CBD-specific Eddington normalizations and TM23 retention prescription. +constexpr double CBD_BH_EDDINGTON_RATE_AT_REFERENCE_MASS = 4.42E-7; +constexpr double CBD_BH_EDDINGTON_REFERENCE_MASS_MSOL = 10.0; +constexpr double CBD_NS_EDDINGTON_REFERENCE_RADIUS_RSOL = 1.4374E-5; +constexpr double CBD_NS_EDDINGTON_RATE_AT_REFERENCE_RADIUS = 3.0E-8; +constexpr double CBD_WD_EDDINGTON_RATE_PER_RADIUS = 2.08E-3 / 1.7; +constexpr double CBD_TM23_RETENTION_COEFFICIENT = 0.6; +constexpr double CBD_TM23_RETENTION_EXPONENT = 0.6; + constexpr double MAXIMUM_MASS_LOSS_FRACTION = 0.001; // Maximum allowable mass loss - 0.1% (of mass) expressed as a fraction constexpr double MAXIMUM_RADIAL_CHANGE = 0.1; // Maximum allowable radial change - 10% (of radius) expressed as a fraction constexpr double MAXIMUM_WIND_MASS_LOSS_RATE = 0.1; // Maximum wind mass loss rates (in solar masses per year) to avoid convergence issues diff --git a/src/typedefs.h b/src/typedefs.h index d8a2ddc6d..278d9b886 100755 --- a/src/typedefs.h +++ b/src/typedefs.h @@ -352,6 +352,18 @@ const COMPASUnorderedMap CASE_BB_ST { CASE_BB_STABILITY_PRESCRIPTION::ALWAYS_UNSTABLE, "ALWAYS_UNSTABLE" } }; + +// circumbinary disk Eddington-limited accretion modes +// CAP applies a user-scaled Eddington cap; GE23 applies a fixed retained +// fraction motivated by Ghodla & Eldridge; TM23 applies the Tuna & Metzger +// compact-object prescription. +enum class CIRCUMBINARY_DISK_EDDINGTON_MODE: int { CAP, GE23, TM23 }; +const COMPASUnorderedMap CIRCUMBINARY_DISK_EDDINGTON_MODE_LABEL = { + { CIRCUMBINARY_DISK_EDDINGTON_MODE::CAP, "CAP" }, + { CIRCUMBINARY_DISK_EDDINGTON_MODE::GE23, "GE23" }, + { CIRCUMBINARY_DISK_EDDINGTON_MODE::TM23, "TM23" } +}; + // common envelope cccretion prescriptions enum class CE_ACCRETION_PRESCRIPTION: int { ZERO, CONSTANT, UNIFORM, MACLEOD, CHEVALIER }; const COMPASUnorderedMap CE_ACCRETION_PRESCRIPTION_LABEL = { @@ -762,6 +774,18 @@ const COMPASUnorderedMap OBJECT_TYPE_LABEL = { { OBJECT_TYPE::BINARY_CONSTITUENT_STAR, "BinaryConstituentStar" } }; + +// post-common-envelope eccentricity prescriptions +// CIRCULAR uses the usual e=0 CE outcome; PRE_RLOF_CAP uses +// min(e_pre_RLOF, eccentricity_cap); PRE_RLOF_FRACTION uses +// eccentricity_fraction * e_pre_RLOF. +enum class POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION: int { CIRCULAR, PRE_RLOF_CAP, PRE_RLOF_FRACTION }; +const COMPASUnorderedMap POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION_LABEL = { + { POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION::CIRCULAR, "CIRCULAR" }, + { POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION::PRE_RLOF_CAP, "PRE_RLOF_CAP" }, + { POST_COMMON_ENVELOPE_ECCENTRICITY_PRESCRIPTION::PRE_RLOF_FRACTION, "PRE_RLOF_FRACTION" } +}; + // program options origin indicator (command line or gridfile line) enum class OPTIONS_ORIGIN: int { CMDLINE, GRIDFILE }; diff --git a/src/yaml.h b/src/yaml.h index e64190982..2ffa15da3 100644 --- a/src/yaml.h +++ b/src/yaml.h @@ -119,6 +119,9 @@ namespace yaml { " --common-envelope-lambda-nanjing-use-rejuvenated-mass", " --revised-energy-formalism-nandez-ivanova", "", + " ### CIRCUMBINARY DISK", + " --circumbinary-disk", + "", " ### SUPERNOVAE, KICKS AND REMNANTS", " --allow-non-stripped-ECSN", " --pair-instability-supernovae", @@ -235,6 +238,17 @@ namespace yaml { " --common-envelope-second-stage-beta", " --common-envelope-slope-kruckow", " --maximum-mass-donor-nandez-ivanova", + " --post-common-envelope-eccentricity-cap", + " --post-common-envelope-eccentricity-fraction", + "", + " ### CIRCUMBINARY DISK", + " --circumbinary-disk-beta", + " --circumbinary-disk-eddington-cap-factor", + " --circumbinary-disk-eddington-ge23-factor", + " --circumbinary-disk-envelope-mass-fraction-supplied", + " --circumbinary-disk-inner-radius-over-separation", + " --circumbinary-disk-lifetime", + " --circumbinary-disk-structure-factor", "", " ### SUPERNOVAE, KICKS AND REMNANTS", " --eddington-accretion-factor # multiplication Factor for eddington accretion onto NS&BH", @@ -341,6 +355,10 @@ namespace yaml { " --common-envelope-lambda-prescription # Xu & Li 2010", " --common-envelope-mass-accretion-prescription", " --common-envelope-second-stage-gamma-prescription", + " --post-common-envelope-eccentricity-prescription", + "", + " ### CIRCUMBINARY DISK", + " --circumbinary-disk-eddington-mode", "", " ### TIDES", " --tides-prescription",