Return to MESA tutorial
run_stars_extra.fto perform calculations not already in MESA
NOTE: Global variables can be found in
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:
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.
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
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
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:
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_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
That's pretty good, but if we were to
./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.
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
./rn this, we will get a model having
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
vals vectors and it will appear in
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.