8000 Fix cosmic by MathieuVenet · Pull Request #694 · COSMIC-PopSynth/COSMIC · GitHub
[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Skip to content

Fix cosmic #694

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Apr 1, 2025
Merged

Fix cosmic #694

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5d22072
Update evolv2.f
MathieuVenet Mar 4, 2025
7a2a2b1
Update evolv2.f
MathieuVenet Mar 12, 2025
b39d935
Add files via upload
MathieuVenet Apr 1, 2025
fed98d2
Update test_evolve.py
MathieuVenet Apr 1, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 18 additions & 9 deletions src/cosmic/src/evolv2.f
View file
Open in desktop
8000
Original file line number Diff line number Diff line change
Expand Up @@ -1824,7 +1824,12 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
*
7 km0 = dtm0*1.0d+03/tb
if(km0.lt.tiny) km0 = 0.5d0
*

* Check for collision at periastron for a stable RLOF
pd = sep*(1.d0 - ecc)
if(pd.lt.(rad(1)+rad(2))) goto 130

*
* Force co-rotation of primary and orbit to ensure that the tides do not
* lead to unstable Roche (not currently used).
*
Expand Down Expand Up @@ -2449,29 +2454,27 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
*
evolve_type = 8.0

age = tphys - epoch(1)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@MathieuVenet can you explain to me what this line is doing and why it was changed, as well as the line that reads in aj(1) in hrdiag?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Originally this line is supposed to define the current evolution time as an input for the hrdiag function. However, when the comenv function is called a new variable aj(1) (and (aj(2)) is defined. This variable is a new evolution time much shorter that the one before the comenv function is called.

When defining the current evolution time, and not this new value, as an input for the hrdiag function, an issue arises.

In such cases, luminosity is computed by the lgbtf function, which relies on this age parameter. If the age is too large, the function encounters a negative value raised to a specific power, which Fortran cannot process, resulting in a NaN. This NaN then propagates through key parameters such as the stellar radius, system separation, and other essential quantities, leading to erroneous results.

To resolve this issue, I removed the redundant redefinition of the evolution time, as it is unnecessary for the time step calculation. Instead, I modified hrdiag to use the time variable aj as defined by the comenv function, ensuring numerical stability and accurate computations.


mc = massc(1)
rc = radc(1)
CALL star(kstar(1),mass0(1),mass(1),tm,tn,tscls,lums,GB,zpars)
CALL hrdiag(mass0(1),age,mass(1),tm,tn,tscls,lums,GB,zpars,
CALL hrdiag(mass0(1),aj(1),mass(1),tm,tn,tscls,lums,GB,zpars,
& rm,lum,kstar(1),mc,rc,me,re,k2,bhspin(1),1)
aj(1) = age
epoch(1) = tphys - age

rad(1) = rm
lumin(1) = lum
massc(1) = mc
radc(1) = rc
menv(1) = me
renv(1) = re

age = tphys - epoch(2)

mc = massc(2)
rc = radc(2)
CALL star(kstar(2),mass0(2),mass(2),tm,tn,tscls,lums,GB,zpars)
CALL hrdiag(mass0(2),age,mass(2),tm,tn,tscls,lums,GB,zpars,
CALL hrdiag(mass0(2),aj(2),mass(2),tm,tn,tscls,lums,GB,zpars,
& rm,lum,kstar(2),mc,rc,me,re,k2,bhspin(2),2)
aj(2) = age
epoch(2) = tphys - age

rad(2) = rm
lumin(2) = lum
massc(2) = mc
Expand Down Expand Up @@ -2520,6 +2523,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
& formation(2),binstate,mergertype,'bpp')
*
epoch(j1) = tphys - aj(j1)
com = .false.
if(coel)then
com = .true.
goto 135
Expand Down Expand Up @@ -3556,6 +3560,11 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars)
CALL hrdiag(m0,age,mt,tm,tn,tscls,lums,GB,zpars,
& rm,lum,kw,mc,rc,me,re,k2,bhspin(k),k)
pd = sep*(1.d0 - ecc)
if(pd.lt.(rad(1)+rad(2))) goto 130



*
* Check for a supernova and correct the semi-major axis if so.
*
Expand Down
Binary file modified src/cosmic/tests/data/kick_initial_conditions.h5
View file
Open in desktop
Binary file not shown.
2 changes: 1 addition & 1 deletion src/cosmic/tests/test_evolve.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,4 +105,4 @@ def test_ejection_velocity_pfahl(self):
EvolvedBinaryBPP, EvolvedBinaryBCM, initCond, kick_info = Evolve.evolve(
initialbinarytable=KICK_INITC)

self.assertAlmostEqual(kick_info['vsys_2_total'].iloc[0], 17.322114, places=5)
self.assertAlmostEqual(kick_info['vsys_2_total'].iloc[0], 482.346136, places=5)
0