diff --git a/.github/workflows/compas-compile-ci.yml b/.github/workflows/compas-compile-ci.yml index 2459bbe91..e860330e4 100644 --- a/.github/workflows/compas-compile-ci.yml +++ b/.github/workflows/compas-compile-ci.yml @@ -12,6 +12,12 @@ name: COMPAS compile test push: branches: - dev + - +# Ensures only the latest workflow run for the same branch is active, canceling any in-progress runs. +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + jobs: compas: env: @@ -25,17 +31,18 @@ jobs: fail-fast: false matrix: os: - - ubuntu-20.04 + - ubuntu-22.04 steps: - - uses: actions/checkout@v3 - - uses: actions/setup-python@v4 + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 with: python-version: '3.9' cache: pip + cache-dependency-path: setup.py - name: Install TeXLive uses: teatimeguest/setup-texlive-action@v3 - name: Install dependencies on ubuntu - if: 'startsWith(matrix.os, ''ubuntu-20'')' + if: 'startsWith(matrix.os, ''ubuntu-2'')' run: | cd misc/cicd-scripts chmod 755 linux-dependencies @@ -66,7 +73,7 @@ jobs: coverage html coverage-badge -o coverage_badge.svg -f - name: Archive code coverage results - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: code-coverage path: | @@ -74,8 +81,15 @@ jobs: coverage_badge.svg - name: Archive COMPAS detailed-evolution plot id: upload - uses: actions/upload-artifact@v3.1.2 + uses: actions/upload-artifact@v4 with: name: '${{ env.ARTIFACT_NAME }}' path: '${{ env.ARTIFACT_PATH }}/${{ env.ARTIFACT_NAME }}' if-no-files-found: error + - name: Test Summary + run: | + echo "### Test Results" >> $GITHUB_STEP_SUMMARY + echo "- Compas Build: Success" >> $GITHUB_STEP_SUMMARY + echo "- Python Utils Installation: Success" >> $GITHUB_STEP_SUMMARY + echo "- Example COMPAS Job: Success" >> $GITHUB_STEP_SUMMARY + echo "- Pytest Execution: Success" >> $GITHUB_STEP_SUMMARY \ No newline at end of file diff --git a/.github/workflows/dockerhub-ci.yml b/.github/workflows/dockerhub-ci.yml index 2fadc4412..f032aa9a7 100755 --- a/.github/workflows/dockerhub-ci.yml +++ b/.github/workflows/dockerhub-ci.yml @@ -5,26 +5,41 @@ on: branches: [ dev ] workflow_dispatch: null +# Ensures only the latest workflow run for the same branch is active, canceling any in-progress runs. +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + jobs: release: name: Build Docker image and push to DockerHub - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 steps: - - uses: actions/checkout@v2 - + - uses: actions/checkout@v4 + - name: Log in to Dockerhub - run: docker login -u ${{ secrets.DOCKER_USERNAME }} -p ${{ secrets.DOCKER_PASSWORD }} - + uses: docker/login-action@v3 + with: + username: ${{ secrets.DOCKER_USERNAME }} + password: ${{ secrets.DOCKER_PASSWORD }} + - name: Get release version run: echo "COMPAS_VERSION=$(sed -n '/const std::string VERSION_STRING/,/^$/p' ./src/changelog.h | sed 's/.*"\(.*\)"[^"]*$/\1/')" >> $GITHUB_ENV - + - name: Print version run: echo $COMPAS_VERSION - - - name: Build and tag Docker image - run: docker build -t teamcompas/compas:$COMPAS_VERSION -t teamcompas/compas:latest . - - - name: Push Docker image - run: | - docker push teamcompas/compas:$COMPAS_VERSION - docker push teamcompas/compas:latest + + - name: Set up Docker Buildx + uses: docker/setup-buildx-action@v3 + + - name: Build and push Docker image + uses: docker/build-push-action@v5 + with: + context: . + push: true + tags: | + teamcompas/compas:${{ env.COMPAS_VERSION }} + teamcompas/compas:latest + cache-from: type=gha + cache-to: type=gha,mode=max + diff --git a/.github/workflows/latex-compile-ci.yml b/.github/workflows/latex-compile-ci.yml deleted file mode 100644 index 46d61c65b..000000000 --- a/.github/workflows/latex-compile-ci.yml +++ /dev/null @@ -1,71 +0,0 @@ -name: LaTeX compile test - -# Test Latex compiles on pull_requests -# Public documentation on push to dev - -on: - pull_request: - paths: ['main/**'] # only changes to src/ trigger this workflow - push: - branches: [dev] - paths: ['main/**'] - - -jobs: - latex: - name: Build LaTeX documentation - runs-on: ubuntu-18.04 - - steps: - - uses: actions/checkout@v2 - - - name: Install dependencies - run: sudo apt install texlive-base - - # Need to do pdflatex, bibtex, pdflatex, pdflatex - - name: pdflatex main - uses: dante-ev/latex-action@latest - with: - working_directory: docs/COMPAS_LaTeX - root_file: main.tex - compiler: pdflatex - args: -interaction=nonstopmode -shell-escape - - - name: bibtex main - uses: dante-ev/latex-action@latest - with: - working_directory: docs/COMPAS_LaTeX - root_file: main.aux - compiler: bibtex - args: - - - name: pdflatex main - uses: dante-ev/latex-action@latest - with: - working_directory: docs/COMPAS_LaTeX - root_file: main.tex - compiler: pdflatex - args: -interaction=nonstopmode -shell-escape - - - name: pdflatex main - uses: dante-ev/latex-action@latest - with: - working_directory: docs/COMPAS_LaTeX - root_file: main.tex - compiler: pdflatex - args: -interaction=nonstopmode -shell-escape - - - # If push to dev, publish the documentation - - name: publish documentation - if: github.event_name == 'push' - run: | - git checkout --orphan Documentation - PDF_FILE=COMPAS_Documentation.pdf - mv docs/COMPAS_LaTeX/main.pdf $PDF_FILE - git rm -rf . - git add -f $PDF_FILE - git config user.name "Team COMPAS" - git config user.email "<>" - git commit -m "Documentation update" - git push -f --set-upstream origin Documentation diff --git a/.github/workflows/pr_artifact_url_commenter.yml b/.github/workflows/pr_artifact_url_commenter.yml index dae1097fd..62bca978d 100644 --- a/.github/workflows/pr_artifact_url_commenter.yml +++ b/.github/workflows/pr_artifact_url_commenter.yml @@ -9,68 +9,67 @@ on: jobs: comment: runs-on: ubuntu-latest - continue-on-error: true + if: ${{ github.event.workflow_run.conclusion == 'success' }} steps: - name: Get Artifact and Pull request info env: - GITHUB_TOKEN: ${{ github.token }} + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} WORKFLOW_RUN_EVENT_OBJ: ${{ toJSON(github.event.workflow_run) }} OWNER: ${{ github.repository_owner }} REPO: ${{ github.event.repository.name }} run: | PREVIOUS_JOB_ID=$(jq -r '.id' <<< "$WORKFLOW_RUN_EVENT_OBJ") - echo "Previous Job ID: $PREVIOUS_JOB_ID" echo "PREVIOUS_JOB_ID=$PREVIOUS_JOB_ID" >> "$GITHUB_ENV" SUITE_ID=$(jq -r '.check_suite_id' <<< "$WORKFLOW_RUN_EVENT_OBJ") - echo "Previous Suite ID: $SUITE_ID" echo "SUITE_ID=$SUITE_ID" >> "$GITHUB_ENV" ARTIFACT_ID=$(gh api "/repos/$OWNER/$REPO/actions/artifacts" \ - --jq ".artifacts.[] | - select(.workflow_run.id==${PREVIOUS_JOB_ID}) | - select(.expired==false) | - .id") - - echo "Artifact ID: $ARTIFACT_ID" + --jq ".artifacts[] | select(.workflow_run.id==$PREVIOUS_JOB_ID and .expired==false and .name=='detailedEvolutionPlot.png') | .id" | head -n1) echo "ARTIFACT_ID=$ARTIFACT_ID" >> "$GITHUB_ENV" - PR_NUMBER=$(jq -r '.pull_requests[0].number' \ - <<< "$WORKFLOW_RUN_EVENT_OBJ") - - echo "Pull request Number: $PR_NUMBER" - echo "PR_NUMBER=$PR_NUMBER" >> "$GITHUB_ENV" + PR_NUMBER=$(gh api "/repos/$OWNER/$REPO/commits/${{ github.event.workflow_run.head_sha }}/pulls" --jq '.[0].number') + echo "PR_NUMBER=${PR_NUMBER:-null}" >> "$GITHUB_ENV" - HEAD_SHA=$(jq -r '.pull_requests[0].head.sha' \ - <<< "$WORKFLOW_RUN_EVENT_OBJ") - - echo "Head SHA: $HEAD_SHA" + HEAD_SHA=${{ github.event.workflow_run.head_sha }} echo "HEAD_SHA=$HEAD_SHA" >> "$GITHUB_ENV" - - name: Find Comment - uses: peter-evans/find-comment@v2 - id: find-comment - with: - issue-number: ${{ env.PR_NUMBER }} - comment-author: 'github-actions[bot]' + + - name: Download artifact + if: env.PR_NUMBER != 'null' && env.ARTIFACT_ID != '' + run: | + gh api /repos/${{ github.repository }}/actions/artifacts/${{ env.ARTIFACT_ID }}/zip -H "Accept: application/vnd.github.v3+json" --output artifact.zip + unzip artifact.zip + rm artifact.zip + - name: Update Comment + if: env.PR_NUMBER != 'null' && env.ARTIFACT_ID != '' env: JOB_PATH: "${{ github.server_url }}/${{ github.repository }}/actions/runs/${{ env.PREVIOUS_JOB_ID }}" - ARTIFACT_URL: "${{ github.server_url }}/${{ github.repository }}/suites/$SUITE_ID/artifacts/$ARTIFACT_ID" - HEAD_SHA: "${{ env.HEAD_SHA }}" - uses: peter-evans/create-or-update-comment@v3 + ARTIFACT_URL: "${{ github.server_url }}/${{ github.repository }}/suites/${{ env.SUITE_ID }}/artifacts/${{ env.ARTIFACT_ID }}" + uses: actions/github-script@v6 with: - issue-number: ${{ env.PR_NUMBER }} - comment-id: ${{ steps.find-comment.outputs.comment-id }} - edit-mode: replace - body: |- - ![badge] - - Build Successful! You can find a link to the downloadable artifact below. - - | Name | Link | - | -------- | ----------------------- | - | Commit | ${{ env.HEAD_SHA }} | - | Logs | ${{ env.JOB_PATH }} | - | Download | ${{ env.ARTIFACT_URL }} | - - [badge]: https://img.shields.io/badge/Build_Success!-0d1117?style=for-the-badge&labelColor=3fb950&logo=data:image/svg%2bxml;base64,PHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHZpZXdCb3g9IjAgMCAyNCAyNCIgd2lkdGg9IjI0IiBoZWlnaHQ9IjI0IiBmaWxsPSIjZmZmZmZmIj48cGF0aCBkPSJNMjEuMDMgNS43MmEuNzUuNzUgMCAwIDEgMCAxLjA2bC0xMS41IDExLjVhLjc0Ny43NDcgMCAwIDEtMS4wNzItLjAxMmwtNS41LTUuNzVhLjc1Ljc1IDAgMSAxIDEuMDg0LTEuMDM2bDQuOTcgNS4xOTVMMTkuOTcgNS43MmEuNzUuNzUgMCAwIDEgMS4wNiAwWiI+PC9wYXRoPjwvc3ZnPg== + github-token: ${{secrets.GITHUB_TOKEN}} + script: | + const fs = require('fs'); + const imageData = fs.readFileSync('detailedEvolutionPlot.png', {encoding: 'base64'}); + github.rest.issues.createComment({ + issue_number: process.env.PR_NUMBER, + owner: context.repo.owner, + repo: context.repo.name, + body: ` + ![badge] + + Build Successful! You can find a link to the downloadable artifact below. + + | Name | Link | + | -------- | ----------------------- | + | Commit | ${process.env.HEAD_SHA} | + | Logs | ${process.env.JOB_PATH} | + | Download | ${process.env.ARTIFACT_URL} | + + ### Detailed Evolution Plot + ![Detailed Evolution Plot](data:image/png;base64,${imageData}) + + [badge]: https://img.shields.io/badge/Build_Success!-0d1117?style=for-the-badge&labelColor=3fb950&logo=data:image/svg%2bxml;base64,PHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHZpZXdCb3g9IjAgMCAyNCAyNCIgd2lkdGg9IjI0IiBoZWlnaHQ9IjI0IiBmaWxsPSIjZmZmZmZmIj48cGF0aCBkPSJNMjEuMDMgNS43MmEuNzUuNzUgMCAwIDEgMCAxLjA2bC0xMS41IDExLjVhLjc0Ny43NDcgMCAwIDEtMS4wNzItLjAxMmwtNS41LTUuNzVhLjc1Ljc1IDAgMSAxIDEuMDg0LTEuMDM2bDQuOTcgNS4xOTVMMTkuOTcgNS43MmEuNzUuNzUgMCAwIDEgMS4wNiAwWiI+PC9wYXRoPjwvc3ZnPg== + ` + }) diff --git a/.github/workflows/precommit-checks.yml b/.github/workflows/precommit-checks.yml index bc7f4762c..66885dc15 100644 --- a/.github/workflows/precommit-checks.yml +++ b/.github/workflows/precommit-checks.yml @@ -1,16 +1,24 @@ -name: precommit checks +name: Pre-commit Checks on: push: + pull_request: + +# Ensures only the latest workflow run for the same branch is active, canceling any in-progress runs. +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true jobs: pre-commit: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 - - uses: actions/setup-python@v4 + - uses: actions/setup-python@v5 with: - python-version: 3.8 + python-version: '3.9' - uses: pre-commit/action@v3.0.0 + - uses: pre-commit-ci/lite-action@v1.0.1 + if: always() diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 79a543313..73a323f2b 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -22,7 +22,6 @@ formats: [] # Set the version of Python and other tools you might need python: install: - - requirements: requirements.txt - method: pip path: . extra_requirements: diff --git a/online-docs/.readthedocs.yaml b/online-docs/.readthedocs.yaml deleted file mode 100644 index d616bd499..000000000 --- a/online-docs/.readthedocs.yaml +++ /dev/null @@ -1,15 +0,0 @@ -version: 1 - -build: - os: "ubuntu-20.04" - tools: - python: "3.12" - -# Build from the docs/ directory with Sphinx -sphinx: - configuration: conf.py - -# Explicitly set the version of Python and its requirements -python: - install: - - requirements: requirements.txt 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 625c25782..e58a9c23a 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 @@ -857,7 +857,7 @@ Default = ISOTROPIC **--mass-transfer-fa** |br| Mass Transfer fraction accreted (beta). |br| -Used when ``--mass-transfer-accretion-efficiency-prescription = FIXED_FRACTION``. |br| +Used when ``--mass-transfer-accretion-efficiency-prescription = FIXED``. |br| Default = 0.5 **--mass-transfer-jloss** |br| diff --git a/online-docs/requirements.txt b/online-docs/requirements.txt deleted file mode 100644 index 576eb06b0..000000000 --- a/online-docs/requirements.txt +++ /dev/null @@ -1,10 +0,0 @@ -sphinx>=7.0.1 -sphinx_rtd_theme>=1.2.2 -readthedocs-sphinx-search>=0.1.1b -sphinx-togglebutton -linuxdoc -sphinx_math_dollar -nbsphinx -numpydoc -sphinx-argparse -sphinx_tabs diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index df55164e3..000000000 --- a/requirements.txt +++ /dev/null @@ -1,49 +0,0 @@ -# File: docs/requirements.txt - -# Defining the exact version will make sure things don't break - -# sphinx is required. readthedocs does that automatically, so it doesn't need -# to be here but you will need it if you want to build the html files locally - -# seems obvious, but python is required. readthedocs does that automatically, -# so it doesn't need to be here but you will need it if you want to build the -# html files locally - -sphinx_rtd_theme==0.5.1 -sphinx-math-dollar>=1.2 - -sphinxcontrib.bibtex - -linuxdoc==20210324 - - -numpydoc -sphinx-tabs -sphinx-argparse - -nbsphinx - -ipykernel -IPython -pandas -numpy>=1.16 -# numpydoc -astropy>=4.0 -# numba>=0.50 -scipy>=1.5.0 -h5py -matplotlib>=3.3.2 -latex -PyYAML - -# The following are listed in the 'extensions []' block in conf.py -# Readthedocs doesn't want them here, but they are required if you want -# to build the html files locally -# -# matplotlib.sphinxext.plot_directive -# sphinx.ext.mathjax -# sphinx.ext.intersphinx -# matplotlib.sphinxext.plot_directive -# IPython.sphinxext.ipython_console_highlighting -# IPython.sphinxext.ipython_directive - diff --git a/setup.py b/setup.py index cfc45a311..8cf5abe20 100644 --- a/setup.py +++ b/setup.py @@ -24,7 +24,7 @@ "Programming Language :: Python :: 3", ] INSTALL_REQUIRES = [ - "numpy", + "numpy>=1.16", "h5py", "argparse", "stroopwafel", @@ -33,10 +33,10 @@ "flake8", "black==22.10.0", "isort", - "matplotlib", + "matplotlib>=3.3.2", "pandas", - "astropy", - "scipy", + "astropy>=4.0", + "scipy>=1.5.0", "latex", "PyYAML", "tqdm", @@ -44,16 +44,19 @@ ] EXTRA_REQUIRE = dict( docs=[ - "sphinx", + "sphinx>=7.0.1", "numpydoc", "nbsphinx", - "sphinx_rtd_theme", + "sphinx_rtd_theme>=1.2.2", "sphinx-tabs", "sphinx-argparse", - "sphinx-math-dollar", + "sphinx-math-dollar>=1.2", "sphinxcontrib.bibtex", "linuxdoc", - "ipython" + "ipython", + "readthedocs-sphinx-search>=0.1.1b", + "sphinx-togglebutton", + "linuxdoc>=20210324" ], dev=[ "pytest-cov", diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index f53d26483..892b15f1b 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -224,9 +224,14 @@ BaseBinaryStar::BaseBinaryStar(const unsigned long int p_Seed, const long int p_ if (!done) error = ERROR::INVALID_INITIAL_ATTRIBUTES; // too many iterations - bad initial conditions - if (error != ERROR::NONE) THROW_ERROR(error); // throw error if necessary - - SetRemainingValues(); // complete the construction of the binary + if (error != ERROR::NONE) { // ok? + m_EvolutionStatus = EVOLUTION_STATUS::BINARY_ERROR; // set evolutionary status + (void)PrintBinarySystemParameters(); // no - print (log) binary system parameters + THROW_ERROR(error); // throw error - can't return it... + } + else { // yes - ok + SetRemainingValues(); // complete the construction of the binary + } } @@ -296,18 +301,14 @@ void BaseBinaryStar::SetRemainingValues() { // check for CHE // - // because we've changed the rotational frequency of the constituent stars we - // have to reset the stellar type - at this stage, based on their rotational - // frequency at birth, they may have already been assigned one of MS_LTE_07, - // MS_GT_07, or CHEMICALLY_HOMOGENEOUS - // - // here we need to change from MS_* -> CH, or from CH->MS* based on the - // newly-assigned rotational frequencies + // here we need to change from MS_* -> CH, or from CH->MS* based on the binary orbital frequency, assuming that the stars are tidally locked + // set the spin to the orbital frequency, unless the user has specified a spin frequency // star 1 - if (utils::Compare(m_Star1->Omega(), m_Star1->OmegaCHE()) >= 0) { // star 1 CH? - if (m_Star1->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) (void)m_Star1->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous - m_Star1->SetOmega(omega); + if (utils::Compare(omega, m_Star1->OmegaCHE()) >= 0) { // star 1 CH? + if (m_Star1->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { + (void)m_Star1->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous + if (OPTIONS->OptionDefaulted("rotational-frequency-1")) m_Star1->SetOmega(omega); } // set spin to orbital frequency unless user specified } else if (m_Star1->MZAMS() <= 0.7) { // no - MS - initial mass determines actual type (don't use utils::Compare() here) if (m_Star1->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star1->SwitchTo(STELLAR_TYPE::MS_LTE_07, true); // MS <= 0.7 Msol - switch if necessary @@ -317,9 +318,10 @@ void BaseBinaryStar::SetRemainingValues() { } // star 2 - if (utils::Compare(m_Star2->Omega(), m_Star2->OmegaCHE()) >= 0) { // star 2 CH? - if (m_Star2->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) (void)m_Star2->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous - m_Star2->SetOmega(omega); + if (utils::Compare(omega, m_Star2->OmegaCHE()) >= 0) { // star 2 CH? + if (m_Star2->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { + (void)m_Star2->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous + if (OPTIONS->OptionDefaulted("rotational-frequency-2")) m_Star2->SetOmega(omega); } // set spin to orbital frequency unless user specified } else if (m_Star2->MZAMS() <= 0.7) { // no - MS - initial mass determines actual type (don't use utils::Compare() here) if (m_Star2->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star2->SwitchTo(STELLAR_TYPE::MS_LTE_07, true); // MS <= 0.0 Msol - switch if necessary @@ -1047,7 +1049,8 @@ double BaseBinaryStar::CalculateDEccentricityTidalDt(const DBL_DBL_DBL_DBL p_ImK double R1_AU = radiusStar * RSOL_TO_AU; double R1_over_a = R1_AU / m_SemiMajorAxis; double R1_over_a_8 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; - + + // No need to ignore quadratic e order terms during (super) synchronous rotation, since this formula is already linear in eccentricity return -(3.0 / 4.0) * (m_Eccentricity / OrbitalAngularVelocity()) * (1.0 + (massCompanion / massStar)) * (G_AU_Msol_yr * massCompanion / R1_AU / R1_AU / R1_AU) * R1_over_a_8 * ((3.0 * ImK10 / 2.0) - (ImK12 / 4.0) - ImK22 + (49.0 * ImK32 / 4.0)); } @@ -1075,8 +1078,11 @@ double BaseBinaryStar::CalculateDOmegaTidalDt(const DBL_DBL_DBL_DBL p_ImKlm, con double R1_AU = radiusStar * RSOL_TO_AU; double R1_over_a = R1_AU / m_SemiMajorAxis; double R1_over_a_6 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; + double e2_spin_term = (m_Eccentricity * m_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)); - return (3.0 / 2.0) * (1.0 / MoIstar) * (G_AU_Msol_yr * massCompanion * massCompanion / R1_AU) * R1_over_a_6 * (ImK22 + ((m_Eccentricity * m_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)))); + // if the star is rotating (super) synchronously AND quadratic 'e' terms cause the star to spin up further, ignore the higher order terms + if ((utils::Compare(p_Star->Omega(), OrbitalAngularVelocity()) > 0) && (utils::Compare((ImK22 + e2_spin_term), 0.0) > 0)){e2_spin_term = 0.0;} + return (3.0 / 2.0) * (1.0 / MoIstar) * (G_AU_Msol_yr * massCompanion * massCompanion / R1_AU) * R1_over_a_6 * (ImK22 + e2_spin_term); } @@ -1103,8 +1109,15 @@ double BaseBinaryStar::CalculateDSemiMajorAxisTidalDt(const DBL_DBL_DBL_DBL p_Im double R1_AU = radiusStar * RSOL_TO_AU; double R1_over_a = R1_AU / m_SemiMajorAxis; double R1_over_a_7 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; + double e2_sma_term = (m_Eccentricity * m_Eccentricity) * ((3.0 * ImK10 / 4.0) + (ImK12 / 8.0) - (5.0 * ImK22) + (147.0 * ImK32 / 8.0)); + + // if the star is rotating (super) synchronously AND quadratic 'e' terms cause the star to spin up further, ignore the higher order terms. + // Note: here we use the SPIN e^2 terms (not the semi-major axis terms) to determine when to ignore the higher order terms in semi-major axis evolution. + // this is to ensure that the higher order terms are always consistently applied/ignored across the tidal evolution equations. + double e2_spin_term = (m_Eccentricity * m_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)); + if ((utils::Compare(p_Star->Omega(), OrbitalAngularVelocity()) > 0) && (utils::Compare((ImK22 + e2_spin_term), 0.0) > 0)){e2_sma_term = 0.0;} - return -(3.0 / OrbitalAngularVelocity()) * (1.0 + (massCompanion / massStar)) * (G_AU_Msol_yr * massCompanion / R1_AU / R1_AU) * R1_over_a_7 * (ImK22 + ((m_Eccentricity * m_Eccentricity) * ((3.0 * ImK10 / 4.0) + (ImK12 / 8.0) - (5.0 * ImK22) + (147.0 * ImK32 / 8.0)))); + return -(3.0 / OrbitalAngularVelocity()) * (1.0 + (massCompanion / massStar)) * (G_AU_Msol_yr * massCompanion / R1_AU / R1_AU) * R1_over_a_7 * (ImK22 + e2_sma_term); } @@ -1307,19 +1320,23 @@ void BaseBinaryStar::ResolveSupernova() { // then the cross product is not well-defined, and we need to account for degeneracy between eccentricity vectors. // Also, if either eccentricity is 0.0, then the eccentricity vector is not well defined. - if ((utils::Compare(m_ThetaE, 0.0) == 0) && // orbitalAngularMomentumVectorPrev parallel to orbitalAngularMomentumVector? - ((utils::Compare(eccentricityPrev, 0.0) > 0) && (utils::Compare(m_Eccentricity, 0.0) > 0))) { // ... and both eccentricityVectorPrev and eccentricityVector well-defined? - - double psiPlusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); // yes - then psi + phi is constant - m_PhiE = _2_PI * RAND->Random(); - m_PsiE = psiPlusPhi - m_PhiE; - } - else if ((utils::Compare(m_ThetaE, M_PI) == 0) && // orbitalAngularMomentumVectorPrev anti-parallel to orbitalAngularMomentumVector? - ((utils::Compare(eccentricityPrev, 0.0) > 0) && (utils::Compare(m_Eccentricity, 0.0) > 0))) { // ... and both eccentricityVectorPrev and eccentricityVector well-defined? - - double psiMinusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); // yes - then psi - phi is constant - m_PhiE = _2_PI * RAND->Random(); - m_PsiE = psiMinusPhi + m_PhiE; + if ((utils::Compare(m_ThetaE, 0.0) == 0) || (utils::Compare(m_ThetaE, M_PI) == 0)) { // orbitalAngularMomentumVectorPrev parallel or anti-parallel to orbitalAngularMomentumVector + if ((utils::Compare(eccentricityPrev, 0.0) == 0) || (utils::Compare(m_Eccentricity, 0.0) == 0)) { // either e_prev or e_now is 0, so eccentricity vector is not well-defined + m_PhiE = _2_PI * RAND->Random(); + m_PsiE = _2_PI * RAND->Random(); + } + else { // both eccentricityVectorPrev and eccentricityVector well-defined + if (utils::Compare(m_ThetaE, 0.0) == 0){ // Orbital AM is parallel ? + double psiPlusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); // yes - then psi + phi is constant + m_PhiE = _2_PI * RAND->Random(); + m_PsiE = psiPlusPhi - m_PhiE; + } + else { + double psiMinusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); // no - then psi - phi is constant + m_PhiE = _2_PI * RAND->Random(); + m_PsiE = psiMinusPhi + m_PhiE; + } + } } else { // neither - the cross product of the orbit normals is well-defined Vector3d orbitalPivotAxis = cross(orbitalAngularMomentumVectorPrev, orbitalAngularMomentumVector); // cross product of the orbit normals @@ -1644,7 +1661,7 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { } } - if (utils::Compare(m_SemiMajorAxis, 0.0) <= 0 || utils::Compare(m_Star1->CalculateRemnantRadius() + m_Star2->CalculateRemnantRadius(), m_SemiMajorAxis * AU_TO_RSOL) > 0) { // catch merger in CE here, do not update stars + if (utils::Compare(m_SemiMajorAxis, 0.0) <= 0 || utils::Compare(m_Star1->CalculateRemnantRadius() + m_Star2->CalculateRemnantRadius(), m_SemiMajorAxis * AU_TO_RSOL) > 0) { // catch merger in CE here, do not update stars m_MassTransferTrackerHistory = MT_TRACKING::MERGER; m_Flags.stellarMerger = true; } @@ -1676,13 +1693,6 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::STELLAR_TYPE_CHANGE_DURING_CEE); // yes - print (log) detailed output } - // if stars are evolving as CHE stars, update their rotational frequency under the assumption of tidal locking if tides are not enabled - if (!m_Flags.stellarMerger && OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::NONE) { - double omega = OrbitalAngularVelocity(); // orbital angular velocity - if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega); - if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega); - } - m_Star1->SetPostCEEValues(); // squirrel away post CEE stellar values for star 1 m_Star2->SetPostCEEValues(); // squirrel away post CEE stellar values for star 2 SetPostCEEValues(m_SemiMajorAxis * AU_TO_RSOL, m_Eccentricity, rRLdfin1Rsol, rRLdfin2Rsol); // squirrel away post CEE binary values (checks for post-CE RLOF, so should be done at end) @@ -1939,8 +1949,8 @@ void BaseBinaryStar::CalculateWindsMassLoss() { double mWinds1 = m_Star1->CalculateMassLossValues(true); // calculate new values assuming mass loss applied double mWinds2 = m_Star2->CalculateMassLossValues(true); // calculate new values assuming mass loss applied - double aWinds = m_SemiMajorAxisPrev / (2.0 - ((m_Star1->MassPrev() + m_Star2->MassPrev()) / (mWinds1 + mWinds2))); // new semi-major axis for circularlised orbit - + double aWinds = m_SemiMajorAxisPrev * (m_Star1->Mass() + m_Star2->Mass()) / (mWinds1 + mWinds2); // new semi-major axis after wind mass loss, integrated to ensure a*M conservation + m_Star1->SetMassLossDiff(mWinds1 - m_Star1->Mass()); // JR: todo: find a better way? m_Star2->SetMassLossDiff(mWinds2 - m_Star2->Mass()); // JR: todo: find a better way? @@ -2403,7 +2413,7 @@ double BaseBinaryStar::CalculateAngularMomentum(const double p_SemiMajorAxis, double Jorb = CalculateOrbitalAngularMomentum(p_Star1Mass, p_Star2Mass, p_SemiMajorAxis, p_Eccentricity); - return (p_Star1MomentOfInertia * p_Star1SpinAngularVelocity) + (p_Star2MomentOfInertia * p_Star2SpinAngularVelocity) + Jorb; + return (p_Star1MomentOfInertia * p_Star1SpinAngularVelocity) + (p_Star2MomentOfInertia * p_Star2SpinAngularVelocity) + Jorb; } @@ -2425,8 +2435,8 @@ void BaseBinaryStar::CalculateEnergyAndAngularMomentum() { m_OrbitalAngularMomentumPrev = m_OrbitalAngularMomentum; m_TotalAngularMomentumPrev = m_TotalAngularMomentum; - double totalMass = m_Star1->Mass() + m_Star2->Mass(); - double reducedMass = (m_Star1->Mass() * m_Star2->Mass()) / totalMass; + double totalMass = m_Star1->Mass() + m_Star2->Mass(); + double reducedMass = (m_Star1->Mass() * m_Star2->Mass()) / totalMass; m_OrbitalEnergy = CalculateOrbitalEnergy(reducedMass, totalMass, m_SemiMajorAxis); m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); @@ -2464,9 +2474,9 @@ void BaseBinaryStar::ResolveMassChanges() { if (utils::Compare(massChange, 0.0) != 0) { // winds/mass transfer changes mass? // yes - calculate new angular momentum; assume accretor is adding angular momentum from a circular orbit at the stellar radius - double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) ? - massChange * sqrt(G_AU_Msol_yr * m_Star1->Mass() * m_Star1->Radius() * RSOL_TO_AU) : - (2.0 / 3.0) * massChange * m_Star1->Radius() * RSOL_TO_AU * m_Star1->Radius() * RSOL_TO_AU * m_Star1->Omega(); + double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) + ? massChange * sqrt(G_AU_Msol_yr * m_Star1->Mass() * m_Star1->Radius() * RSOL_TO_AU) + : (2.0 / 3.0) * massChange * m_Star1->Radius() * RSOL_TO_AU * m_Star1->Radius() * RSOL_TO_AU * m_Star1->Omega(); // update mass of star according to mass loss and mass transfer, then update age accordingly (void)m_Star1->UpdateAttributes(massChange, 0.0); // update mass for star m_Star1->UpdateInitialMass(); // update effective initial mass of star (MS, HG & HeMS) @@ -2487,9 +2497,9 @@ void BaseBinaryStar::ResolveMassChanges() { double massChange = m_Star2->MassLossDiff() + m_Star2->MassTransferDiff(); // mass change due to winds and mass transfer if (utils::Compare(massChange, 0.0) != 0) { // winds/mass transfer changes mass? // yes - calculate new angular momentum; assume accretor is adding angular momentum from a circular orbit at the stellar radius - double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) ? - massChange * sqrt(G_AU_Msol_yr * m_Star2->Mass() * m_Star2->Radius() * RSOL_TO_AU) : - (2.0 / 3.0) * massChange * m_Star2->Radius() * RSOL_TO_AU * m_Star2->Radius() * RSOL_TO_AU * m_Star2->Omega(); + double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) + ? massChange * sqrt(G_AU_Msol_yr * m_Star2->Mass() * m_Star2->Radius() * RSOL_TO_AU) + : (2.0 / 3.0) * massChange * m_Star2->Radius() * RSOL_TO_AU * m_Star2->Radius() * RSOL_TO_AU * m_Star2->Omega(); // update mass of star according to mass loss and mass transfer, then update age accordingly (void)m_Star2->UpdateAttributes(massChange, 0.0); // update mass for star m_Star2->UpdateInitialMass(); // update effective initial mass of star (MS, HG & HeMS) @@ -2515,14 +2525,6 @@ void BaseBinaryStar::ResolveMassChanges() { m_Star2->ResetEnvelopeExpulsationByPulsations(); } - // if CHE enabled, update rotational frequency for constituent stars - assume tidally locked if tidal mode is none (otherwise, let tides do the work) - if(OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::NONE) { - double omega = OrbitalAngularVelocity(); // orbital angular velocity - if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega); - if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega); - } - - CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations if ((m_Star1->StellarType() != stellarType1) || (m_Star2->StellarType() != stellarType2)) { // stellar type change? @@ -2543,15 +2545,45 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { if (!m_Unbound) { // binary bound? // yes - process tides if enabled - double omega = OrbitalAngularVelocity(); switch (OPTIONS->TidesPrescription()) { // which tides prescription? case TIDES_PRESCRIPTION::NONE: { // NONE - tides not enabled - // do nothing, except for CHE stars which are allowed to remain CHE even though this does not preserve angular momentum - if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega); - if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega); + // do nothing, except for CHE stars which are allowed to remain CHE + + // if at least one star is CHE, then circularize the binary and synchronize only the CHE stars conserving total angular momentum + if (OPTIONS->CHEMode() != CHE_MODE::NONE && HasOneOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS})) { // one CHE star? + double che_I1 = 0.0; + double che_I2 = 0.0; + double che_Ltot = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); + + if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { + che_I1 = m_Star1->CalculateMomentOfInertiaAU(); + che_Ltot += m_Star1->AngularMomentum(); + } + + if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { + che_I2 = m_Star2->CalculateMomentOfInertiaAU(); + che_Ltot += m_Star2->AngularMomentum(); + } + + double omega_sync = OmegaAfterSynchronisation(m_Star1->Mass(), m_Star2->Mass(), che_I1, che_I2, che_Ltot, omega); + if (omega_sync >= 0.0) { // root found? (don't use utils::Compare() here) + // yes + if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star1->SetOmega(omega_sync);} + if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star2->SetOmega(omega_sync);} + + m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega_sync / omega_sync); // re-calculate semi-major axis + m_Eccentricity = 0.0; // circularise + m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum + m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); + } + else { // no (real) root found, synchronize CHE stars ignoring angular momentum conservation + if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star1->SetOmega(omega);} + if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star2->SetOmega(omega);} + } + } } break; @@ -2577,6 +2609,7 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { m_SemiMajorAxis = m_SemiMajorAxis + ((DSemiMajorAxis1Dt + DSemiMajorAxis2Dt) * p_Dt * MYR_TO_YEAR); // evolve separation m_Eccentricity = m_Eccentricity + ((DEccentricity1Dt + DEccentricity2Dt) * p_Dt * MYR_TO_YEAR); // evolve eccentricity m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum + m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); } break; @@ -2592,9 +2625,10 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { m_Star1->SetOmega(omega); // synchronise star 1 m_Star2->SetOmega(omega); // synchronise star 2 - m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega / omega); // re-calculate semi-major axis - m_Eccentricity = 0.0; // circularise - m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum + m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega / omega); // re-calculate semi-major axis + m_Eccentricity = 0.0; // circularise + m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum + m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); } else { // no (real) root found @@ -2637,6 +2671,7 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { } } + /* * Root solver to determine rotational frequency after synchronisation for tides * @@ -2740,7 +2775,6 @@ double BaseBinaryStar::OmegaAfterSynchronisation(const double p_M1, const double } - /* * Calculate and emit gravitational radiation. * @@ -2916,8 +2950,8 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { } } - if (!StellarMerger()) { // stellar merger? - // no - continue evolution + 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())) { 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 @@ -3171,23 +3205,25 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { } } - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::PRE_STELLAR_TIMESTEP); // print (log) detailed output - if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? // yes + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::PRE_STELLAR_TIMESTEP); // print (log) detailed output + error = EvolveOneTimestep(dt); // evolve the binary system one timestep if (error != ERROR::NONE) { // SSE error for either constituent star? evolutionStatus = EVOLUTION_STATUS::SSE_ERROR; // yes - stop evolution } - else { // continue evolution + } - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::TIMESTEP_COMPLETED); // print (log) detailed output: this is after all changes made in the timestep + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::TIMESTEP_COMPLETED); // print (log) detailed output: this is after all changes made in the timestep - if (stepNum >= OPTIONS->MaxNumberOfTimestepIterations()) evolutionStatus = EVOLUTION_STATUS::STEPS_UP; // number of timesteps for evolution exceeds maximum - else if (evolutionStatus == EVOLUTION_STATUS::CONTINUE && usingProvidedTimesteps && stepNum >= timesteps.size()) { // using user-provided timesteps and all consumed - evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED; // yes - set status - SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED); // show warning - } + + if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? + // yes + if (stepNum >= OPTIONS->MaxNumberOfTimestepIterations()) evolutionStatus = EVOLUTION_STATUS::STEPS_UP; // number of timesteps for evolution exceeds maximum + else if (evolutionStatus == EVOLUTION_STATUS::CONTINUE && usingProvidedTimesteps && stepNum >= timesteps.size()) { // using user-provided timesteps and all consumed + evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED; // yes - set status + SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED); // show warning } if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 6152004c4..590cf1727 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -4533,7 +4533,7 @@ STELLAR_TYPE BaseStar::UpdateAttributesAndAgeOneTimestep(const double p_DeltaMas stellarType = STELLAR_TYPE::MASSLESS_REMNANT; } else { - stellarType = ResolveSupernova(); // handle supernova + stellarType = ResolveSupernova(); // handle supernova if (stellarType == m_StellarType) { // still on phase? UpdateAttributesAndAgeOneTimestepPreamble(p_DeltaMass, p_DeltaMass0, p_DeltaTime); // apply mass changes and save current values if required diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index fb48a2e0b..df159baad 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -1991,7 +1991,7 @@ STELLAR_TYPE GiantBranch::ResolveCoreCollapseSN() { /* - * Resolve Electron capture Supernova + * Resolve Electron Capture Supernova * * Calculate the mass of the remnant and set remnant type - always a Neutron Star * Updates attributes of star; sets SN flags @@ -2007,39 +2007,16 @@ STELLAR_TYPE GiantBranch::ResolveCoreCollapseSN() { * @return Stellar type of remnant (always STELLAR_TYPE::NEUTRON_STAR) */ STELLAR_TYPE GiantBranch::ResolveElectronCaptureSN() { - - STELLAR_TYPE stellarType = m_StellarType; // remnant stellar type - - if (!m_MassTransferDonorHistory.empty() || (OPTIONS->AllowNonStrippedECSN())) { // if progenitor has never been a MT donor, is it allowed to ECSN? - // yes - m_Mass = MECS_REM; // defined in constants.h - m_CoreMass = m_Mass; - m_HeCoreMass = m_Mass; - m_COCoreMass = m_Mass; - m_Mass0 = m_Mass; - - stellarType = STELLAR_TYPE::NEUTRON_STAR; - - SetSNCurrentEvent(SN_EVENT::ECSN); // electron capture SN happening now - SetSNPastEvent(SN_EVENT::ECSN); // ... and will be a past event - } - else { // not allowed to ECSN, treat as ONeWD + m_Mass = MECS_REM; // defined in constants.h + m_CoreMass = m_Mass; + m_HeCoreMass = m_Mass; + m_COCoreMass = m_Mass; + m_Mass0 = m_Mass; - if (utils::Compare(m_COCoreMass, MCH) > 0) { - SHOW_WARN(ERROR::WHITE_DWARF_TOO_MASSIVE, "Setting mass to Chandraskhar mass."); - } - m_Mass = std::min(m_COCoreMass, MCH); // no WD masses above Chandrasekhar mass - m_CoreMass = m_Mass; - m_HeCoreMass = m_Mass; - m_COCoreMass = m_Mass; - m_Mass0 = m_Mass; - m_Radius = WhiteDwarfs::CalculateRadiusOnPhase_Static(m_Mass); // radius is defined equivalently for all WDs - m_Luminosity = ONeWD::CalculateLuminosityOnPhase_Static(m_Mass, m_Time, m_Metallicity); // need to set the luminosity for ONeWD specifically - - stellarType = STELLAR_TYPE::OXYGEN_NEON_WHITE_DWARF; - } + SetSNCurrentEvent(SN_EVENT::ECSN); // electron capture SN happening now + SetSNPastEvent(SN_EVENT::ECSN); // ... and will be a past event - return stellarType; + return STELLAR_TYPE::NEUTRON_STAR; } @@ -2252,8 +2229,8 @@ STELLAR_TYPE GiantBranch::ResolveSupernova() { stellarType = ResolvePairInstabilitySN(); // MR } - else if (utils::Compare(snMass, MCBUR2) < 0) { // Electron Capture Supernova - stellarType = ResolveElectronCaptureSN(); // NS or ONeWD + else if (utils::Compare(snMass, MCBUR2) < 0 && (!m_MassTransferDonorHistory.empty() || OPTIONS->AllowNonStrippedECSN())) { + stellarType = ResolveElectronCaptureSN(); // electron capture SN; requires progenitor to have been a MT donor unless non-stripped ECSN are allowed; forms NS } else { // Core Collapse Supernova stellarType = ResolveCoreCollapseSN(); // BH or NS diff --git a/src/TPAGB.cpp b/src/TPAGB.cpp index fb137a6f7..a72465c10 100755 --- a/src/TPAGB.cpp +++ b/src/TPAGB.cpp @@ -926,7 +926,7 @@ STELLAR_TYPE TPAGB::ResolveEnvelopeLoss(bool p_Force) { if (ShouldEnvelopeBeExpelledByPulsations()) m_EnvelopeJustExpelledByPulsations = true; if (p_Force || (utils::Compare(m_CoreMass, m_Mass)) >= 0 || m_EnvelopeJustExpelledByPulsations) { // envelope loss - + m_Mass = std::min(m_CoreMass, m_Mass); m_CoreMass = m_Mass; m_HeCoreMass = m_Mass; @@ -967,3 +967,19 @@ double TPAGB::ChooseTimestep(const double p_Time) const { #undef timescales } + +/* + * Determine if star should continue evolution as a Supernova + * + * + * bool IsSupernova() + * + * @return Boolean flag: true if star has gone Supernova, false if not + */ +bool TPAGB::IsSupernova() const { + if(utils::SNEventType(m_SupernovaDetails.events.current) != SN_EVENT::NONE) + return true; // already labeled as going through a SN right now + double snMass = CalculateInitialSupernovaMass(); // calculate SN initial mass + return ( utils::Compare(m_COCoreMass, m_GBParams[static_cast(GBP::McSN)]) >=0 && utils::Compare(snMass, OPTIONS->MCBUR1()) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0 ); + // no supernova if CO core mass is too low or helium core mass is too low at base of AGB or the envelope has already been removed +} diff --git a/src/TPAGB.h b/src/TPAGB.h index 59b8f9f67..d979ff9fc 100755 --- a/src/TPAGB.h +++ b/src/TPAGB.h @@ -97,7 +97,7 @@ class TPAGB: virtual public BaseStar, public EAGB { STELLAR_TYPE EvolveToNextPhase() { return m_StellarType; } // NO-OP bool IsEndOfPhase() const { return !ShouldEvolveOnPhase(); } // Phase ends when envelope loss or going supernova - bool IsSupernova() const { return (utils::Compare(m_COCoreMass, m_GBParams[static_cast(GBP::McSN)]) >= 0 && utils::Compare(CalculateInitialSupernovaMass(), OPTIONS->MCBUR1()) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0); } // Going supernova if still has envelope and core mass large enough + bool IsSupernova() const; STELLAR_TYPE ResolveEnvelopeLoss(bool p_Force = false); void ResolveHeliumFlash() { } // NO-OP diff --git a/src/changelog.h b/src/changelog.h index 0bc9f88e4..d977012a1 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1416,12 +1416,32 @@ // - fixed issue with updating helium giants that manifested as supernovae with nan core mass (see #1245) // - added check for exceeding Chandrasekhar mass when computing white dwarf radius (resolves issue #1264) // - added check to only compute McBGB for stars with mass above MHeF, following text above Eq. 44 in Hurley+, 2000 (resolves issue #1256) -// 03.11.00 AB - Dec 04, 2024 - Enhancement: +// 03.10.02 IM - Dec 13, 2024 - Defect repair: +// - if the Hurley supernova criteria are met yet ECSN criteria based on mass transfer history are not met, a normal CCSN ensues as opposed to an ONeWD +// - exactly preserve the product of semi-major axis * total mass on wind mass loss +// 03.10.03 JR - Dec 16, 2024 - Defect repair: +// - fix for issue #1310 - run terminates prematurely if error in grid file +// 03.10.04 RTW - Nov 27, 2024 - Defect repair: +// - fix for issue #1247 - SN Euler angles had incomplete logic, leading to a div by zero in some cases +// 03.10.05 JR - Jan 08, 2025 - Defect repair: +// - fix for issue #1317 - SN events not always logged in BSE SN file when evolving MS merger products +// - added code to ensure final BSE detailed output file TIMESTEP_COMPLETED record is always logged +// (may duplicate FINAL_STATE record, but logging TIMESTEP_COMPLETED is consistent, and it's what most people look for) +// 03.10.06 VK - Jan 13, 2025 - Enhancement: +// - Modified the KAPIL2024 tides to ignore quadratic 'e' terms (for spin and separation evolution) if they spin up an already synchronized star. +// 03.11.00 VK - Jan 14, 2025 - Enhancement, Defect repair: +// - Fix for issue #1303 - Reduction in production of BHBH from CHE, other CHE-related improvements. +// - Stars that have sufficiently rapid angular frequencies at ZAMS are now initialized as CHE stars, regardless of the tidal prescription. +// - At CHE initialization, stellar spin is set to orbital frequency, unless rotational frequency has been specified by user. This process does not conserve angular momentum (implicitly assuming spin-up in the pre-ZAMS phase). +// - When checking for CHE, compare threshold frequency against orbit rather than stellar spin, in case the star has zero frequency (no tides, no user-specified value). +// - Moved all CHE rotation related code to ProcessTides(), ensuring that any spin up during binary evolution conserves total angular momentum. +// 03.12.00 AB - Jan 16, 2025 - Enhancement: // - Added Shikauchi et al. (2024) core mass prescription, describing convective core evolution under mass loss/gain // - New options: --main-sequence-core-mass-prescription SHIKAUCHI (new prescription), MANDEL (replaces --retain-core-mass-during-caseA-mass-transfer), -// NONE (no core mass treatment) -// - Added new luminosity prescription from Shikauchi et al. (2024) -// - Added treatment for rejuvenation of main sequence accretors -const std::string VERSION_STRING = "03.11.00"; +// ZERO (main sequence core mass set to zero, no treatment) +// - Added new luminosity prescription for main sequence stars from Shikauchi et al. (2024) +// - Added treatment for rejuvenation of main sequence accretors when the new prescription is used + +const std::string VERSION_STRING = "03.12.00"; # endif // __changelog_h__ diff --git a/src/main.cpp b/src/main.cpp index b688cdd95..08c1cee95 100755 --- a/src/main.cpp +++ b/src/main.cpp @@ -537,230 +537,247 @@ std::tuple EvolveBinaryStars() { EVOLUTION_STATUS evolutionStatus = EVOLUTION_STATUS::CONTINUE; - auto wallStart = std::chrono::system_clock::now(); // start wall timer - clock_t clockStart = clock(); // start CPU timer + auto wallStart = std::chrono::system_clock::now(); // start wall timer + clock_t clockStart = clock(); // start CPU timer std::time_t timeStart = std::chrono::system_clock::to_time_t(wallStart); SAY("Start generating binaries at " << std::ctime(&timeStart)); BinaryStar* binary = nullptr; - bool usingGrid = !OPTIONS->GridFilename().empty(); // using grid file? - size_t index = 0; // which binary + bool usingGrid = !OPTIONS->GridFilename().empty(); // using grid file? + size_t index = 0; // which binary - signal(SIGUSR1, SIGUSR1handler); // install SIGUSR1 signal handler + signal(SIGUSR1, SIGUSR1handler); // install SIGUSR1 signal handler // The options specified by the user at the commandline are set to their initial values. // OPTIONS->AdvanceCmdLineOptionValues(), called at the end of the loop, advances the // options specified by the user at the commandline to their next variation (if necessary, // based on any ranges and/or sets specified by the user). - while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // while all ok and not done + while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // while all ok and not done try { - // generate and evolve binaries + // generate and evolve binaries - int gridLineVariation = 0; // grid line variation number - bool doneGridFile = false; // flags we're done with the grid file (for this commandline variation) - bool processingGridLine = false; // processing a gridfile line? - while (!doneGridFile && evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // for each binary to be evolved + int gridLineVariation = 0; // grid line variation number + bool doneGridFile = false; // flags we're done with the grid file (for this commandline variation) + bool processingGridLine = false; // processing a gridfile line? + while (!doneGridFile && evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // for each binary to be evolved - if (OPTIONS->FPErrorMode() != FP_ERROR_MODE::OFF) { // floating-point error handling mode on/debug? - (void)feclearexcept(FE_ALL_EXCEPT); // yes - clear all FE traps - (void)feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); // enable FE traps (don't trap FE_INEXACT - would trap on almost all FP operations...) - } - else (void)fedisableexcept(FE_ALL_EXCEPT); // disable all FE traps - - evolvingBinaryStar = NULL; // unset global pointer to evolving binary (for BSE Switch Log) - evolvingBinaryStarValid = false; // indicate that the global pointer is not (yet) valid (for BSE Switch log) - - bool doneGridLine = false; // flags we're done with this grid file line (if using a grid file) - if (usingGrid) { // using grid file? - gridLineVariation = 0; // yes - first variation of this grid line - int gridResult = OPTIONS->ApplyNextGridLine(); // set options according to specified values in grid file - std::string errStr = "Accessing grid file '" + OPTIONS->GridFilename() + "'"; // common error string - switch (gridResult) { // handle result of grid file read - case -1: // error - unexpected end of file grid file read - case -2: { // error - read error for grid file - ERROR error = gridResult == -1 ? ERROR::UNEXPECTED_END_OF_FILE : ERROR::FILE_READ_ERROR; // set error - THROW_ERROR_STATIC(error, errStr); // throw error - } break; - case 0: { // end of file - doneGridLine = true; // flag we're done with this grid line - doneGridFile = true; // flag we're done with the grid file - ERROR error = OPTIONS->RewindGridFile(); // ready for next commandline options variation - if (error != ERROR::NONE) { // rewind ok? - THROW_ERROR_STATIC(error, errStr); // throw error - } - } break; - case 1: // grid line read - processingGridLine = true; // not done yet... - break; - default: // unexpected error - THROW_ERROR_STATIC(ERROR::ERROR, errStr); // throw error - break; + if (OPTIONS->FPErrorMode() != FP_ERROR_MODE::OFF) { // floating-point error handling mode on/debug? + (void)feclearexcept(FE_ALL_EXCEPT); // yes - clear all FE traps + (void)feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); // enable FE traps (don't trap FE_INEXACT - would trap on almost all FP operations...) } - } - else { // no, not using a grid file - doneGridFile = true; // flag we're done with the grid file - } - - // The options specified by the user in the grid file line are set to their initial values. - // (by OPTIONS->ApplyNextGridLine()). - // OPTIONS->AdvanceGridLineOptionValues(), called at the end of the loop, advances the - // options specified by the user in the grid file line to their next variation (if necessary, - // based on any ranges and/or sets specified by the user). - // Note that `doneGridLine` may be a proxy for the command line here (it is set FALSE even if - // there is no grid file - so the loop executes once to process the command line) - - while (!doneGridLine && (evolutionStatus == EVOLUTION_STATUS::CONTINUE)) { // while all ok and not done - - // we only need to pass the index number to the binary - we let the BinaryStar class do the work - // wrt setting the parameters for each of the constituent stars - // (The index is really only needed for legacy comparison, so can probably be removed at any time) - - // create the binary - - // Binary stars (in BSE) are provided with a random seed that is used to seed the random - // number generator. The random number generator is re-seeded for each binary. Here we - // generate the seed for the binary being evolved - by this point we have picked up the - // option value from either the commandline or the grid file. - // - // there are three scenarios: - // - // if the user did not specify a random seed, either on the commandline or in a grid file - // record, we use a randomly chosen seed, based on the system time. - // - // if the user specified a random seed on the commandline, and not in the grid file for - // the current binary, the random seed specified on the commandline is used - and the offset - // applied (the index of the binary being evolved). The index of the binary being evolved - // starts at 0 for the first binary, and increments by 1 for each subsequent binary evolved - // (so the base random seed specified by the user is also the initial random seed - the - // random seed of the first binary evolved) - // - // if the user specified a random seed in the grid file for the current binary, regardless of - // whether a random seed was specified on the commandline, the random seed from the grid - // file is used, and an offset is added if the grid line also specified ranges or sets for - // and options (if no ranges or sets were specified on the grid line then no offset is added - // (i.e. the random seed specified is used as is)). Note that in this scenario it is the - // user's responsibility to ensure that there is no duplication of seeds. - - unsigned long int thisId = index + gridLineVariation; // set the id for the binary - - std::string errorStr; // error string - unsigned long int randomSeed = 0l; // random seed - OPTIONS_ORIGIN optsOrigin = processingGridLine ? OPTIONS_ORIGIN::GRIDFILE : OPTIONS_ORIGIN::CMDLINE; // indicate which set of program options we're using - if (OPTIONS->FixedRandomSeedGridLine()) { // user specified a random seed in the grid file for this binary? - // yes - use it (indexed) - randomSeed = OPTIONS->RandomSeedGridLine() + (unsigned long int)gridLineVariation; // random seed - errorStr = OPTIONS->SetRandomSeed(randomSeed, optsOrigin); // set it - if (!errorStr.empty()) { // ok? - THROW_ERROR_STATIC(ERROR::ERROR_PROCESSING_GRIDLINE_OPTIONS, errorStr); // throw error - } - } - else if (OPTIONS->FixedRandomSeedCmdLine()) { // no - user specified a random seed on the commandline? - // yes - use it (indexed) - randomSeed = OPTIONS->RandomSeedCmdLine() + (unsigned long int)index + (unsigned long int)gridLineVariation; // random seed - errorStr = OPTIONS->SetRandomSeed(randomSeed, optsOrigin); // set it - if (!errorStr.empty()) { // ok? - THROW_ERROR_STATIC(ERROR::ERROR_PROCESSING_CMDLINE_OPTIONS, errorStr); // throw error + else (void)fedisableexcept(FE_ALL_EXCEPT); // disable all FE traps + + evolvingBinaryStar = NULL; // unset global pointer to evolving binary (for BSE Switch Log) + evolvingBinaryStarValid = false; // indicate that the global pointer is not (yet) valid (for BSE Switch log) + + bool doneGridLine = false; // flags we're done with this grid file line (if using a grid file) + if (usingGrid) { // using grid file? + gridLineVariation = 0; // yes - first variation of this grid line + int gridResult = OPTIONS->ApplyNextGridLine(); // set options according to specified values in grid file + std::string errStr = "Accessing grid file '" + OPTIONS->GridFilename() + "'"; // common error string + switch (gridResult) { // handle result of grid file read + case -1: // error - unexpected end of file grid file read + case -2: { // error - read error for grid file + ERROR error = gridResult == -1 ? ERROR::UNEXPECTED_END_OF_FILE : ERROR::FILE_READ_ERROR; // set error + THROW_ERROR_STATIC(error, errStr); // throw error + } break; + case 0: { // end of file + doneGridLine = true; // flag we're done with this grid line + doneGridFile = true; // flag we're done with the grid file + ERROR error = OPTIONS->RewindGridFile(); // ready for next commandline options variation + if (error != ERROR::NONE) { // rewind ok? + THROW_ERROR_STATIC(error, errStr); // throw error + } + } break; + case 1: // grid line read + processingGridLine = true; // not done yet... + break; + default: // unexpected error + THROW_ERROR_STATIC(ERROR::ERROR, errStr); // throw error + break; } } - else { // no - // use default seed (based on system time) + id (index) - randomSeed = RAND->DefaultSeed() + (unsigned long int)index + (unsigned long int)gridLineVariation; // random seed - errorStr = OPTIONS->SetRandomSeed(randomSeed, optsOrigin); // set it - if (!errorStr.empty()) { // ok? - THROW_ERROR_STATIC(ERROR::ERROR_PROCESSING_CMDLINE_OPTIONS, errorStr); // throw error - } + else { // no, not using a grid file + doneGridFile = true; // flag we're done with the grid file } - if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // ok? - // yes - continue - randomSeed = RAND->CurrentSeed(); // current random seed - to pass to binary object + // The options specified by the user in the grid file line are set to their initial values. + // (by OPTIONS->ApplyNextGridLine()). + // OPTIONS->AdvanceGridLineOptionValues(), called at the end of the loop, advances the + // options specified by the user in the grid file line to their next variation (if necessary, + // based on any ranges and/or sets specified by the user). + // Note that `doneGridLine` may be a proxy for the command line here (it is set FALSE even if + // there is no grid file - so the loop executes once to process the command line) + + while (!doneGridLine && (evolutionStatus == EVOLUTION_STATUS::CONTINUE)) { // while all ok and not done - delete binary; binary = nullptr; // so we don't leak - binary = new BinaryStar(randomSeed, thisId); // generate binary according to the user options + // we only need to pass the index number to the binary - we let the BinaryStar class do the work + // wrt setting the parameters for each of the constituent stars + // (The index is really only needed for legacy comparison, so can probably be removed at any time) - evolvingBinaryStar = binary; // set global pointer to evolving binary (for BSE Switch Log) - evolvingBinaryStarValid = true; // indicate that the global pointer is now valid (for BSE Switch Log) + // create the binary - EVOLUTION_STATUS thisBinaryStatus = binary->Evolve(); // evolve the binary - - // announce result of evolving the binary - if (!OPTIONS->Quiet()) { // quiet mode? - // no - announce result of evolving the binary - SAY(thisId << ": " << - EVOLUTION_STATUS_LABEL.at(thisBinaryStatus) << ": (" << - STELLAR_TYPE_LABEL.at(binary->Star1InitialType()) << " -> " << - STELLAR_TYPE_LABEL.at(binary->Star1Type()) << ") + (" << - STELLAR_TYPE_LABEL.at(binary->Star2InitialType()) << " -> " << - STELLAR_TYPE_LABEL.at(binary->Star2Type()) << ")" - ); - } + // Binary stars (in BSE) are provided with a random seed that is used to seed the random + // number generator. The random number generator is re-seeded for each binary. Here we + // generate the seed for the binary being evolved - by this point we have picked up the + // option value from either the commandline or the grid file. + // + // there are three scenarios: + // + // if the user did not specify a random seed, either on the commandline or in a grid file + // record, we use a randomly chosen seed, based on the system time. + // + // if the user specified a random seed on the commandline, and not in the grid file for + // the current binary, the random seed specified on the commandline is used - and the offset + // applied (the index of the binary being evolved). The index of the binary being evolved + // starts at 0 for the first binary, and increments by 1 for each subsequent binary evolved + // (so the base random seed specified by the user is also the initial random seed - the + // random seed of the first binary evolved) + // + // if the user specified a random seed in the grid file for the current binary, regardless of + // whether a random seed was specified on the commandline, the random seed from the grid + // file is used, and an offset is added if the grid line also specified ranges or sets for + // and options (if no ranges or sets were specified on the grid line then no offset is added + // (i.e. the random seed specified is used as is)). Note that in this scenario it is the + // user's responsibility to ensure that there is no duplication of seeds. + + unsigned long int thisId = index + gridLineVariation; // set the id for the binary - if (!LOGGING->CloseStandardFile(LOGFILE::BSE_DETAILED_OUTPUT)) { // close detailed output file if necessary - THROW_ERROR_STATIC(ERROR::FILE_NOT_CLOSED); // close failed - this will cause problems later - throw error + std::string errorStr; // error string + unsigned long int randomSeed = 0l; // random seed + OPTIONS_ORIGIN optsOrigin = processingGridLine ? OPTIONS_ORIGIN::GRIDFILE : OPTIONS_ORIGIN::CMDLINE; // indicate which set of program options we're using + if (OPTIONS->FixedRandomSeedGridLine()) { // user specified a random seed in the grid file for this binary? + // yes - use it (indexed) + randomSeed = OPTIONS->RandomSeedGridLine() + (unsigned long int)gridLineVariation; // random seed + errorStr = OPTIONS->SetRandomSeed(randomSeed, optsOrigin); // set it + if (!errorStr.empty()) { // ok? + THROW_ERROR_STATIC(ERROR::ERROR_PROCESSING_GRIDLINE_OPTIONS, errorStr); // throw error + } + } + else if (OPTIONS->FixedRandomSeedCmdLine()) { // no - user specified a random seed on the commandline? + // yes - use it (indexed) + randomSeed = OPTIONS->RandomSeedCmdLine() + (unsigned long int)index + (unsigned long int)gridLineVariation; // random seed + errorStr = OPTIONS->SetRandomSeed(randomSeed, optsOrigin); // set it + if (!errorStr.empty()) { // ok? + THROW_ERROR_STATIC(ERROR::ERROR_PROCESSING_CMDLINE_OPTIONS, errorStr); // throw error + } + } + else { // no + // use default seed (based on system time) + id (index) + randomSeed = RAND->DefaultSeed() + (unsigned long int)index + (unsigned long int)gridLineVariation; // random seed + errorStr = OPTIONS->SetRandomSeed(randomSeed, optsOrigin); // set it + if (!errorStr.empty()) { // ok? + THROW_ERROR_STATIC(ERROR::ERROR_PROCESSING_CMDLINE_OPTIONS, errorStr); // throw error + } } - ERRORS->Clean(); // clean the dynamic error catalog + if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // ok? + // yes - continue + randomSeed = RAND->CurrentSeed(); // current random seed - to pass to binary object - if (usingGrid) { // using grid file? - gridLineVariation++; // yes - increment grid line variation number - int optionsStatus = OPTIONS->AdvanceGridLineOptionValues(); // apply next grid file options (ranges/sets) - if (optionsStatus < 0) { // ok? - THROW_ERROR_STATIC(ERROR::ERROR_PROCESSING_GRIDLINE_OPTIONS); // no - throw error + EVOLUTION_STATUS thisBinaryStatus; // status of this binary + bool haveBinary = true; // binary constructed ok? + try { + delete binary; binary = nullptr; // so we don't leak + binary = new BinaryStar(randomSeed, thisId); // generate binary according to the user options } - else if (optionsStatus == 0) { // end of grid file options variations? - doneGridLine = true; // yes - we're done + catch (...) { // catchall + // if we catch an error here it happened while creating the binary - in the class constructor + // we don't need to halt the program here - just prevent evolution of the binary + // any error caught here has already been displayed to the user, so we just catch all + haveBinary = false; // binary not constructed + thisBinaryStatus = EVOLUTION_STATUS::BINARY_ERROR; // error with binary + + // announce result + if (!OPTIONS->Quiet()) { // quiet mode? + SAY(thisId << ": " << EVOLUTION_STATUS_LABEL.at(thisBinaryStatus) << ": not evolved"); // no - announce result of (not) evolving the binary + } } - } - else doneGridLine = true; // not using grid file - done + + if (haveBinary) { // binary constructed ok? + evolvingBinaryStar = binary; // yes - set global pointer to evolving binary (for BSE Switch Log) + evolvingBinaryStarValid = true; // indicate that the global pointer is now valid (for BSE Switch Log) + + thisBinaryStatus = binary->Evolve(); // evolve the binary + + // announce result + if (!OPTIONS->Quiet()) { // quiet mode? + // no - announce result of evolving the binary + SAY(thisId << ": " << + EVOLUTION_STATUS_LABEL.at(thisBinaryStatus) << ": (" << + STELLAR_TYPE_LABEL.at(binary->Star1InitialType()) << " -> " << + STELLAR_TYPE_LABEL.at(binary->Star1Type()) << ") + (" << + STELLAR_TYPE_LABEL.at(binary->Star2InitialType()) << " -> " << + STELLAR_TYPE_LABEL.at(binary->Star2Type()) << ")" + ); + } + + if (!LOGGING->CloseStandardFile(LOGFILE::BSE_DETAILED_OUTPUT)) { // close detailed output file if necessary + THROW_ERROR_STATIC(ERROR::FILE_NOT_CLOSED); // close failed - this will cause problems later - throw error + } + } + + ERRORS->Clean(); // clean the dynamic error catalog + + if (usingGrid) { // using grid file? + gridLineVariation++; // yes - increment grid line variation number + int optionsStatus = OPTIONS->AdvanceGridLineOptionValues(); // apply next grid file options (ranges/sets) + if (optionsStatus < 0) { // ok? + THROW_ERROR_STATIC(ERROR::ERROR_PROCESSING_GRIDLINE_OPTIONS); // no - throw error + } + else if (optionsStatus == 0) { // end of grid file options variations? + doneGridLine = true; // yes - we're done + } + } + else doneGridLine = true; // not using grid file - done - if (doneGridLine) index = thisId + 1; // increment index + if (doneGridLine) index = thisId + 1; // increment index + } } } - } - delete binary; binary = nullptr; + delete binary; binary = nullptr; - if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // ok? - int optionsStatus = OPTIONS->AdvanceCmdLineOptionValues(); // apply next commandline options (ranges/sets) - if (optionsStatus < 0) { // ok? - evolutionStatus = EVOLUTION_STATUS::ERROR; // no - set status - THROW_ERROR_STATIC(ERROR::ERROR_PROCESSING_CMDLINE_OPTIONS); // throw error - } - else if (optionsStatus == 0) { // end of options variations? - if (usingGrid || OPTIONS->CommandLineGrid() || (!usingGrid && index >= OPTIONS->nObjectsToEvolve())) { // created required number of stars? - evolutionStatus = EVOLUTION_STATUS::DONE; // yes - we're done + if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // ok? + int optionsStatus = OPTIONS->AdvanceCmdLineOptionValues(); // apply next commandline options (ranges/sets) + if (optionsStatus < 0) { // ok? + evolutionStatus = EVOLUTION_STATUS::ERROR; // no - set status + THROW_ERROR_STATIC(ERROR::ERROR_PROCESSING_CMDLINE_OPTIONS); // throw error + } + else if (optionsStatus == 0) { // end of options variations? + if (usingGrid || OPTIONS->CommandLineGrid() || (!usingGrid && index >= OPTIONS->nObjectsToEvolve())) { // created required number of stars? + evolutionStatus = EVOLUTION_STATUS::DONE; // yes - we're done + } } } } - } - // if we catch an error here it happened outside the evolution of a star - meaning that there is a problem // in recording the results of the evolution of the last star evolved, or in setting up the next star to be // evolved - either way we should halt the program here because this could result in undefined behaviour and // results that can't be trusted. So here we just report the error, report what we got up to (how many stars // were evolved before this happened), and terminate the program. - catch (const std::runtime_error& e) { // catch runtime exceptions + catch (const std::runtime_error& e) { // catch runtime exceptions // anything we catch here should not already have been displayed to the user, // so display the error and flag termination (do not rethrow the error) - if (std::string(e.what()) == "FPE") SHOW_ERROR_STATIC(ERROR::FLOATING_POINT_ERROR) // floating-point error - else SHOW_ERROR_STATIC(ERROR::ERROR) // unspecified error - evolutionStatus = EVOLUTION_STATUS::ERROR; // set status + if (std::string(e.what()) == "FPE") SHOW_ERROR_STATIC(ERROR::FLOATING_POINT_ERROR) // floating-point error + else SHOW_ERROR_STATIC(ERROR::ERROR) // unspecified error + evolutionStatus = EVOLUTION_STATUS::ERROR; // set status } - catch (int e) { // catch errors thrown + catch (int e) { // catch errors thrown // anything we catch here should already have been displayed to the user, // so just flag termination - evolutionStatus = EVOLUTION_STATUS::ERROR; // evolution terminated + evolutionStatus = EVOLUTION_STATUS::ERROR; // evolution terminated } - catch (...) { // catchall + catch (...) { // catchall // anything we catch here should not already have been displayed to the user, // so display the error and flag termination (do not rethrow the error) - SHOW_ERROR_STATIC(ERROR::ERROR); // unspecified error - evolutionStatus = EVOLUTION_STATUS::ERROR; // evolution terminated + SHOW_ERROR_STATIC(ERROR::ERROR); // unspecified error + evolutionStatus = EVOLUTION_STATUS::ERROR; // evolution terminated } } @@ -770,11 +787,11 @@ std::tuple EvolveBinaryStars() { // announce result if (!OPTIONS->Quiet()) { - if (evolutionStatus != EVOLUTION_STATUS::CONTINUE) { // shouldn't be... + if (evolutionStatus != EVOLUTION_STATUS::CONTINUE) { // shouldn't be... SAY("\n" << EVOLUTION_STATUS_LABEL.at(evolutionStatus)); } else { - SHOW_WARN_STATIC(ERROR::BINARY_SIMULATION_STOPPED, EVOLUTION_STATUS_LABEL.at(EVOLUTION_STATUS::ERROR)); // show warning + SHOW_WARN_STATIC(ERROR::BINARY_SIMULATION_STOPPED, EVOLUTION_STATUS_LABEL.at(EVOLUTION_STATUS::ERROR)); // show warning } } @@ -782,20 +799,20 @@ std::tuple EvolveBinaryStars() { // don't check result here - let log system handle it (void)LOGGING->CloseAllStandardFiles(); - double cpuSeconds = (clock() - clockStart) / (double) CLOCKS_PER_SEC; // stop CPU timer and calculate seconds + double cpuSeconds = (clock() - clockStart) / (double) CLOCKS_PER_SEC; // stop CPU timer and calculate seconds - auto wallEnd = std::chrono::system_clock::now(); // stop wall timer - std::time_t timeEnd = std::chrono::system_clock::to_time_t(wallEnd); // get end time and date + auto wallEnd = std::chrono::system_clock::now(); // stop wall timer + std::time_t timeEnd = std::chrono::system_clock::to_time_t(wallEnd); // get end time and date SAY("\nEnd generating binaries at " << std::ctime(&timeEnd)); SAY("Clock time = " << cpuSeconds << " CPU seconds"); - std::chrono::duration wallSeconds = wallEnd - wallStart; // elapsed seconds + std::chrono::duration wallSeconds = wallEnd - wallStart; // elapsed seconds - int wallHH = (int)(wallSeconds.count() / 3600.0); // hours - int wallMM = (int)((wallSeconds.count() - ((double)wallHH * 3600.0)) / 60.0); // minutes - int wallSS = (int)(wallSeconds.count() - ((double)wallHH * 3600.0) - ((double)wallMM * 60.0)); // seconds + int wallHH = (int)(wallSeconds.count() / 3600.0); // hours + int wallMM = (int)((wallSeconds.count() - ((double)wallHH * 3600.0)) / 60.0); // minutes + int wallSS = (int)(wallSeconds.count() - ((double)wallHH * 3600.0) - ((double)wallMM * 60.0)); // seconds SAY("Wall time = " << std::setfill('0') << std::setw(4) << wallHH << ":" << std::setfill('0') << std::setw(2) << wallMM << ":" <<