Skip to main content

Extending MESA Skip to project menu

Return to MESA tutorial

Learning outcomes

  • Learn how to write your own run_stars_extra.f to perform calculations not already in MESA
  • Pass arguments through the inlist to custom routines
  • Add custom defined variables to history log file

NOTE: Global variables can be found in $MESA_DIR/const/public/const_def.f90.

The default run_stars_extras.f file, which can be found in /src from a MESA work directory, looks as follows:

      module run_star_extras

      use star_lib
      use star_def
      use const_def
      use crlibm_lib

      implicit none

      ! these routines are called by the standard run_star check_model
      contains

      include 'standard_run_star_extras.inc'

      end module run_star_extras

As can be seen, all of the heavy lifting is found in standard_run_star_extras.inc. This file can be found in $MESA_DIR/include/standard_run_star_extras.inc. To create your own run_stars_extra.f routines, you will need to first paste the contents of that file here.

The execution of MESA proceeds as follows:


(Courtesy Josiah Schwab, MESA Summer School 2016)

All of the boxes to the right are places in run_stars_extra.f where control can be taken. The gray take_step box is where all of the heavy machinery of MESA takes place.

An example extension

We will now modify MESA to calculate the blackbody temperature of the Earth. We will further make MESA stop if it exceeds a certain value. That temperature can be calculated with

$$ \frac{T_\oplus}{T_\odot} = \sqrt{\frac{R_\odot}{2\text{AU}}}. $$

In order to insert a stopping condition, we must modify extras_finish_step. We need to declare the variable we want to use (Tearth) and then return a command to terminate if our condition is met. We furthermore want the user to be able to specify the stopping value. We can do that as follows:

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

         real(dp) :: Tearth

         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_finish_step = keep_going
         call store_extra_info(s)

         ! calculate blackbody temperature of earth
         Tearth = s% Teff * sqrt(s% photosphere_r * Rsun / (2.0 * AU))
         write(*,*) "Tearth =", Tearth

         ! stop if it exceeds 300 K
         if (Tearth > s% x_ctrl(1)) extras_finish_step = terminate

         ! to save a profile,
            ! s% need_to_save_profiles_now = .true.
         ! to update the star log,
            ! s% need_to_update_history_now = .true.

         ! see extras_check_model for information about custom termination codes
         ! by default, indicate where (in the code) MESA terminated
         if (extras_finish_step == terminate) s% termination_code = t_extras_finish_step
      end function extras_finish_step

NOTE 1: Structure variables in Fortran can be accessed with %, e.g. s% Teff accesses the effective temperature of variable s. This is analogous to s.Teff in C-like languages.

NOTE 2: Constants (i.e. values defined in $MESA_DIR/const) do not need s% to be accessed; e.g. Rsun can be accessed directly.

NOTE 3: The comments before the header of the function call says the possible return values: keep_going and terminate.

NOTE 4: x_ctrl is a vector of real numbers that the user can set in the inlist. In this case, we have designated x_ctrl(1) to denote the maximum blackbody temperature of the Earth. One can use as many as desired. Also available are x_integer_ctrl and x_logical_ctrl which are vectors for integers and booleans, respectively.

Now we can modify our inlist to have, for example,

x_ctrl(1) = 300 ! K; Tearth cut-off value

inside of the &controls section.

That's pretty good, but if we were to ./mk and ./rn this, we would obtain a model whose temperature exceeds the value we set. Instead, it might be better to iterate until we get one exactly on the value that we want.

Iteratively finding a specified value

For this, we can modify extras_check_model and retry until we get the exact value -- or at least one within a user-specified tolerance:

      ! returns either keep_going, retry, backup, or terminate.
      integer function extras_check_model(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s

         real(dp) :: Tearth, Tearth_cutoff, Tearth_tol

         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_check_model = keep_going
         if (.false. .and. s% star_mass_h1 < 0.35d0) then
            ! stop when star hydrogen mass drops to specified level
            extras_check_model = terminate
            write(*, *) 'have reached desired hydrogen mass'
            return
         end if

         Tearth_cutoff = s% x_ctrl(1)
         Tearth_tol = s% x_ctrl(2)
         Tearth = s% Teff * sqrt(s% photosphere_r * Rsun / (2.0 * AU))

         if (Tearth < Tearth_cutoff - Tearth_tol) then
             extras_check_model = keep_going
             write(*, *) 'Tearth too low, keep going'
         end if
         if (Tearth > Tearth_cutoff + Tearth_tol) then
             extras_check_model = retry
             s% dt = s% dt / 10
             write(*, *) 'Tearth exceeds limit, backing up'
         end if
         if (Tearth >= Tearth_cutoff - Tearth_tol .and. &
             Tearth <= Tearth_cutoff + Tearth_tol) then
             extras_check_model = terminate
             write(*, *) 'Bullseye!'
         end if

         ! if you want to check multiple conditions, it can be useful
         ! to set a different termination code depending on which
         ! condition was triggered.  MESA provides 9 customizeable
         ! termination codes, named t_xtra1 .. t_xtra9.  You can
         ! customize the messages that will be printed upon exit by
         ! setting the corresponding termination_code_str value.
         ! termination_code_str(t_xtra1) = 'my termination condition'

         ! by default, indicate where (in the code) MESA terminated
         if (extras_check_model == terminate) s% termination_code = t_extras_check_model
      end function extras_check_model

Now we can modify our inlist to have, for example,

x_ctrl(1) = 300 ! K; Tearth cut-off value
x_ctrl(2) = 1d-6 ! K; Tearth tolerance

inside of the &controls section. When we ./mk and ./rn this, we will get a model having $T_\oplus$ within a millionth of the specified value of 300 K.

Adding a history column

Finally, we want these values to be reported in the LOGS/history.data file. We must first tell MESA that we want another column:

      integer function how_many_extra_history_columns(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_history_columns = 1
      end function how_many_extra_history_columns

Now we need only to add the calculation to the names and vals vectors and it will appear in LOGS/history.data:

      subroutine data_for_extra_history_columns(id, id_extra, n, names, vals, ierr)
         integer, intent(in) :: id, id_extra, n
         character (len=maxlen_history_column_name) :: names(n)
         real(dp) :: vals(n)
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return

         !note: do NOT add the extras 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.

         names(1) = "Tearth"
         vals(1) = s% Teff * sqrt(s% photosphere_r * Rsun / (2.0 * AU))

      end subroutine data_for_extra_history_columns

Of course we must ./mk again, and then if we ./rn we will obtain these numbers at each time step.


Figure: Blackbody temperature evolution over the history of planet Earth
Creative Commons License
Unless otherwise specified, original content from Earl Bellinger's Ph.D. is licensed under a
Creative Commons Attribution 4.0 International License.