Skip to content

Bug in reverse rates mass exponent#975

Merged
Debraheem merged 8 commits into
release/r26.4.1from
bugfix/reverse_rate_mass_r26.4.1
Apr 23, 2026
Merged

Bug in reverse rates mass exponent#975
Debraheem merged 8 commits into
release/r26.4.1from
bugfix/reverse_rate_mass_r26.4.1

Conversation

@Debraheem
Copy link
Copy Markdown
Member

This critical bugfix addresses a mistake in the mass exponent term in the reverse reaction rates, raised by @haakoan here: #974.

I have added an additional unit test which tests changes to these reverse ratios in the future. I believe this is important enough to warrant being merged into the r26.4.1 release.

@Debraheem Debraheem requested a review from fxt44 as a code owner April 21, 2026 22:02
@Debraheem Debraheem added bug Something isn't working release-blocker Things that should be fixed before the next release rates-net Rates and net modules labels Apr 21, 2026
@Debraheem Debraheem linked an issue Apr 21, 2026 that may be closed by this pull request
@Debraheem Debraheem changed the title fix mass exponent, add inverse coefficients to rates testing Bug in reverse rates mass exponent Apr 21, 2026
@Debraheem Debraheem requested a review from rhdtownsend as a code owner April 21, 2026 22:24
Comment thread net/test/test_output
Comment thread docs/source/changelog.rst Outdated
Comment thread rates/private/reaclib_support.f90 Outdated
@Debraheem Debraheem requested review from aurimontem and mathren April 22, 2026 01:18
@Debraheem
Copy link
Copy Markdown
Member Author

Any further reviews are appreciated as I would like this pr merged and tested asap, so we can correct it in this release.

Comment thread docs/source/changelog.rst Outdated
@matthiasfabry
Copy link
Copy Markdown
Contributor

I am by no means an expert in nuclear physics, but from your comments it seems both formulations exist in the literature? Is Reichert's Eq. 21 simply a typo? If so, I'd clarify that in a code comment so someone doesn't mistakenly reintroduce it by seeing the discrepancy again.

@Debraheem
Copy link
Copy Markdown
Member Author

I am by no means an expert in nuclear physics, but from your comments it seems both formulations exist in the literature? Is Reichert's Eq. 21 simply a typo? If so, I'd clarify that in a code comment so someone doesn't mistakenly reintroduce it by seeing the discrepancy again.

yes, agreed, i believe it is a typo in their manuscript but not in their code. I made comments on this in reaclib_support.f90:

! log(prefactor of reverse_ratio)
tmp = product(mass(1:Ni))/product(mass(Ni+1:Nt))
! The reactant/product mass ratio always enters with a fixed 3/2 power.
! This follows standard detailed-balance derivations such as
! Fowler et al. 1967 and Smith Clark et al. 2023 Appendix B.
rates% inverse_coefficients(1,i) = pow(tmp,1.5d0)*(product(g(1:Ni))/product(g(Ni+1:Nt)))
! -Q/(kB*10**9)
sum1 = sum(winvn% binding_energy(ps(1:Ni)))
sum2 = sum(winvn% binding_energy(ps(Ni+1:Nt)))
rates% inverse_coefficients(2,i) = (sum1-sum2)*conv/kB/1d9
! Reichert et al. 2023 Eq. 21 is the source for the fac^(Ni-No) term,
! but the fixed 3/2 mass exponent above follows Fowler et al. 1967 and
! Smith Clark et al. 2023 Appendix B. The current WinNet implementation
! also uses the fixed 3/2 mass exponent.
! delta_reactants/delta_products is handled in net.
! fac == (mu * kb * T / (2 * pi * hbar^2))^3/2 term without a T^3/2.
! The n-dependent phase-space factor enters as fac^(Ni-No), where
! n = Ni - No. The mass ratio above remains a separate fixed 3/2 term.
! The T makes its way back into our expression inside
! the subroutine compute_some_inverse_lambdas, in reaclib_eval.f90.
! It appears in log form as 1.5d0*rates% inverse_exp(i)*lnT9, where,
! rates% inverse_exp(i) = Ni-No, so,
! 1.5d0*rates% inverse_exp(i)*lnT9 == ln(T^(3n/2)), where n = Ni-No.
rates% inverse_exp(i) = Ni - No ! We use this term in a log() expression in reaclib_eval

@Debraheem
Copy link
Copy Markdown
Member Author

I compared the original and corrected reverse-rate implementations with BBQ hydrostatic one-zone burns using six MESA networks: approx21_plus_co56, mesa_80, mesa_151, mesa_206, mesa_330, and mesa_495. The sweep covered three physical test setups, expanded into five concrete burn cases: pure he4 at 0.1 GK and 0.2 GK with logRho = 4.0, each evolved for 100 Myr; a 50/50 by-mass c12/o16 mixture at 1 GK with logRho = 6.0, evolved for 1000 yr; and Z-split compositions at 3 GK and 6 GK with logRho = 8.0, each evolved for 1e5 s, where each network was initialized with equal total mass fraction in each nuclear charge Z present in that network, and that Z fraction was then split equally among the isotopes with that Z. These were fixed-duration burns, not runs stopped by an equilibrium criterion. From these runs I generated three plot sets in plots/bbq_sweep: by_temperature, which compares the final isotope mass fractions across networks for a single source state (new or old) and includes both grouped bar plots and multi-panel pie summaries of the top 15 final isotopes; compare_labels, which compares the final isotope distributions for new versus old for the same network and burn case; and eps_nuc, which compares the corresponding eps_nuc(t) histories versus age for new and old for the same network and burn case.

Because i made hundreds of plots I share a zip, but please if others will download and look at the results, and provide comment it would be appreciated!
bbq_sweep.zip

when looking at the plots: old = r26.04.1 and new = r26.04.1 + mass_exponent_fix

A coupe notable things, It seems low temperatures are relatively unaffected in abundance or eps_nuc, but higher temperatures near NSE conditions at 3GK and 6GK are most affected. This is good news in my opinion, as it means most of H/He/C/O burning is unaffected T<1GK. Only when reverse rates become important does this correction come into play.

As for the test_suite, I've returned the old flame speed and width targets to the conductive_flame test case, as the intial bugfix for the reverse rates which introduced this max_exponent bug is what caused us to change the flame targets in the first place.

@Debraheem
Copy link
Copy Markdown
Member Author

495_comparison.pdf
approx_comparison.pdf
I share pdfs of two of the comparisons for approx21 and a 495 isotope network, under the conditions described above.

Expanded explanation of the incorrect mass exponent issue and its relevance in advanced burning stages. Added reference to related issue gh-974 for more context.
Copy link
Copy Markdown
Contributor

@mathren mathren left a comment

Choose a reason for hiding this comment

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

Great catch @haakoan and @Debraheem!

The one-zone-burner tests you shared are pretty convincing that this is a small, but noticeable effect for late evolution. In the interest of not blocking the release, I am approving this PR and I think it could be merged as is.

I would still run in parallel tests with full stars (maybe one mass that goes to core-collapse) which may show the effect is even smaller (not pushing too hard the models, if they fail because of any other reason, then this should become a science project on its own, and interested users should always run their own tests anyways).

Based on the bbq tests, I don't expect tests with full stars to surprise us or introduce astrophysically significant issues.

@Debraheem
Copy link
Copy Markdown
Member Author

testing seems fine excluding the conductive_flame test. https://testhub.mesastar.org/bugfix%2Freverse_rate_mass_r26.4.1/commits/5eb8b12

I ran an optional test locally and found that ppisn now fails on my arm machine on this branch. I'm not surprised given how sensitive it is in its current form. I will try to address this, and run another test.

@haakoan
Copy link
Copy Markdown

haakoan commented Apr 23, 2026

Great catch @haakoan and @Debraheem!

The one-zone-burner tests you shared are pretty convincing that this is a small, but noticeable effect for late evolution. In the interest of not blocking the release, I am approving this PR and I think it could be merged as is.

I would still run in parallel tests with full stars (maybe one mass that goes to core-collapse) which may show the effect is even smaller (not pushing too hard the models, if they fail because of any other reason, then this should become a science project on its own, and interested users should always run their own tests anyways).

Based on the bbq tests, I don't expect tests with full stars to surprise us or introduce astrophysically significant issues.

I can add that in my testing against FLASH, I found differences around of 10-50%, in terms of mass fractions. However, it is not a clean test. We ran MESA models until core-collapse, went back 10-15 minutes and then started 3D FLASH simulations from the MESA profiles. Thus, differences in the convective activity can account for some of that difference. The rates used for approx21 in FLASH are also older than the MESA rates and I do see some differences in the raw rates too. To complicate the issue further, the discrepancies are model dependent.

It is also not clear to me how the approx21 results would transfer to larger networks.

@Debraheem Debraheem merged commit c8c2026 into release/r26.4.1 Apr 23, 2026
4 of 5 checks passed
@Debraheem
Copy link
Copy Markdown
Member Author

@VincentVanlaer VincentVanlaer deleted the bugfix/reverse_rate_mass_r26.4.1 branch April 23, 2026 21:12
Debraheem added a commit that referenced this pull request Apr 24, 2026
* fix mass exponent, add inverse coefficients to rates testing

* [ci skip] update changelog and known bugs

* [ci skip] update docs

* [ci skip] Fix Haakan's last name

* return conductive flame properties to old values

* [ci skip] update changelog

* [ci skip] Clarify known bug details in reverse rates calculation

Expanded explanation of the incorrect mass exponent issue and its relevance in advanced burning stages. Added reference to related issue gh-974 for more context.

* [ci optional] tweak ppisn, correct mlt_vc interpolation in split_merge
VincentVanlaer pushed a commit that referenced this pull request Apr 29, 2026
* fix mass exponent, add inverse coefficients to rates testing

* [ci skip] update changelog and known bugs

* [ci skip] update docs

* [ci skip] Fix Haakan's last name

* return conductive flame properties to old values

* [ci skip] update changelog

* [ci skip] Clarify known bug details in reverse rates calculation

Expanded explanation of the incorrect mass exponent issue and its relevance in advanced burning stages. Added reference to related issue gh-974 for more context.

* [ci optional] tweak ppisn, correct mlt_vc interpolation in split_merge
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bug Something isn't working rates-net Rates and net modules release-blocker Things that should be fixed before the next release

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Question about compute_rev_ratio

5 participants