[mesa-users] Change reaction rate of the triple alpha process

Bill Paxton paxton at kitp.ucsb.edu
Mon Aug 1 11:53:34 EDT 2016


I 2nd Frank's comment about changing the code -- it should not be necessary.  There are simpler ways.

The scaling factors Frank mentioned are temperature independent.  If you want something the depends on temperature, then you should create a table of (T,rate) pairs for mesa to use instead of the standard.  If you want to try that, there are lots of folks (well, a few anyway) who can help -- just send another message to the list.

-B


On Aug 1, 2016, at 8:39 AM, Francis Timmes wrote:

> changing code is probably not what you want; 
> rate multipliers can be done via the inlist:
> 
> http://mesa.sourceforge.net/star_job_defaults.html#num_special_rate_factors
> 
> the overall scaling factors
> 
> http://mesa.sourceforge.net/controls_defaults.html#eps_nuc_factor
> http://mesa.sourceforge.net/controls_defaults.html#dxdt_nuc_factor
> 
> may also help you reach your goal.
> 
> fxt
> 
> 
> 
> 
>> On Aug 1, 2016, at 3:15 AM, Byeong-Chan Park <pbc1918 at gmail.com> wrote:
>> 
>> Dears.
>> 
>> I need to change the rate to investigate energy generation process in a star. Hence I need to find the code that means nuclear reaction rates as a function of temperature. I checked 'ratelib.f90', 'rates_def_mic.f90', and 'rates_def.f90' in the mesa/rates/privates and mesa/rates/public directory.
>> 
>> However I can't distinguish what code says the reaction rate formular of the triple alpha process in these 3 files. I think that the next in the mesa/rates/private/ratelib.f90 includes reaction formular;
>> 
>> 
>> ! r3a, triple alpha
>> 
>> 
>>      subroutine rate_tripalf_jina(tf, temp, fr, rr)
>>         type (T_Factors), pointer :: tf
>>         real(dp), intent(in) :: temp
>>         real(dp), intent(out) :: fr, rr
>>      real(dp) :: fr1, rr1
>>      include 'formats.dek'
>> 
>> !       he4  he4  he4  c12                  fy05r     7.27500d+00          
>>         call jina_reaclib_3_1(ihe4, ihe4, ihe4, ic12, tf, fr, rr, 'rate_tripalf_jina')
>> 
>>         return 
>>         call rate_tripalf_reaclib(tf, temp, fr1, rr1)
>> 
>>         write(*,1) 'fr', fr
>>         write(*,1) 'fr1', fr1
>>         write(*,1) 'rr', rr
>>         write(*,1) 'rr1', rr1
>>         write(*,*)
>>         stop 'rate_tripalf_jina' 
>>      end subroutine rate_tripalf_jina
>> 
>> 
>> 
>>      subroutine rate_tripalf_reaclib(tf, temp, fr, rr)
>>      !      HE4(2A,G)C12    reaclib JINA - Fynbo et al. 2005 Nature 433, 136-139
>>      type (T_Factors), pointer :: tf
>>      real(dp), intent(in) :: temp
>>      real(dp), intent(out) :: fr, rr
>>      real(dp) :: fr1, fr2, fr3, rev
>>      rr = 0
>>      if (temp < 1d7) then
>>         fr = 0; return
>>      end if
>>      call do_reaclib(tf,   &
>>                  -9.710520d-01,  0.000000d+00, -3.706000d+01,  2.934930d+01,                        &
>>                  -1.155070d+02, -1.000000d+01, -1.333330d+00,                                     &
>>                  fr1)
>>      call do_reaclib(tf,   &
>>                  -2.435050d+01, -4.126560d+00, -1.349000d+01,  2.142590d+01,                        &
>>                  -1.347690d+00,  8.798160d-02, -1.316530d+01,                                     &
>>                  fr2)
>>      call do_reaclib(tf,   &
>>                  -1.178840d+01, -1.024460d+00, -2.357000d+01,  2.048860d+01,                        &
>>                  -1.298820d+01, -2.000000d+01, -2.166670d+00,                                     &
>>                  fr3)
>>         fr = fr1 + fr2 + fr3         
>>         ! use the fxt reverse rate term
>>         rev    = 2.00d+20*(tf% t93)*exp_cr(-84.424*(tf% t9i))
>>         rr = fr * rev         
>>      end subroutine rate_tripalf_reaclib
>> 
>> 
>> 
>>      subroutine rate_tripalf_nacre(tf, temp, fr, rr)
>>         type (T_Factors), pointer :: tf
>>         real(dp), intent(in) :: temp
>>         real(dp), intent(out) :: fr, rr
>>         real(dp) :: r2abe, rbeac, bb, term, rev
>>         ! he4(a, g)be8
>>         ! a0 T9i23 exp_cr(-a1 T9i13 - (T9*a2)^2) 
>>         !     * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95) 
>>         ! + c0 T9i32 exp_cr(-c1/T9)
>>         ! + d0 T9i32 exp_cr(-d1/T9)
>>         ! + e0 T9^e1 exp_cr(-e2/T9)
>> 
>>         if (tf% t9 < lowT9_cutoff) then
>>            fr = 0; rr = 0; return
>>         end if 
>> 
>>         call rnacre(tf,  &
>>               2.43d9, 13.490d0, 1d0/0.15d0, & ! a0, a1, a2
>>               74.5d0, 0d0, 0d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
>>               6.09d5, 1.054d0, &  ! c0, c1
>>               0d0, 0d0, & ! d0, d1
>>               0d0, 0d0, 0d0, & ! e0, e1, e2
>>               r2abe)         
>>        ! be8(a, g)c12
>>         ! a0 T9i23 exp_cr(-a1 T9i13 - (T9*a2)^2) 
>>         !     * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95) 
>>         ! + c0 T9i32 exp_cr(-c1/T9)
>>         ! + d0 T9i32 exp_cr(-d1/T9)
>>         ! + e0 T9^e1 exp_cr(-e2/T9)
>>         call rnacre(tf,  &
>>               2.76d7, 23.570d0, 1d0/0.4d0, & ! a0, a1, a2
>>               5.47d0, 326d0, 0d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
>>               130.7d0, 3.338d0, &  ! c0, c1
>>               2.51d4, 20.307d0, & ! d0, d1
>>               0d0, 0d0, 0d0, & ! e0, e1, e2
>>               rbeac)               
>>         if (tf% T9 <= 0.03d0) then
>>            bb    = 3.07d-16*(1 - 29.1d0*(tf% T9) + 1308d0*(tf% T92))
>>            if (bb < 0) then
>>               bb = 0
>>            end if      
>>         else
>>            bb    = 3.44d-16*(1 +     0.0158d0*pow_cr(tf% T9,-0.65d0))
>>         end if      
>>         term    = r2abe * rbeac * bb
>>         call rnacre_rev(tf, &  ! a0 T932 exp_cr(-a1/T9)
>>            2.003d20, 84.415d0, &  ! a0, a1
>>            rev)     
>>         fr    = term
>>         rr    = rev * term
>>      end subroutine rate_tripalf_nacre
>> 
>> 
>>      subroutine rate_tripalf_fxt(tf, temp, fr, rr)
>>         type (T_Factors), pointer :: tf
>>         real(dp), intent(in) :: temp
>>         real(dp), intent(out) :: fr, rr
>> 
>>      real(dp) term, dtermdt, rev, drevdt, r2abe, dr2abedt, rbeac,  &
>>                       drbeacdt, aa, daa, bb, dbb, cc, dcc, dd, ddd, ee, dee,  &
>>                       ff, dff, xx, dxx, yy, dyy, zz, dzz, uu, vv, f1, df1, rc28,  &
>>                       q1, q2
>>      parameter        (rc28   = 0.1d0,  &
>>                        q1     = 1.0d0/0.009604d0,  &
>>                        q2     = 1.0d0/0.055225d0) 
>> 
>> 
>>         if (tf% t9 < lowT9_cutoff) then
>>            fr = 0; rr = 0; return
>>         end if 
>> 
>> ! this is a(a, g)be8
>>         aa    = 7.40d+05 * (tf% t9i32) * exp_cr(-1.0663*(tf% t9i)) 
>> 
>>         bb    = 4.164d+09 * (tf% t9i23) * exp_cr(-13.49*(tf% t9i13) - (tf% t92)*q1)
>> 
>>         cc    = 1.0d0 + 0.031*(tf% t913) + 8.009*(tf% t923) + 1.732*(tf% t9)   &
>>              + 49.883*(tf% t943) + 27.426*(tf% t953)
>> 
>>         r2abe    = aa + bb * cc
>> 
>> ! this is be8(a, g)c12
>>         dd    = 130.0d0 * (tf% t9i32) * exp_cr(-3.3364*(tf% t9i))  
>> 
>>         ee    = 2.510d+07 * (tf% t9i23) * exp_cr(-23.57*(tf% t9i13) - (tf% t92)*q2)
>> 
>>         ff    = 1.0d0 + 0.018*(tf% t913) + 5.249*(tf% t923) + 0.650*(tf% t9) +  &
>>              19.176*(tf% t943) + 6.034*(tf% t953)
>> 
>>         rbeac    = dd + ee * ff
>> 
>> ! a factor
>>         xx    = rc28 * 1.35d-07 * (tf% t9i32) * exp_cr(-24.811*(tf% t9i))
>> 
>> 
>> ! high temperature rate
>>         if ((tf% t9).gt.0.08) then
>>          term    = 2.90d-16 * r2abe * rbeac + xx
>> 
>> 
>> ! low temperature rate
>>         else
>>          uu   = 0.8d0*exp_cr(-pow_cr(0.025*(tf% t9i),3.263d0)) 
>>          yy   = 0.2d0 + uu  ! fixes a typo in Frank's original
>>          vv   = 4.0d0*exp_cr(-pow_cr((tf% t9)/0.025d0,9.227d0)) 
>>          zz   = 1.0d0 + vv
>>          aa   = 1.0d0/zz
>>          f1   = 0.01 + yy * aa  ! fixes a typo in Frank's original
>>          term = 2.90d-16 * r2abe * rbeac * f1 +  xx 
>>         end if
>> 
>>         rev    = 2.00d+20*(tf% t93)*exp_cr(-84.424d0*(tf% t9i))
>> 
>> !      term    = 1.2d0 * term
>> !      dtermdt = 1.2d0 * term
>> 
>>      fr    = term
>> 
>>      rr    = rev * term
>> 
>>      end subroutine rate_tripalf_fxt
>> 
>> 
>> but I don't know what code that I should change. In addition, it makes me confused that form of the 'subroutine rate_tripalf_jina' is different with the others.
>> 
>> I need your help.
>> 
>> Best regards
>> Park.
>> ------------------------------------------------------------------------------
>> _______________________________________________
>> mesa-users mailing list
>> mesa-users at lists.sourceforge.net
>> https://lists.sourceforge.net/lists/listinfo/mesa-users
> 
> 
> ------------------------------------------------------------------------------
> _______________________________________________
> mesa-users mailing list
> mesa-users at lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/mesa-users





More information about the Mesa-users mailing list