Run Star Extras

1. Preliminaries
2. Adding History & Profile Output
3. Adding a custom stopping condition
4. Using the "other" hooks
- 4.1. Writing a routine
- 4.2. Instruct MESA to use your routine
5. Overall Program Flow

Please note that this guide was written for r5271 of MESA as part of the 2013 MESA summer school. It is a gradual introduction to run_star_extras.f which tries to introduce some other MESA features and best practices along the way.

1 Preliminaries

1.1 Make a New MESA Working Directory

As you should do each time you want to start a MESA project, make a new copy of the star/work directory.

cp -r $MESA_DIR/star/work rse_tutorial

I've called mine rse_tutorial; you can call yours whatever you want.

If you're familiar with git, there is a git repository which contains snapshots of this MESA working directory at each stage of the tutorial.

git clone https://github.com/jschwab/rse_tutorial

The git tag names do not exactly correspond to the section titles, but hopefully they are intuitively enough named.

1.2 Unpack Configuration Files

Now go ahead an unpack the provided tarball. We have provided a very simple inlist for evolving a 15 Msun star with Z=0.02 and no rotation, no mass loss, etc beginning from the pre-main sequence.

1.3 Activate run_star_extras.f

Navigate to the src directory

cd rse_tutorial/src

and open up run_star_extras.f in your text editor of choice. The stock version of run_star_extras.f is quite boring. It "includes" another file which holds the default set of routines.

include 'standard_run_star_extras.inc'

The routines included in this file are the ones we will want to customize. Because we want these modifications to apply only to this working copy of MESA, and not to MESA as a whole, we want to replace this include statement with the contents of the included file.

Delete the aforementioned include line and insert the contents of $MESA_DIR/include/standard_run_star_extras.inc. (The command to insert the contents of a file in emacs is C-x f <filename> or in vim :r <filename>.)

Before we make any changes, we should check that the code compiles.

cd ..
./mk

If it doesn't compile, make sure that you cleanly inserted the file and removed the include line.

2 Adding History & Profile Output

2.1 Preexisting Output

MESA already knows how to output a tremendous amount of information. The two key file types are history files, which store the value of scalar quantities (e.g. mass, luminosity) at different timesteps and profile files which store the value of spatially varying quantities (e.g. density, pressure) at a single timestep.

The default output is set by the files

$MESA_DIR/star/defaults/history_columns.list
$MESA_DIR/star/defaults/profile_columns.list

In order to customize the output, copy these files to your work directory:

cp $MESA_DIR/star/defaults/history_columns.list .
cp $MESA_DIR/star/defaults/profile_columns.list .

Then, open up history_columns.list or profile_columns.list in a text editor and comment/uncomment any lines to add/remove the columns of interest ('!' is the comment character.)

During our evolution of massive stars we will be interested in the properties of the core. Open up history_columns.list and go to line 256. You should see a section describing properties at the c12_boundary. Uncomment the first four lines.

c12_boundary_mass
   ! c12 boundary is first location going inward from he4 boundary
   ! with c12 abundance <= c12_boundary_limit
c12_boundary_radius
c12_boundary_lgT
c12_boundary_lgRho

As the comment indicates, this boundary marks the first place where the c12 abundance falls below some threshold: this will give us information about the O/Ne (& beyond) core that forms late in the evolution. Take a look at star/controls.defaults to learn more about this variable.

c12_boundary_limit = 1d-4 ! outer boundary of carbon depleted core
   ! c12 boundary is first location going inward from he4 boundary
   ! with c12 abundance <= c12_boundary_limit

Likewise, one can edit your local copy of profile_columns.list to track other quantities. Let's track some more details about the non-nuclear neutrino losses. Jump to line 138 in profile_columns.list and uncomment the following lines.

non_nuc_neu ! non-nuclear-reaction neutrino losses
nonnucneu_plas ! plasmon neutrinos (for collective reactions like gamma_plasmon => nu_e + nubar_e)
nonnucneu_brem ! bremsstrahlung (for reactions like e- + (z,a) => e- + (z,a) + nu + nubar)
nonnucneu_phot ! photon neutrinos (for reactions like e- + gamma => e- + nu_e + nubar_e)
nonnucneu_pair ! pair production (for reactions like e+ + e- => nu_e + nubar_e)

Sometimes one wants to know some details about the interior structure of the star as a function of time. For the most common information – the locations of regions of burning and regions of mixing – MESA can write this to the history file. In other cases, this might require you to output profiles frequently, or to make use of the custom history columns that we'll get to soon.

Open history_columns.list and jump to around line 92, where you'll find a section that looks like this:

!mixing_regions <integer> ! note: this includes regions where the mixing type is no_mixing.

   ! the <integer> is the number of regions to report
   ! there will be 2*<integer> columns for this in the log file, 2 for each region.
   ! the first column for a region gives the mixing type as defined in mlt/public/mlt_def.

   ! the second column for a region gives the m/mstar location of the top of the region
   ! entries for extra columns after the last region in the star will have an invalid mixing_type value of -1.
   ! mstar is the total mass of the star, so these locations range from 0 to 1
   ! all regions are include starting from the center, so the bottom of one region
   ! is the top of the previous one.  since we start at the center, the bottom of the 1st region is 0.

   ! the columns in the log file will have names like 'mix_type_1' and 'mix_qtop_1'

   ! if the star has too many regions to report them all,
   ! the smallest regions will be merged with neighbors for reporting purposes only.

Uncommment the first line and set the <integer> to be 16 which give us plenty of detail.

Find the part of the history_columns.list that deals with burning regions and do the same thing.

But what if you want a history or profile quantity that isn't on the list?

2.2 New History Output

If you already looked at run_star_extras.f you may have noticed the function `how_many_extra_history_columns' and the subroutine `data_for_extra_history_columns'. Adding new data to the history file is as easy as editing these functions.

Perhaps we're interested in quantifying how centrally concentrated the burning is. Let's add columns to our history file that track the Lagrangian mass and physical radius interior to which 90% of the nuclear energy generation takes place.

Let's indicate that we want to add two columns to our history file.

integer function how_many_extra_history_columns(s, id, id_extra)
   type (star_info), pointer :: s
   integer, intent(in) :: id, id_extra
   how_many_extra_history_columns = 2
end function how_many_extra_history_columns

Now let's calculate the information we want to output. This subroutine has access to the star_info pointer, so you can make use of any of the quantities defined in star/private/star_data.inc to compute additional information about the star. Here's a subroutine that calculates what we want to know. There are a few "best practices" shown in this routine that are worth remembering, so read through the code and pay attention to the comments.

subroutine data_for_extra_history_columns(s, id, id_extra, n, names, vals, ierr)

   ! MESA provides routines for taking logarithms which will handle
   ! cases where you pass them zero, etc.  Take advantage of this!
   use num_lib, only : safe_log10

   type (star_info), pointer :: s
   integer, intent(in) :: id, id_extra, n
   character (len=maxlen_history_column_name) :: names(n)
   real(dp) :: vals(n)
   integer, intent(out) :: ierr

   real(dp), parameter :: frac = 0.90
   integer :: i
   real(dp) :: edot, edot_partial

   ! calculate the total nuclear energy release rate by integrating
   ! the specific rate (eps_nuc) over the star.  using the dot_product
   ! intrinsic is a common idiom for calculating integrated quantities.  
   ! note that one needs to explicitly limit the range of the arrays.
   ! NEVER assume that the array size is the same as the number of zones.
   edot = dot_product(s% dm(1:s% nz), s% eps_nuc(1:s% nz))

   ! the center of the star is at i = s% nz and the surface at i = 1 .
   ! so go from the center outward until 90% of the integrated eps_nuc
   ! is enclosed.  exit and then i will contain the desired cell index.
   edot_partial = 0
   do i = s% nz, 1, -1
      edot_partial = edot_partial + s% dm(i) * s% eps_nuc(i)
      if (edot_partial .ge. (frac * edot)) exit
   end do

   ! note: do NOT add these names to history_columns.list
   ! the history_columns.list is only for the built-in log column options.
   ! it must not include the new column names you are adding here.

   ! column 1
   names(1) = "m90"
   vals(1) = s% q(i) * s% star_mass  ! in solar masses

   ! column 2
   names(2) = "log_R90"
   vals(2) = safe_log10(s% R(i) / rsol) ! in solar radii

   ierr = 0
end subroutine data_for_extra_history_columns

2.3 New Profile Output

Adding new output to the profiles, proceeds by close analogy to the history, using the function `how_many_extra_profile_columns' and the subroutine `data_for_extra_profile_columns'.

Here I've just uncommented the stock example.

integer function how_many_extra_profile_columns(s, id, id_extra)
   type (star_info), pointer :: s
   integer, intent(in) :: id, id_extra
   how_many_extra_profile_columns = 1
end function how_many_extra_profile_columns

subroutine data_for_extra_profile_columns(s, id, id_extra, n, nz, names, vals, ierr)
   type (star_info), pointer :: s
   integer, intent(in) :: id, id_extra, n, nz
   character (len=maxlen_profile_column_name) :: names(n)
   real(dp) :: vals(nz,n)
   integer, intent(out) :: ierr
   integer :: k
   ierr = 0

   ! note: do NOT add these names to profile_columns.list
   ! the profile_columns.list is only for the built-in profile column options.
   ! it must not include the new column names you are adding here.

   ! here is an example for adding a profile column
   if (n /= 1) stop 'data_for_extra_profile_columns'
   names(1) = 'beta'
   do k = 1, nz
     vals(k,1) = s% Pgas(k)/s% P(k)
   end do

end subroutine data_for_extra_profile_columns

2.4 Choosing when to output history or profiles

Frequently, you want to output a profile or update the history at a set of milestones along the way. This is especially true for the profiles, when it's not possible (or even desirable) to save and store a profile every few steps.

The star_info structure has a couple of flags – need_to_save_profiles_now and need_to_update_history_now – that let you tell MESA that it is time to output data.

The place to set these is in the function extras_finish_step which will get called at the end of each step.

For our massive star, let's dump a profile at logarithmically-spaced intervals in central density. We will let the user specify the number of divisions per decade. We need to be careful, because the trajectory in rho-T space is not guaranteed to be monotonic.

Here's some code to do just that. Again, read though it as the comments discuss some useful MESA features.

! returns either keep_going or terminate.
! note: cannot request retry or backup; extras_check_model can do that.
integer function extras_finish_step(s, id, id_extra)
   type (star_info), pointer :: s
   integer, intent(in) :: id, id_extra
   integer :: ierr
   integer :: f
   extras_finish_step = keep_going
   call store_extra_info(s)

   ! MESA provides a number of variables that make it easy to get user input.
   ! these are part of the star_info structure and are named
   ! x_character_ctrl, x_integer_ctrl, x_logical_ctrl, and x_ctrl.
   ! by default there are num_x_ctrls, which defaults to 100, of each.
   ! they can be specified in the controls section of your inlist.

   f = s% x_integer_ctrl(1)

   ! MESA also provides a number variables that are useful for implementing
   ! algorithms which require a state. if you just use these variables
   ! restarts, retries, and backups will work without doing anything special.
   ! they are named xtra1 .. xtra30, ixtra1 .. ixtra30, and lxtra1 .. lxtra30.
   ! they are automatically versioned, that is if you set s% xtra1, then
   ! s% xtra1_old will contains the value of s% xtra1 from the previous step
   ! and s% xtra1_older contains the one from two steps ago.

   s% xtra1 = s% log_center_density

   ! this expression will evaluate to true if f times the log center density
   ! has crossed an integer during the last step.  If f = 5, then we will get
   ! output at log center density = {... 1.0, 1.2, 1.4, 1.6, 1.8, 2.0 ... }
   if ((floor(f * s% xtra1_old) - floor(f * s% xtra1) .ne. 0)) then 

      ! save a profile & update the history
      s% need_to_update_history_now = .true.
      s% need_to_save_profiles_now = .true.

      ! by default the priority is 1; you can change that if you'd like
      ! s% save_profiles_model_priority = ?

   endif

end function extras_finish_step

To prevent you from filling up your disk, MESA will only save a limited number of profiles. The default is 100. Here's what Bill has to say in star/defaults/controls.defaults.

max_num_profile_models = 100  ! < 0 means no limit
   ! maximum number of saved profiles
   ! if there's no limit on the number of profiles saved,
   ! you can fill up your disk -- I've done it.
   ! so it's a good idea to set this limit to a reasonable number such as 20 or 30.
   ! once that many have been saved during a run, old ones will be discarded
   ! to make room for new ones.
   ! profiles that were saved for key events are given priority
   ! and aren't removed as long as there
   ! is a lower priority profile that can be discarded instead.

If you want to be sure the profiles that you're triggering in extras_finish_step stick around – perhaps you're making a movie – you should set max_num_profile_models to be greater that the number of profiles you anticipate generating. You might also want to crank up the priority which with they are saved by setting save_profiles_model_priority to be 10. This will prevent MESA from discarding them in lieu of other automatically saved profiles.

3 Adding a custom stopping condition

MESA provides many ways to choose when to terminate a run. If your condition can be expressed as "stop when quantity X rises above or falls below some limit", there's a good chance that you can choose to stop simply by setting a few existing flags. Take a look at the "when to stop" section of star/defaults/controls.defaults, which starts around line 230.

If your condition isn't there or is a more complicated logical combination of conditions, then you will probably need to implement it yourself.

The place to do this in the subroutine extras_check_model. (You can also use this routine to trigger backups or retries.)

Here's a routine that stops when the star's luminosity is dominated by neon burning.

! returns either keep_going, retry, backup, or terminate.
integer function extras_check_model(s, id, id_extra)

   use chem_def, only : i_burn_ne, category_name

   type (star_info), pointer :: s
   integer, intent(in) :: id, id_extra
   integer :: i_burn_max
   extras_check_model = keep_going         

   ! determine the category of maximum burning
   i_burn_max = maxloc(s% L_by_category,1)

   ! stop if the luminosity is dominated by neon burning
   if ( i_burn_max .eq. i_burn_ne) then
      extras_check_model = terminate
      write(*, *) 'terminating because neon burning is dominant'
      return
   end if
end function extras_check_model

4 Using the "other" hooks

MESA provides a way to override most of the physics routines with no need to modify anything more than run_star_extras. There are two main steps needed to take advantage of this functionality. In the following example, we will add controls that allow us to control the various non-nuclear neutrino losses (e.g. plasmon, bremsstrahlung) in our massive star.

4.1 Writing a routine

Navigate to $MESA_DIR/star/other, where you will see a set of files named with the pattern other_*.f, where the wildcard match describes some sort of MESA physics (e.g. other_neu.f, other_wind.f ).

Find the one corresponding to the physics that you want to alter. Open it up and read through the contents. Many of the files have helpful comments and examples.

Note that we do not want to directly edit these files. Instead we want to copy the template routines that these files provide into our working directory copy of run_star_extras.f and then further modify them there. The template routines are named either null_other_* or default_other_*.

For our neutrino example, copy the following subroutine into run_star_extras.f.

subroutine null_other_neu(  &
      id, k, T, log10_T, Rho, log10_Rho, abar, zbar, z2bar, log10_Tlim, flags, &
      loss, sources, ierr)
   use neu_lib, only: neu_get
   use neu_def
   integer, intent(in) :: id ! id for star         
   integer, intent(in) :: k ! cell number or 0 if not for a particular cell         
   real(dp), intent(in) :: T ! temperature
   real(dp), intent(in) :: log10_T ! log10 of temperature
   real(dp), intent(in) :: Rho ! density
   real(dp), intent(in) :: log10_Rho ! log10 of density
   real(dp), intent(in) :: abar ! mean atomic weight
   real(dp), intent(in) :: zbar ! mean charge
   real(dp), intent(in) :: z2bar ! mean charge squared
   real(dp), intent(in) :: log10_Tlim 
   logical, intent(inout) :: flags(num_neu_types) ! true if should include the type of loss
   real(dp), intent(out) :: loss(num_neu_rvs) ! total from all sources
   real(dp), intent(out) :: sources(num_neu_types, num_neu_rvs)
   integer, intent(out) :: ierr
   call neu_get(  &
      T, log10_T, Rho, log10_Rho, abar, zbar, z2bar, log10_Tlim, flags, &
      loss, sources, ierr)
end subroutine null_other_neu

This template routine illustrates the interface and as is, it will produce exactly the same results as the default MESA routine. This is useful because the most frequent sorts of modifications that one wants to make are based on modifying information that MESA already calculates as opposed to a wholesale re-implementation of the physics routines.

In this case, we want to add controls that will allow us to toggle whether specific kinds of non-nuclear neutrino losses are included.

Let's rename the subroutine and add the functionality that we want.

subroutine tutorial_other_neu(  &
     id, k, T, log10_T, Rho, log10_Rho, abar, zbar, z2bar, log10_Tlim, flags, &
     loss, sources, ierr)
   use neu_lib, only: neu_get
   use neu_def
   integer, intent(in) :: id ! id for star         
   integer, intent(in) :: k ! cell number or 0 if not for a particular cell         
   real(dp), intent(in) :: T ! temperature
   real(dp), intent(in) :: log10_T ! log10 of temperature
   real(dp), intent(in) :: Rho ! density
   real(dp), intent(in) :: log10_Rho ! log10 of density
   real(dp), intent(in) :: abar ! mean atomic weight
   real(dp), intent(in) :: zbar ! mean charge
   real(dp), intent(in) :: z2bar ! mean charge squared
   real(dp), intent(in) :: log10_Tlim 
   logical, intent(inout) :: flags(num_neu_types) ! true if should include the type of loss
   real(dp), intent(out) :: loss(num_neu_rvs) ! total from all sources
   real(dp), intent(out) :: sources(num_neu_types, num_neu_rvs)
   integer, intent(out) :: ierr

   ! before we can use controls associated with the star we need to get access 
   type (star_info), pointer :: s
   call star_ptr(id, s, ierr)
   if (ierr /= 0) then ! OOPS
      return
   end if

   ! separately control whether each type of neutrino loss is included
   flags(pair_neu_type) = s% x_logical_ctrl(1)
   flags(plas_neu_type) = s% x_logical_ctrl(2)
   flags(phot_neu_type) = s% x_logical_ctrl(3)
   flags(brem_neu_type) = s% x_logical_ctrl(4)
   flags(reco_neu_type) = s% x_logical_ctrl(5)

   ! the is the normal routine that MESA provides
   call neu_get(  &
       T, log10_T, Rho, log10_Rho, abar, zbar, z2bar, log10_Tlim, flags, &
       loss, sources, ierr)

end subroutine tutorial_other_neu

4.2 Instruct MESA to use your routine

There are two things that you must to in order to have MESA execute your other_* routine. Failure to do both of these is the most common problem people encounter when using the other_* hooks.

First, edit the controls section of your inlist to set the appropriate use_other_* flag to .true. . In our example, this means adding the line

use_other_neu = .true.

Second, edit the extras_controls routine in run_star_extras.f to point the other_neu at the routine you want to be executed.

subroutine extras_controls(s, ierr)
   type (star_info), pointer :: s
   integer, intent(out) :: ierr
   ierr = 0

   ! this is the place to set any procedure pointers you want to change
   ! e.g., other_wind, other_mixing, other_energy  (see star_data.inc)
   s% other_neu => tutorial_other_neu

end subroutine extras_controls

Now, recompile your working directory

./mk

5 Overall Program Flow

This is an overview of when the functions defined in run_star_extras.f get called during MESA execution. This is quite high-level. If you're interested in how MESA goes about solving the stellar structure equations – that is, what happens in star_evolve_step – you should read the appendix of Paxton et al. (2013).

If this is your first time dealing with run_star_extras, it's probably best to skip over this section for now and come back to it later. It'll make more sense once you've see some of the functions in action.

If you're at the other extreme and interested in the deep guts of MESA, you can look at the real version of run1_star in star/job/run_star_support.f.

Here's the flow in fortran-ish pseudocode

subroutine run1_star

   call star_setup(id, inlist_fname, ierr)

   call extras_controls(s, ierr)

   call do_star_job_controls_before(id, s, restart, ierr)
   call do_load1_star(id, s, restart, 'restart_photo', ierr)
   call do_star_job_controls_after(id, s, restart, ierr)

   id_extra = extras_startup(s, id, restart, ierr)

   evolve_loop: do while(continue_evolve_loop) ! evolve one step per loop

      step_loop: do ! may need to repeat this loop

         result = star_evolve_step(id, first_try, just_did_backup)
         if (result == keep_going) result = star_check_model(id)
         if (result == keep_going) result = extras_check_model(s, id, id_extra)
         if (result == keep_going) exit step_loop

      end do step_loop

      if (result == keep_going) then

         ! if you have data that needs to be saved and restored for restarts,
         ! save it in s% extra_iwork and s% extra_work
         ! before calling star_finish_step
         result = extras_finish_step(s, id, id_extra)
         if (result == terminate) exit evolve_loop

         result = star_finish_step(id, id_extra, .false., &
            how_many_extra_profile_columns, data_for_extra_profile_columns, &
            how_many_extra_history_columns, data_for_extra_history_columns, ierr)
         if (result /= keep_going) exit evolve_loop

      else if (result == terminate) then

         if (result_reason == result_reason_normal) then

            call save_profile(id, id_extra, &
               how_many_extra_profile_columns, data_for_extra_profile_columns, &
               3, ierr)
            s% need_to_save_profiles_now = .false.
            s% need_to_update_history_now = .true.
            result = star_finish_step(id, id_extra, save_photo_when_terminate, &
               how_many_extra_profile_columns, data_for_extra_profile_columns, &
               how_many_extra_history_columns, data_for_extra_history_columns, ierr)
            call do_saves( &
               id, id_extra, s, &
               how_many_extra_history_columns, &
               data_for_extra_history_columns, &
               how_many_extra_profile_columns, &
               data_for_extra_profile_columns)
         end if
         exit evolve_loop
      end if

      call do_saves( &
         id, id_extra, s, &
         how_many_extra_history_columns, &
         data_for_extra_history_columns, &
         how_many_extra_profile_columns, &
         data_for_extra_profile_columns)

   end do evolve_loop

   call extras_after_evolve(s, id, id_extra, ierr)

   call star_after_evolve(id, ierr)

end subroutine run1_star