Skip to content

Commit

Permalink
Merge pull request #632 from MESAHub/EbF/fix_rates_detailed_balance
Browse files Browse the repository at this point in the history 
fix rates detailed balance
  • Loading branch information
@Debraheem
Debraheem committed May 6, 2024
2 parents 550a702 + ec9283c commit 7cff4c9
Show file tree
Hide file tree
Showing 4 changed files with 294 additions and 263 deletions.
16 changes: 16 additions & 0 deletions docs/source/changelog.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,22 @@ New Features
Bug Fixes
---------

Rates
~~~~~

There has been a bug present in the rates module due to the incorrect
phase space factors for reverse reaction rates involving greater than 2 reactants or
products. This bug resulted in inconsistent equillibrium compositions when the network
evolves into nuclear statistical equillibrium (NSE), at temperatures exceeding 4 GK.
This bug effects users who evolve models into NSE using large reaction networks. This
includes evolving massive stars to core-collapse. Smaller networks such as the ``approx21``
networks are less affected. We strongly recommend that users update to the latest MESA release.

See `gh-575 <https://github.com/MESAHub/mesa/issues/575>`_

``max_allowed_nz``
~~~~~~~~~~~~~~~~~~

MESA no longer produces a segmentation fault if it tries to increase
the number of cells beyond ``max_allowed_nz``.

Expand Down
61 changes: 31 additions & 30 deletions net/test/test_output
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,27 +105,27 @@
log temp 9.600000000
log rho 6.000000000

log(eps_nuc) 24.178055700
log(eps_nuc) 24.177999295

eps_nuc 0.150680031E+25
eps_nuc 0.150660462E+25

d_lneps_dlnT -0.251306070
d_lneps_dlnRho 1.008425906
d_lneps_dlnT -0.253424018
d_lneps_dlnRho 1.008418067


energy generation by category

category log rate (erg/g/sec)

fe_co_ni 24.178013 0.150665E+25
photo 20.161509 0.145047E+21
si_alpha 18.180628 0.151575E+19
s_alpha 18.053619 0.113141E+19
ar_alpha 17.551212 0.355805E+18
mg_alpha 17.035498 0.108517E+18
ti_alpha 16.444832 0.278504E+17
ca_alpha 16.021612 0.105102E+17
cr_alpha 15.767744 0.585792E+16
fe_co_ni 24.177953 0.150644E+25
photo 20.161465 0.145032E+21
s_alpha 18.844132 0.698444E+19
si_alpha 18.653887 0.450699E+19
ar_alpha 18.392651 0.246974E+19
mg_alpha 17.802202 0.634164E+18
ti_alpha 17.300311 0.199669E+18
ca_alpha 16.259105 0.181596E+17
cr_alpha 15.768197 0.586404E+16
o_alpha 15.265751 0.184396E+16
ne_alpha 14.736951 0.545696E+15
c_alpha 14.073434 0.118422E+15
Expand All @@ -147,27 +147,27 @@
log temp 9.600000000
log rho 6.000000000

log(eps_nuc) 24.178055700
log(eps_nuc) 24.177999295

eps_nuc 0.150680031E+25
eps_nuc 0.150660462E+25

d_lneps_dlnT -0.251306070
d_lneps_dlnRho 1.008425906
d_lneps_dlnT -0.253424018
d_lneps_dlnRho 1.008418067


energy generation by category

category log rate (erg/g/sec)

fe_co_ni 24.178013 0.150665E+25
photo 20.161509 0.145047E+21
si_alpha 18.180628 0.151575E+19
s_alpha 18.053619 0.113141E+19
ar_alpha 17.551212 0.355805E+18
mg_alpha 17.035498 0.108517E+18
ti_alpha 16.444832 0.278504E+17
ca_alpha 16.021612 0.105102E+17
cr_alpha 15.767744 0.585792E+16
fe_co_ni 24.177953 0.150644E+25
photo 20.161465 0.145032E+21
s_alpha 18.844132 0.698444E+19
si_alpha 18.653887 0.450699E+19
ar_alpha 18.392651 0.246974E+19
mg_alpha 17.802202 0.634164E+18
ti_alpha 17.300311 0.199669E+18
ca_alpha 16.259105 0.181596E+17
cr_alpha 15.768197 0.586404E+16
o_alpha 15.265751 0.184396E+16
ne_alpha 14.736951 0.545696E+15
c_alpha 14.073434 0.118422E+15
Expand Down Expand Up @@ -203,10 +203,11 @@
test_one_zone_burn_const_P
number of species 21
large final abundances 1.0000000000000000D-02
o16 1 6.9572247281582056D-01
si28 2 2.0457753929379313D-01
s32 3 4.4284801800242163D-02
mg24 4 3.7924075365987756D-02
o16 1 7.1798810044269901D-01
si28 2 1.6142410416068431D-01
s32 3 6.0439378324441878D-02
ar36 4 2.9692711889885499D-02
ca40 5 1.6560528697491497D-02

xsum 1.0000000000000000D+00

Expand Down
48 changes: 31 additions & 17 deletions rates/private/reaclib_support.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,29 +217,43 @@ subroutine compute_rev_ratio(rates,winvn)

! log(prefactor of reverse_ratio)
tmp = product(mass(1:Ni))/product(mass(Ni+1:Nt))
rates% inverse_coefficients(1,i) = pow(tmp,1.5d0)*(product(g(1:Ni))/product(g(Ni+1:Nt)))
rates% inverse_coefficients(1,i) = pow(tmp,1.5d0*(Ni-No))*(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

! This should be 0 for non-photo-disintegration reverse rates and 1 for photos's in the reverse channel
if (No==1) then
rates% inverse_exp(i) = 1
if(rates% inverse_coefficients(2,i)<0) then
! negative values denote endothermic photodisintegrations
! We want rate_photo/rate_forward
rates% inverse_coefficients(1,i) = rates% inverse_coefficients(1,i) * fac
else
! positive values denote exothermic photodisintegrations
! We divide by fac and invert the T^3/2 as we want to compute
! rate_reverse/rate_photo
rates% inverse_coefficients(1,i) = rates% inverse_coefficients(1,i) / fac
rates% inverse_exp(i) = -1 ! We us this term in a log() expression
end if
else
rates% inverse_exp(i) = 0
! see equation 21 in Reichert et al. 2023 (https://doi.org/10.3847/1538-4365/acf033)
! delta_reactants/delta_products is handled in net.
! fac == (mu * kb * T / (2 * pi * hbar^2))^3/2 term with out a T^3/2.
! fac shows up as fac^(Ni-No) in rates% inverse_coefficients(1,i)
! so rates% inverse_coefficients(1,i) contains terms for
! fac^(n) == fac^(Ni-No), where n = Ni - No.

! 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

! Ni-No should be 0 for non-photo-disintegration reverse rates
! and >=1 for photos's in the reverse channel
if (Ni-No .ne. 0) then
! whether endothermic or exothermic,
! Ni-No handles the sign of Q from rates% inverse_coefficients(2,i)
rates% inverse_coefficients(1,i) = rates% inverse_coefficients(1,i) * pow(fac, Ni-No)
! negative values of Q denote endothermic photodisintegrations
! We multiply by fac^|Ni-No| as we want to compute
! rate_photo/rate_forward

! positive values of Q denote exothermic photodisintegrations
! We DIVIDE by fac^|Ni-No| as we want to compute
! rate_reverse/rate_photo
else ! Ni - No = 0
rates% inverse_exp(i) = 0 ! ensure this is 0.
rates% inverse_coefficients(1,i) = rates% inverse_coefficients(1,i) ! no point calling pow(fac,0)
end if
rates% inverse_coefficients(1,i) = log(rates% inverse_coefficients(1,i))

Expand Down

0 comments on commit 7cff4c9

Please sign in to comment.