[mesa-users] Questions about making stellar isochrone and changing triple alpha reaction rates

[Bill Paxton]

[Dear other readers, 
If you aren't a rates junky, you might want to skip most of this, but at least skim the paragraph including "This is silly!!!!". ;p]

On Oct 18, 2013, at 10:32 AM, 이현택 wrote:
> Hi
> 
> I have one more questions about this topic.
>  I examined what you suggested and now I clearly understand the structure of the source code.
>  But it seems that the additional temprature(T8) and reaction rate pair can only have some discrete values. 
>  The continuous curve that I want can be easily calculated by using subroutine reaclib, which is declared in ~/mesa/rates/private/ratelib.f. 


Hi,

Good question -- to answer it we need to go even deeper into the mesa/rates module.

As you mention, the implementation of mesa/rates includes routines to evaluate reaclib rates as well as other analytic approximations such as CF88 and NACRE.  Rather than execute those expressions repeatedly, we automatically make a table for interpolation and cache that for future use.   You'll find those files in mesa/data/rates_data/cache.  The file names start with the name of the reaction and end with ".bin" indicating that they are binary files.  Since there can be multiple rates for a given reaction, the file name also includes a number indicating the choice of rate.  For example, in my current rates cache I find 3 files for triple alpha with the names r_he4_he4_he4_to_c12_1.bin, r_he4_he4_he4_to_c12_2.bin, and r_he4_he4_he4_to_c12_4.bin.  (Apparently, I haven't use the rate option "3"!)

The cached rates have values for many temperatures covering the full range supported.   The number of points and the range are determined by values in the rates_def module (see mesa/rates/public/rates_def.f).   You can change the defaults, but that's probably not necessary.
      
      ! for tabular evaluation of the raw reaction rates
      real(dp) :: rattab_thi != 10.30102999566398d0 ! log10(highest temp = 2e10)
      real(dp) :: rattab_tlo != 5.30102999566398d0 ! log10(lowest temp = 2e5)
      integer :: rattab_points_per_decade = 2000

With 2000 points per decade and 5 decades, we have 10,000 points for doing interpolation; that seems to be adequate.  The interpolation is done using Steffen's piecewise monotonic cubic routine "interp_pm" from mesa/interp_1d.

That's how it works for rates given by fitting expressions from sources such as reaclib or NACRE: we sample the analytic expressions for 10,000 temperatures and then do monotonicity preserving cubic interpolation using those values.  

Now to get to your question: for rates given as a table (i.e., for reactions included in mesa/data/rates_data/rate_tables/rate_list.txt), the mesa/rates routines get the values for the 10,000 temperatures to go in the cache by interpolating the given table.   Then those values are cached and used exactly the same way as built-in rates.   If you create a table with all of the 10,000 temperatures that will go into the cache, then you'll completely determine the rates used by mesa.

This is a general solution to the need for specifying reaction rates.   Hopefully someday all rates will be provided in this form.  Think about it --- currently we start from a set of points for rates at different temperatures, then someone fits those points to a complex formula and publishes the best fit coefficients.  Then we take those coefficients and recreate a table of points.  This is silly!!!!   Once in the long-past old days, computers were much more limited and storing large tables of rates wasn't an option.  Times have changed and now we are happy with large tables, but because of historical inertia, we still distribute rates as approximate fits.  In a rational world, we'd directly distribute the points without introducing the inaccuracy of fitting to a preset formula.

Until that time, I'm going to ask you to take your new reaclib coefficients and write a small routine to create a table by evaluating the analytic approximation at a large number of temperatures.    (and to answer your other question, yes, T8 = T*10^-8)

Can that work for you?

Bill



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20131018/192b35e5/attachment.html>