Skip to content

Commit 7108ae4

Browse files
rjfarmerrjfarmer
authored
Merge pull request #494 from MESAHub/rf/missing_rates
Various fixes for missing reaction rates as well as handling the special rate factors
2 parents 5f984d2 + 1406c51 commit 7108ae4

File tree

10 files changed

+85
-14
lines changed

10 files changed

+85
-14
lines changed

‎data/rates_data/reactions.list‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -355,6 +355,7 @@ r_h2_pg_he3 1 h1 + 1 h2 => 1 he3
355355
r_he3_ag_be7 1 he3 + 1 he4 => 1 be7 2 he3 he4 pp he4(he3,g)be7
356356
r_he3_be7_to_h1_h1_he4_he4 1 he3 + 1 be7 => 2 h1 + 2 he4 2 he3 be7 pp he3 + be7 => 2 p + 2 he4
357357
r_li7_pa_he4 1 li7 + 1 h1 => 2 he4 2 h1 li7 pp li7(p,a)he4
358+
r_he4_ap_li7 1 he4 + 1 he4 => 1 li7 + 1 h1 2 he4 he4 other he4(a,p)li7
358359
r_n13_ap_o16 1 n13 + 1 he4 => 1 h1 + 1 o16 2 he4 n13 cno n13(a,p)o16
359360
r_n13_gp_c12 1 n13 => 1 h1 + 1 c12 1 photo n13(g,p)c12
360361
r_n13_pg_o14 1 n13 + 1 h1 => 1 o14 2 h1 n13 cno n13(p,g)o14

‎docs/source/known_bugs.rst‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,23 @@ a newer version of MESA. Note this list is NOT comprehensive, users should check
1111
issue but it may not be complete.
1212

1313

14+
r22.11.1
15+
========
16+
17+
Rates
18+
=====
19+
20+
There is a bug in the rate selection code that certain endothermic weak reactions are not added to the nuclear network. These are
21+
r_be10_wk-minus_b10, r_ni66_wk-minus_cu66, and r_h3_wk-minus_he3. Other weak reactions with heavier parents may also be affected.
22+
23+
A separate issue also meant we are missing the rate r_he4_ap_li7 as the reverse rate of r_li7_pa_he4.
24+
25+
Both issues will effect previous versions of MESA as well.
26+
27+
Both issues have been fixed in the git main branch.
28+
29+
30+
1431
r22.05.1
1532
========
1633

‎rates/private/rates_names.f90‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -87,6 +87,7 @@ subroutine set_reaction_names
8787
reaction_Name(ir_h1_h1_he4_to_he3_he3) = 'r_h1_h1_he4_to_he3_he3'
8888
reaction_Name(ir_he4_he4_he4_to_c12) = 'r_he4_he4_he4_to_c12'
8989
reaction_Name(ir_li7_pa_he4) = 'r_li7_pa_he4'
90+
reaction_Name(ir_he4_ap_li7) = 'r_he4_ap_li7'
9091
reaction_Name(ir_mg24_ag_si28) = 'r_mg24_ag_si28'
9192
reaction_Name(ir_mg24_ap_al27) = 'r_mg24_ap_al27'
9293
reaction_Name(ir_mg24_ga_ne20) = 'r_mg24_ga_ne20'

‎rates/private/rates_reverses.f90‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -168,6 +168,9 @@ subroutine set_reaction_reverses
168168
reverse_reaction_id(irfe52neut_to_fe54) = irfe54g_to_fe52
169169
reverse_reaction_id(irfe54g_to_fe52) = irfe52neut_to_fe54
170170

171+
reverse_reaction_id(ir_li7_pa_he4) = ir_he4_ap_li7
172+
reverse_reaction_id(ir_he4_ap_li7) = ir_li7_pa_he4
173+
171174
end subroutine set_reaction_reverses
172175

173176

‎rates/private/raw_rates.f90‎

Lines changed: 19 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -224,6 +224,9 @@ subroutine set_raw_rate(ir, temp, tf, raw_rate, ierr)
224224
case(ir_li7_pa_he4) ! li7(p, a)he4
225225
call do1(rate_li7pa_nacre)
226226

227+
case(ir_he4_ap_li7) ! li7(p, a)he4
228+
call do1_reverse(rate_li7pa_nacre)
229+
227230
case(ir_n13_gp_c12) ! n13(g, p)c12
228231
call do1_reverse(rate_c12pg_nacre)
229232

@@ -906,7 +909,10 @@ subroutine eval_table_reverse(ir, rir, tf, temp, fr, rr, ierr)
906909
dlambda_dlnT = 0d0
907910

908911
call do_reaclib_indices_for_reaction(reaction_name(rir), reaclib_rates, lo, hi, ierr)
909-
if(ierr/=0) return
912+
if(ierr/=0) then
913+
write(*,*) "Error: Could not get reaclib index for rate ",rir, ' ',trim(reaction_name(rir))
914+
return
915+
end if
910916

911917
inv_lambda = 0d0
912918
dinv_lambda_dlnT = 0d0
@@ -945,14 +951,21 @@ subroutine get_interp_table(f_name, nT8s, T8s_out, f1_out, ierr)
945951
ierr = 0
946952
vec => vec_ary
947953

948-
rate_file = trim(rates_table_dir) // '/' // trim(f_name)
949-
!write(*,*) 'load table ' // trim(rate_file)
954+
955+
! Look for the file based on its name first
956+
957+
rate_file = trim(f_name)
950958

951959
open(newunit=iounit,file=trim(rate_file),action='read',status='old',iostat=ierr)
952960
if (ierr /= 0) then
953-
write(*,*) 'ERROR: cannot open rate info file ' // trim(rate_file)
954-
!return
955-
call mesa_error(__FILE__,__LINE__)
961+
! Look in rates_table_dir
962+
rate_file = trim(rates_table_dir) // '/' // trim(f_name)
963+
open(newunit=iounit,file=trim(rate_file),action='read',status='old',iostat=ierr)
964+
if (ierr /= 0) then
965+
write(*,*) 'ERROR: cannot open rate info file ' // trim(rate_file)
966+
!return
967+
call mesa_error(__FILE__,__LINE__)
968+
end if
956969
end if
957970

958971
do ! read until reach line starting with an integer (nT8s)

‎rates/private/reaclib_input.f90‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -248,7 +248,7 @@ subroutine extract_rates_from_reaclib(reaclib,nreaclib,nuclides,rates,set,use_we
248248
integer :: i,j,l,count,loc_count,nt,indx,cat, &
249249
weaklib_count,chapter,num_in,num_out,num_skip_for_weaklib,num_from_reaclib, &
250250
max_lhs_Z, min_Z, i1, i2, cid_in, cid_out
251-
logical :: include_this_rate, already_included_from_weaklib, found_it
251+
logical :: include_this_rate, already_included_from_weaklib, found_it, is_weak
252252
integer, dimension(max_species_per_reaction) :: pspecies
253253
character(len=max_id_length) :: handle
254254
character(len=iso_name_length) :: name_i, name_j, name1, name2, label
@@ -329,12 +329,16 @@ subroutine extract_rates_from_reaclib(reaclib,nreaclib,nuclides,rates,set,use_we
329329

330330
end do loop_over_nuclides
331331

332+
is_weak = (use_weaklib .and. num_in == 1 .and. num_out == 1)
333+
332334
! only include forward rates
333335
! Define the reverse rate as being the endothermic reaction, always
334336
! Some photo disintegrations can be exothermic, so dont trust what REACLIB calls a reverse rate
335-
if(include_this_rate) then
337+
if(include_this_rate .and. .not. is_weak) then
336338
if(sum(nuclides% binding_energy(pspecies(1:num_in))) - &
337-
sum(nuclides% binding_energy(pspecies(num_in+1:Nt))) > 0 ) cycle loop_over_rates
339+
sum(nuclides% binding_energy(pspecies(num_in+1:Nt))) > 0 ) then
340+
cycle loop_over_rates
341+
end if
338342
end if
339343
if (include_this_rate) then
340344
if (use_weaklib .and. num_in == 1 .and. num_out == 1) then

‎rates/private/reaclib_support.f90‎

Lines changed: 16 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -636,8 +636,22 @@ end function two_one
636636

637637
logical function two_two()
638638
two_two = .true.
639-
if (in1 == 0 .or. in2 == 0 .or. out1 == 0 .or. out2 == 0 .or. &
640-
in1 == in2 .or. out1 == out2) return
639+
640+
if (in1 == 0 .or. in2 == 0 .or. out1 == 0 .or. out2 == 0) return
641+
642+
! Special case r_li7_pa_he4, this must come first otherwise the out1==out2
643+
! check will label this rate as a _to_ reaction instead of an _ap reaction
644+
if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
645+
nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
646+
nuclides% Z(out2) == 2 .and. nuclides% N(out2) == 2 .and. &
647+
nuclides% Z(in2) == 3 .and. nuclides% N(in2) == 4) then
648+
handle = trim(handle) // '_pa'
649+
two_two = .false.
650+
return
651+
end if
652+
653+
if(in1==in2 .or. out1==out2) return
654+
641655
two_two = .false.
642656
if (nuclides% Z(in1) == 2 .and. nuclides% N(in1) == 2 .and. &
643657
nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &

‎rates/public/rates_def.f90‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -362,7 +362,8 @@ end subroutine interpolate_weak_table
362362
integer, parameter :: ir_h1_h1_he4_to_he3_he3 = ir_he3_he3_to_h1_h1_he4+1
363363
integer, parameter :: ir_he4_he4_he4_to_c12 = ir_h1_h1_he4_to_he3_he3+1
364364
integer, parameter :: ir_li7_pa_he4 = ir_he4_he4_he4_to_c12+1
365-
integer, parameter :: ir_mg24_ag_si28 = ir_li7_pa_he4+1
365+
integer, parameter :: ir_he4_ap_li7 = ir_li7_pa_he4+1
366+
integer, parameter :: ir_mg24_ag_si28 = ir_he4_ap_li7+1
366367
integer, parameter :: ir_mg24_ap_al27 = ir_mg24_ag_si28+1
367368
integer, parameter :: ir_mg24_ga_ne20 = ir_mg24_ap_al27+1
368369
integer, parameter :: ir_n13_ap_o16 = ir_mg24_ga_ne20+1

‎rates/test/test_output‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -417,6 +417,7 @@
417417
r_h1_h1_he4_to_he3_he3 0.0000000000000000D+00
418418
r_he4_he4_he4_to_c12 0.0000000000000000D+00
419419
r_li7_pa_he4 1.1481972998876799D-26
420+
r_he4_ap_li7 0.0000000000000000D+00
420421
r_mg24_ag_si28 0.0000000000000000D+00
421422
r_mg24_ap_al27 0.0000000000000000D+00
422423
r_mg24_ga_ne20 0.0000000000000000D+00
@@ -727,6 +728,7 @@
727728
r_h1_h1_he4_to_he3_he3 0.0000000000000000D+00
728729
r_he4_he4_he4_to_c12 0.0000000000000000D+00
729730
r_li7_pa_he4 2.8536849590579373D-15
731+
r_he4_ap_li7 0.0000000000000000D+00
730732
r_mg24_ag_si28 0.0000000000000000D+00
731733
r_mg24_ap_al27 0.0000000000000000D+00
732734
r_mg24_ga_ne20 0.0000000000000000D+00
@@ -1037,6 +1039,7 @@
10371039
r_h1_h1_he4_to_he3_he3 0.0000000000000000D+00
10381040
r_he4_he4_he4_to_c12 6.0425276466603135D-71
10391041
r_li7_pa_he4 1.3036488512272843D-07
1042+
r_he4_ap_li7 0.0000000000000000D+00
10401043
r_mg24_ag_si28 0.0000000000000000D+00
10411044
r_mg24_ap_al27 0.0000000000000000D+00
10421045
r_mg24_ga_ne20 0.0000000000000000D+00
@@ -1347,6 +1350,7 @@
13471350
r_h1_h1_he4_to_he3_he3 0.0000000000000000D+00
13481351
r_he4_he4_he4_to_c12 2.1189548985111926D-46
13491352
r_li7_pa_he4 1.7169556787367157D-02
1353+
r_he4_ap_li7 0.0000000000000000D+00
13501354
r_mg24_ag_si28 3.4094273819191269-162
13511355
r_mg24_ap_al27 9.1263133445433935-274
13521356
r_mg24_ga_ne20 0.0000000000000000D+00
@@ -1657,6 +1661,7 @@
16571661
r_h1_h1_he4_to_he3_he3 0.0000000000000000D+00
16581662
r_he4_he4_he4_to_c12 2.0403192412842943D-24
16591663
r_li7_pa_he4 4.3342837593883942D+01
1664+
r_he4_ap_li7 0.0000000000000000D+00
16601665
r_mg24_ag_si28 1.7640238996806857D-49
16611666
r_mg24_ap_al27 2.1146812954972443D-91
16621667
r_mg24_ga_ne20 0.0000000000000000D+00
@@ -1967,6 +1972,7 @@
19671972
r_h1_h1_he4_to_he3_he3 1.1772013589467652-212
19681973
r_he4_he4_he4_to_c12 7.9590873420121651D-13
19691974
r_li7_pa_he4 7.8177297732306224D+03
1975+
r_he4_ap_li7 2.2679834766347775-273
19701976
r_mg24_ag_si28 2.6469697259346481D-15
19711977
r_mg24_ap_al27 8.0770504592160291D-29
19721978
r_mg24_ga_ne20 3.1486672973522443-150
@@ -2277,6 +2283,7 @@
22772283
r_h1_h1_he4_to_he3_he3 1.2168141674471324D-69
22782284
r_he4_he4_he4_to_c12 3.4041066123320119D-10
22792285
r_li7_pa_he4 2.2984785945795634D+05
2286+
r_he4_ap_li7 4.0535520962769382D-82
22802287
r_mg24_ag_si28 8.0870253378865792D-04
22812288
r_mg24_ap_al27 3.1845928303516034D-06
22822289
r_mg24_ga_ne20 1.0831093286103609D-39
@@ -2587,6 +2594,7 @@
25872594
r_h1_h1_he4_to_he3_he3 2.9036145678610530D-23
25882595
r_he4_he4_he4_to_c12 2.2461964456422845D-10
25892596
r_li7_pa_he4 2.1991028257065220D+06
2597+
r_he4_ap_li7 1.3073398877915164D-20
25902598
r_mg24_ag_si28 1.4662488667811619D+01
25912599
r_mg24_ap_al27 9.0039206583785949D+02
25922600
r_mg24_ga_ne20 1.6781713325994766D-02
@@ -2897,6 +2905,7 @@
28972905
r_h1_h1_he4_to_he3_he3 1.2872091694058283D-07
28982906
r_he4_he4_he4_to_c12 2.1320831532048911D-11
28992907
r_li7_pa_he4 1.6571622519109823D+07
2908+
r_he4_ap_li7 4.4349794705459109D+00
29002909
r_mg24_ag_si28 6.9791363732282798D+02
29012910
r_mg24_ap_al27 3.0549748551998041D+06
29022911
r_mg24_ga_ne20 3.2362267060968735D+10

‎star/job/run_star_support.f90‎

Lines changed: 10 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1582,6 +1582,14 @@ subroutine set_rate_factors(id, ierr)
15821582

15831583
s% rate_factors(:) = 1
15841584
if (s% job% num_special_rate_factors <= 0) return
1585+
1586+
! Dont error if we are changing net
1587+
if ((s% job% change_initial_net .or. s% job% change_net) .and. &
1588+
trim(s% job% new_net_name)/=trim(s% net_name)) then
1589+
!write(*,*) "Not changing special rates untill net change"
1590+
return
1591+
end if
1592+
15851593

15861594
call get_net_reaction_table_ptr(s% net_handle, net_reaction_ptr, ierr)
15871595
if (ierr /= 0) return
@@ -1592,14 +1600,14 @@ subroutine set_rate_factors(id, ierr)
15921600
j = 0
15931601
if (ir > 0) j = net_reaction_ptr(ir)
15941602
if (j <= 0) then
1595-
write(*,2) 'Failed to find reaction_for_special_factor ' // &
1603+
write(*,*) 'Failed to find reaction_for_special_factor ' // &
15961604
trim(s% job% reaction_for_special_factor(i)), &
15971605
j, s% job% special_rate_factor(i)
15981606
error = .true.
15991607
cycle
16001608
end if
16011609
s% rate_factors(j) = s% job% special_rate_factor(i)
1602-
write(*,2) 'set special rate factor for ' // &
1610+
write(*,*) 'set special rate factor for ' // &
16031611
trim(s% job% reaction_for_special_factor(i)), &
16041612
j, s% job% special_rate_factor(i)
16051613
end do

0 commit comments

Comments
 (0)