MESA Summer School 2017: Lecture 1
Table of Contents
Part 0: Overview
This guide was written as part of the 2017 MESA summer school. It is
an introduction to MESA/star and MESA/binary, with a particular focus
on using run_star_extras.f and run_binary_extras.f. It assumes you
are using r9793 of MESA.
If you're new to Fortran, we prepared a short document with some examples. Don't let yourself get hung up by the Fortran; quickly ask your classmates and the TAs for help!
There is a version of this document available with solutions. The git repository hosting this document contains the full source code used in each task, which you can see by looking at the appropriately named tag (i.e. part2-task2).
Part 1: Running and controlling MESA
If you've used MESA before much of this should be familiar.
Part 1a: Getting started
Each time you want to start a MESA project, you should make a new copy
of the star/work directory.
cp -r $MESA_DIR/star/work lecture1
In this case, we have prepared and provided a work directory for you. Download, unpack, and enter this work directory.
unzip lecture1.zip cd lecture1
Task 1: Compile and run the provided work directory
This directory evolves a solar mass star from the middle of the main sequence to hydrogen exhaustion. Confirm that you can compile and run it. A window with a few plots should appear.
Answer
./clean ./mk ./rn
Part 1b: Using inlists
MESA/star has three inlist sections. Each section contains the options for a different aspect of MESA.
- star_job
- options for the program that evolves the star
- controls
- options for the MESA star module
- pgstar
- options for on-screen plotting
The distinction between star_job and controls can be a little
subtle. We won't discuss pgstar in this lecture, but Frank will
later this morning.
star_job contains options that answer questions like:
- how should MESA obtain the initial model?
- are there any changes MESA should make to the initial model?
- what microphysics data should MESA read?
- where should MESA store its output?
controls contains options that answer questions like:
- when should MESA stop evolving the model?
- which angular momentum transport processes should MESA consider?
- what numerical tolerances should MESA's solvers use?
MESA's many inlist options are documented in the files
$MESA_DIR/star/defaults/star_job.defaults$MESA_DIR/star/defaults/controls.defaults$MESA_DIR/star/defaults/pgstar.defaults
They are roughly sorted into groups of related options. When you're searching for an option, see if it seems to match any of the section headings and then look there first. If that fails, try searching for some key words.
Note that inlists can point to other inlists. This can be useful for
keeping things organized. In the lecture1 directory, inlist
points MESA to inlist_project and inlist_pgstar.
Task 2: Read the documentation
Use the documentation to learn about the stopping condition that we
are using (xa_central_lower_limit and
xa_central_lower_limit_species). Look near this option to see some
of the other mass fraction based stopping conditions that are
available.
Answer
Look at the version of controls.defaults on the web or included in MESA. Note that the main webpage has the documentation for the latest release. You should consult the documentation appropriate to your version of MESA. Either use the online documentation archive or should consult the defaults files included in your version of MESA.
The documentation for xa_central_lower_limit and xa_central_lower_limit_species says:
Lower limits on central mass fractions. Stop when central abundance drops below this limit. Can have up to num_xa_central_limits of these (see star_def.inc for value). xa_central_lower_limit_species contains an isotope name as defined in chem_def.f. xa_central_lower_limit contains the lower limit value.
Nearby you will find similar controls with upper/lower limits on average and surface mass fractions.
Part 1c: Controlling 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 contents of MESA's output files is not directly controlled via inlists. 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, you 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.)
You can use run_star_extras.f to define your own history and/or
profile columns. This capability is covered on the MESA website; we
will not cover it today.
Task 3: Add some output
Look at LOGS/history.data and LOGS/profile1.data to see what
information is included by default. In our later exercises, we will
be setting the variable extra_heat, which is defined at each cell in
the star. Add this quantity to the output. Run MESA and confirm that
the column you want is there (its value should be zero).
Answer
Uncomment the following lines in profile_columns.list
extra_heat
and then run MESA
./rn
Part 2: Using run_star_extras.f
To activate run_star_extras.f, navigate to the lecture1/src
directory and open 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 defined in the included 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 i <filename>, vim :r
<filename>, or you can just copy and paste.)
Before we make any changes, we should check that the code compiles.
cd .. ./mk
If it doesn't compile, double check that you cleanly inserted the file and removed the include line.
The two most important things that one needs to know in order to use
run_star_extras.f effectively are (1) the control flow of a MESA run
and (2) the contents of the star_info structure.
The different run_star_extras.f routines get called at different
points during MESA execution. Here is a high-level overview of a MESA
run, written in Fortran-ish pseudocode.
subroutine run1_star(...) ! star is initialized here ! before evolve loop calls: ! extras_controls ! extras_startup call before_evolve_loop(...) ! evolve one step per loop evolve_loop: do while(continue_evolve_loop) call before_step_loop(...) step_loop: do ! may need to repeat this loop if (stop_is_requested(s)) then continue_evolve_loop = .false. result = terminate exit end if result = star_evolve_step(...) if (result == keep_going) result = star_check_model(...) if (result == keep_going) result = extras_check_model(...) if (result == keep_going) result = star_pick_next_timestep(...) if (result == keep_going) exit step_loop ! redo, retry, or backup must be done inside the step_loop if (result == redo) then result = star_prepare_to_redo(...) end if if (result == retry) then result = star_prepare_to_retry(...) end if if (result == backup) then result = star_do1_backup(...) just_did_backup = .true. else just_did_backup = .false. end if if (result == terminate) then continue_evolve_loop = .false. exit step_loop end if end do step_loop ! once we get here, the only options are keep_going or terminate. ! after_step_loop calls: ! extras_finish_step call after_step_loop(...) if (result /= keep_going) then exit evolve_loop end if ! write out data ! ! do_saves calls: ! how_many_extra_history_columns ! data_for_extra_history_columns ! how_many_extra_profile_columns ! data_for_extra_profile_columns call do_saves(...) end do evolve_loop ! after_evolve_loop calls: ! extras_after_evolve call after_evolve_loop(...) end subroutine run1_star
In even more distilled terms, here is a flowchart summarizing this.
The heart of MESA is the grey "take step" box, which contains all of the machinery by which MESA evaluates and solves the equations of stellar structure.
The star_info structure contains all the information about the star
that is being evolved. By convention, the variable name s is used
throughout run_star_extras.f to refer to this structure. In
Fortran, the percent (%) operator is used to access the components of
the structure. (So you can read s% x = 3 in the same way that you
would read s.x = 3 in C.)
The star_info structure contains the stellar model itself (i.e.,
zoning information, thermodynamic profile, composition profile).
These components are listed in the file
$MESA_DIR/star/public/star_data.inc. In addition, star_info
contains the values for the parameters that you set in your controls
inlist (i.e., initial_mass, xa_central_lower_limit). Recall that
the list of controls is located in $MESA_DIR/star/defaults/controls.defaults.
There is one set of controls that will prove useful time and time
again when using run_star_extras.f and that is x_ctrl,
x_integer_ctrl, and x_logical_ctrl. These are arrays (of length
100 by default) of double precision, integer, and boolean values. You
can set the elements in your inlists
&controls x_ctrl(1) = 3.14 x_ctrl(2) = 2.78 x_integer_ctrl(1) = 42 x_logical_ctrl(1) = .true. / ! end of controls inlist
and access them later on as part of the star structure (i.e., s%
x_ctrl(1), etc.).
Part 2a: Monitoring your models
Task 0 (Example): Add a stopping condition
If you assume that the Earth is a perfect blackbody, its equilibrium temperature is given by
\begin{equation*} T_\oplus = T_\odot \left(\frac{R_\odot}{2\,\rm AU}\right)^{1/2} \end{equation*}
Suppose the stellar model we're evolving represents the Sun and I want
to stop my calculation when the Earth would reach a given temperature.
A look through controls.defaults seems to indicate that such a
condition doesn't already exist. How do I do this?
First, look at how the routines in run_star_extras.f fit into a MESA
run. To decide whether to stop, I want to check the value of the
Earth's temperature after each step. Thus, I want the subroutine that
is called after each step, which is extras_finish_step.
Now, I need to figure out how to access information about the
conditions at the stellar photosphere. I open up
star/public/star_data.inc and start looking around. If I search for
the word photosphere, I can find what I'm looking for photosphere_r
and Teff.
MESA uses cgs units unless otherwise noted. The most common non-cgs
units are solar units. MESA defines its constants in
$MESA_DIR/const/public/const_def.f. Since the run_star_extras
module includes the line use const_def, we will be able to access
these values. Using the built in constants lets us make sure we're
using exactly the same definitions as MESA. The constant with the
value of the solar radius (in cm) is named Rsun. Note the other
constants that are defined.
! 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 > 300) 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
Now, recompile your working directory
./mk
You will need to do this step each and every time you edit
run_star_extras.f.
Now start the model again from the beginning
./rn
This run should halt around step 19.
Task 1: Allow the user to specify a temperature
Stop when the temperature of Earth exceeds a given value. Allow the user to specify this value in the inlist. A good value to specify is 330K, as it should take your model about 100 steps to reach this value.
Edit your inlist_project and comment out the central hydrogen
abundance stopping condition that was included. We won't use it
again.
You can receive valuable MESA bonus points if your routine stops when the temperature of the Earth is within one part in a million of the specified temperature. For its built-in stopping conditions, MESA provides the ability to control the absolute and/or relative error between the model and the specified stopping target. You can look at these controls, when_to_stop_rtol and when_to_stop_atol, for inspiration on how to do this.
Answer
Modify the termination condition in extras_finish_step to be
! stop if it exceeds a user-specified temperature (in K) if (Tearth > s% x_ctrl(1)) extras_finish_step = terminate
and then specify x_ctrl in your inlist
! stop when Tearth > this value x_ctrl(1) = 330
Bonus Answer
The basic code stops at the first timestep where the Earth temperature exceeds the threshold. In order to stop very near the threshold, we can ask MESA to "redo" any step that causes us to exceed the threshold by an amount greater than some tolerance.
This gives us the opportunity to use extras_check_model. If the
temperature is greater than the target temperature by more than our
threshold, we reject the step and redo it with half of the previous
timestep.
We define a few variables
real(dp) :: T, dT, delta real(dp), parameter :: epsilon = 1d-6
and then the logic itself straightforward
T = s% Teff * sqrt(s% photosphere_r * Rsun / (2.0 * AU)) dT = T - s% x_ctrl(1) delta = dT / s% x_ctrl(1) if (delta > 0) then if (delta > epsilon) then extras_check_model = redo s% dt = 0.5d0 * s% dt endif endif
You could also do this with retry instead of redo. In a retry, MESA
will do the step again with the timestep scaled down by
timestep_factor_for_retries and then not allow any timestep
increases for retry_hold steps afterwards.
Part 2b: Changing input physics
MESA provides hooks to override or modify many of its built-in routines. (These routines mostly affect things that occur within "take step" box of the flowchart.) These are referred to as "other" routines. There are two main steps needed to take advantage of this functionality: (1) writing the other routine and (2) instructing MESA to use this routine.
Navigate to $MESA_DIR/star/other, where you will see a set of files
named with the pattern other_*.f. In general, find the one
corresponding to the physics (or numerics) that you want to alter.
Open it up and read through it. Many of the files contain comments
and examples.
Note that we do not want to directly edit these files. Instead we
want to copy the template routine into our copy of run_star_extras.f
and then further modify it there. The template routines are usually
named either null_other_* or default_other_*.
In this example, we will focus on other_energy.f. Open up this file.
Copy the subroutine default_other_energy and paste it into your
run_star_extras.f. It should be at the same "level" as the other
subroutines in that file (that is, contained within the
run_star_extras module.).
subroutine default_other_energy(id, ierr) use const_def, only: Rsun integer, intent(in) :: id integer, intent(out) :: ierr type (star_info), pointer :: s integer :: k ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return s% extra_heat(:) = s% extra_power_source return ! here is an example of calculating extra_heat for each cell. do k = 1, s% nz if (s% r(k) > 0.7*Rsun .and. s% r(k) < 0.9*Rsun) then s% extra_heat(k) = 1d3*exp(-10*(s% r(k) - 0.8*Rsun)**2) end if end do end subroutine default_other_energy
The variable s% extra_heat is an additional specific (per mass)
heating rate that will be included. Note that this routine already
does something; default_other_energy is responsible for making the
extra_power_source control work. Go ahead and remove that bit, the
existing example (kudos if you spot the error), and rename it to
lecture1_other_energy.
In Fortran, you can write expressions that operate on the whole array
at once (like s% extra_heat(:) = s% extra_power_source). However,
it is often simplest to explicitly set the value of extra_heat (or
some other array) one value at a time, by using a loop. While we're
looking code with a loop, it is a good time to mention that in MESA,
the outermost zone is at k=1 and the innermost zone is at k=s% nz.
subroutine lecture1_other_energy(id, ierr) integer, intent(in) :: id integer, intent(out) :: ierr type (star_info), pointer :: s integer :: k ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return end subroutine lecture1_other_energy
If you read the comments in other_energy.f (and you should), you can
see that the file tells us how to have MESA use our other_* routine.
Perform these steps (hint: you will need to edit both your
run_star_extras.f and your inlists).
Task 2: Add an extra energy source to the core of the star
Use the other_energy routine to add a heating term
where \(M_r\) is the enclosed mass. Good values are \(\Delta M = 0.05 M_\odot\) and \(L_{\mathrm{extra}} = 0.1 L_\odot\).
The lower left panel in the PGSTAR plots displays the value of s%
extra_heat, so you should be able to easily check if it looks OK.
You can receive valuable MESA bonus points if your routine allows for user-specified values of \(\Delta M\) and \(L_{\mathrm{extra}}\).
Answer
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_energy = .true.
Second, edit the extras_controls routine in run_star_extras.f to
point s% other_energy 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_energy => lecture1_other_energy ... end subroutine extras_controls
Failure to do perform both of these is the most common problem people encounter when using the other_* hooks.
subroutine lecture1_other_energy(id, ierr) integer, intent(in) :: id integer, intent(out) :: ierr type (star_info), pointer :: s integer :: k real(dp) :: dM, Mr ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return dM = s% x_ctrl(3) * Msun do k = 1, s% nz ! m(k) is the enclosed mass at the outer cell edge ! so the mass coordinate at the middle of the cell is Mr = s% m(k) - 0.5 * s% dm(k) s% extra_heat(k) = s% x_ctrl(2) * Lsun * exp(-Mr/dM)/dM end do write(*,*) "Added ", dot_product(s% extra_heat(1:s%nz), s% dm(1:s%nz))/Lsun, & s% x_ctrl(2) end subroutine lecture1_other_energy
Part 2c: Analyzing your models
It is often useful to do some of your analysis in run_star_extras.
At runtime, you have access to more information about the star than
will be in the history and profile columns.
Task 3: Track the total amount of extra energy added
At the end of the run, print the total amount of energy added due to
the other_energy routine.
You can receive valuable MESA bonus points if your routine works even
if you do a restart (e.g., ./re x050). Hint: there are a number of
routines/variables in run_star_extras.f that contain the string
extra_info. Look at where/when each of these are called and what
they do.
Answer
First, define a module-level variable (the declaration goes between
the implicit none and the contains) to keep track of the total
amount of energy added.
real(dp) :: total_extra_energy
We want the initial value of total_extra_energy to be zero.
Therefore, we add the line
total_extra_energy = 0
our extras_startup routine.
In order to track the total amount of energy added, we can add the
amount of energy deposited in the current step in
extras_finish_step.
total_extra_energy = total_extra_energy + dot_product(s% dm(1:s% nz), s% extra_heat(1:s% nz)) * s% dt
Now, at the end of the run we want to write out this value. In
extras_after_evolve we can add the line
write(*,*) 'Energy added (ergs): ', total_extra_energy
Bonus Answer
In order to ensure that a variable is preserved across restarts, add a
call to move_dbl in move_extra_info.
i = 1 call move_dbl(total_extra_energy)
There is analogous code for integers (move_int) and for logicals
(move_flg).
The routine move_extra_info is called in three different ways.
In extras_startup there is the code
if (.not. restart) then call alloc_extra_info(s) else ! it is a restart call unpack_extra_info(s) end if
When you start a run, there is a call to alloc_extra_info. This
ensures that the extra storage space for these variables is allocated.
When you restart a run, there is a call to unpack_extra_info. This
loads the data from the photo into the variables that you're using.
In extras_finish_step there is a call to store_extra_info. This
is responsible for moving the values of your variables into the extra
storage space. (And then whatever is in the extra storage space is
written in the photo file when it is saved.)
Part 2d: Changing controls
Recall that star_info contains the values for the parameters that
you set in your controls inlist. That also means that you can set
the value of these parameters by modifying the star_info structure.
Since run_star_extras gives us hooks to access to the star_info at
each step, that means we can modify parameters as the run proceeds.
This often saves us the hassle of stopping, saving a model, editing
the inlist, and restarting.
Task 4: Turn on other_energy at central hydrogen exhaustion
Instead of having the other energy routine always on, activate it only after central hydrogen exhaustion has occurred.
Answer
In order to activate the other_energy routine, we add the following
line to extras_finish_step
! activate other_energy after central hydrogen depletion if (s% center_h1 < 0.01 ) s% use_other_energy = .true.
So after the end of the first step where the central hydrogen abundance falls below 0.01, the extra heating will occur.
Part 3: Evolving binary stars
The binary capabilities have been a major focus of MESA development in recent years.
Part 3a: Getting started
Each time you want to start a MESA binary project, you should make a new copy of the binary/work directory.
cp -r $MESA_DIR/binary/work lecture1-binary
In this case, we have prepared and provided a work directory for you. Download, unpack, and enter this work directory.
unzip lecture1-binary.zip cd lecture1-binary
The contents of the binary/work directory should look similar to a
standard star work directory, only doubled. There are now two inlists
(called inlist1 and inlist2) and two LOGS directories (called
LOGS1 and LOGS2).
The file inlist_project now contains the binary namelists
&binary_job and &binary_controls. These have analogous roles to
the &star_job and &controls namelists in a regular star/work
directory.
There are only a few possible controls in &binary_job, which are
documented in $MESA_DIR/binary/defaults/binary_job.defaults. This is
where you specify the inlists for each of the stellar models (via
inlist_names). Those files will be regular inlists for MESA/star,
where you can load/save models, change nuclear networks, do what
you're used to doing in MESA/star. You also choose whether to evolve
both stars (one star could be a point mass) and whether to follow
Roche lobe overflow. In this example, our binary_job namelist is
&binary_job ! each star has its own inlists inlist_names(1) = 'inlist1' inlist_names(2) = 'inlist2' ! in this example, we will evolve stellar models for both stars evolve_both_stars = .true. / ! end of binary_job namelist
The options in the &binary_controls inlists, which may be more
unfamiliar, are documented in
$MESA_DIR/binary/defaults/binary_controls.defaults. Most
importantly, this is where we set the initial orbit of the binary and
determine how mass and angular momentum are transferred.
In this example, our binary_controls namelist is
&binary_controls ! since we load stellar models, we do not need to specify masses ! if we were using a point mass we would need to set its mass ! we need to set the initial orbital properties of the binary initial_period_in_days = 0.027777777777777778d0 ! 40min ! choose the types of angular momentum losses we consider do_jdot_mb = .false. do_jdot_gr = .true. do_jdot_ml = .true. do_jdot_ls = .true. ! make frequent terminal output terminal_interval = 5 write_header_frequency = 1 / ! end of binary_controls namelist
Again, binary behaves much like star, so before you can run you
must issue the command
./mk
Task 1: Evolve a binary (two stars and their orbit)
The provided directory evolves an 0.4 \(M_\odot\) helium star in an initially 40 min orbit with an 0.8 \(M_\odot\) white dwarf. Look at the terminal output and identify what belongs to one star, what belongs to other star, and what belongs to the orbit. (You may want to save the terminal output by redirecting stdout to a file, using tee, etc.)
Answer
The terminal output is
DATE: 2017-08-01
TIME: 17:23:29
read inlist_project
version_number 9793
The terminal output contains the following information
'step' is the number of steps since the start of the run,
'lg_dt' is log10 timestep in years,
'age_yr' is the simulated years since the start run,
'lg_Tcntr' is log10 center temperature (K),
'lg_Dcntr' is log10 center density (g/cm^3),
'lg_Pcntr' is log10 center pressure (ergs/cm^3),
'Teff' is the surface temperature (K),
'lg_R' is log10 surface radius (Rsun),
'lg_L' is log10 surface luminosity (Lsun),
'lg_LH' is log10 total PP and CNO hydrogen burning power (Lsun),
'lg_L3a' is log10 total triple-alpha helium burning power (Lsun),
'lg_LZ' is log10 total burning power excluding LH and L3a and photodisintegrations (Lsun),
'lg_LNuc' is log10 nuclear power excluding photodisintegration (Lsun),
'lg_LNeu' is log10 total neutrino power (Lsun),
'lg_Psurf' is log10 surface pressure (gas + radiation),
'Mass' is the total stellar baryonic mass (Msun),
'lg_Mdot' is log10 magnitude of rate of change of mass (Msun/year),
'lg_Dsurf' is log10 surface density (g/cm^3),
'H_env' is the amount of mass where H is the most abundant iso,
'He_core' is the largest mass where He is most abundant iso.
'C_core' is the largest mass where C is most abundant iso.
'H_cntr' is the center H1 mass fraction,
'He_cntr' is the center He4 mass fraction,
'C_cntr' is the center C12 mass fraction,
'N_cntr' is the center N14 mass fraction,
'O_cntr' is the center O16 mass fraction,
'Ne_cntr' is the center Ne20 mass fraction,
'X_avg' is the star average hydrogen mass fraction,
'Y_avg' is the star average helium mass fraction,
'Z_avg' is the star average metallicity,
'gam_cntr' is the center plasma interaction parameter,
'eta_cntr' is the center electron degeneracy parameter,
'zones' is the number of zones in the current model,
'iters' is the number of newton iterations for the current step,
'retry' is the number of step retries required during the run,
'bckup' is the number of step backups required during the run,
'dt_limit' is an indication of what limited the timestep.
All this and more are saved in the LOGS directory during the run.
load saved model HeStar_0.4Msun.mod
set_initial_age 0.0000000000000000D+00
set_initial_model_number 0
net name basic.net
extra_terminal_output_file: star1.log
kappa_file_prefix gs98
kappa_lowT_prefix lowT_fa05_gs98
eos_file_prefix mesa
OMP_NUM_THREADS 2
The terminal output contains the following information
'step' is the number of steps since the start of the run,
'lg_dt' is log10 timestep in years,
'age_yr' is the simulated years since the start run,
'lg_Tcntr' is log10 center temperature (K),
'lg_Dcntr' is log10 center density (g/cm^3),
'lg_Pcntr' is log10 center pressure (ergs/cm^3),
'Teff' is the surface temperature (K),
'lg_R' is log10 surface radius (Rsun),
'lg_L' is log10 surface luminosity (Lsun),
'lg_LH' is log10 total PP and CNO hydrogen burning power (Lsun),
'lg_L3a' is log10 total triple-alpha helium burning power (Lsun),
'lg_LZ' is log10 total burning power excluding LH and L3a and photodisintegrations (Lsun),
'lg_LNuc' is log10 nuclear power excluding photodisintegration (Lsun),
'lg_LNeu' is log10 total neutrino power (Lsun),
'lg_Psurf' is log10 surface pressure (gas + radiation),
'Mass' is the total stellar baryonic mass (Msun),
'lg_Mdot' is log10 magnitude of rate of change of mass (Msun/year),
'lg_Dsurf' is log10 surface density (g/cm^3),
'H_env' is the amount of mass where H is the most abundant iso,
'He_core' is the largest mass where He is most abundant iso.
'C_core' is the largest mass where C is most abundant iso.
'H_cntr' is the center H1 mass fraction,
'He_cntr' is the center He4 mass fraction,
'C_cntr' is the center C12 mass fraction,
'N_cntr' is the center N14 mass fraction,
'O_cntr' is the center O16 mass fraction,
'Ne_cntr' is the center Ne20 mass fraction,
'X_avg' is the star average hydrogen mass fraction,
'Y_avg' is the star average helium mass fraction,
'Z_avg' is the star average metallicity,
'gam_cntr' is the center plasma interaction parameter,
'eta_cntr' is the center electron degeneracy parameter,
'zones' is the number of zones in the current model,
'iters' is the number of newton iterations for the current step,
'retry' is the number of step retries required during the run,
'bckup' is the number of step backups required during the run,
'dt_limit' is an indication of what limited the timestep.
All this and more are saved in the LOGS directory during the run.
set_tau_factor 3.0000000000000000D+02
load saved model wd_0.8Msun.mod
set_initial_age 0.0000000000000000D+00
set_initial_model_number 0
change to co_burn.net
number of species 9
new_v_flag T
net name co_burn.net
extra_terminal_output_file: star2.log
v_flag T
tau_factor 3.0000000000000000D+02
kappa_file_prefix gs98
kappa_lowT_prefix lowT_fa05_gs98
eos_file_prefix mesa
OMP_NUM_THREADS 2
m2 8.0000000000000004D-01
m1 1.0000000000000000D+00
initial_period_in_days 2.7777777777777780D-02
initial_separation_in_Rsun 4.1000206340324563D-01
jdot_multiplier 1.0000000000000000D+00
fr 1.0000000000000001D-01
The binary terminal output contains the following information
'step' is the number of steps since the start of the run,
'lg_dt' is log10 timestep in years,
'age_yr' is the simulated years since the start run,
'M1+M2' is the total mass of the system (Msun),
'M1' is the mass of the primary (Msun)
'M2' is the mass of the secondary (Msun)
'separ' is the semi-major axis of the orbit (Rsun),
'R1' is the radius of the primary (Rsun)
'R2' is the radius of the secondary (Rsun)
'Porb' is the orbital period (days),
'P1' is the rotation period of star 1 (days, zero if not modeling rotation),
'P2' is the rotation period of star 2 (days, zero if not modeling rotation),
'e' orbital eccentricity,
'dot_e' time derivative of e (1/yr),
'Eorb' orbital energy G*M1*M2/2*separation (ergs),
'M2/M1' mass ratio,
'vorb1' orbital velocity of star 1 (km/s),
'vorb2' orbital velocity of star 2 (km/s),
'pm_i' index of star evolved as point mass, zero if both stars are modeled,
'RL1' Roche lobe radius of star 1 (Rsun),
'Rl2' Roche lobe radius of star 2 (Rsun),
'donor_i' index of star taken as donor,
'RL_gap1' (R1-Rl1)/Rl1,
'RL_gap2' (R2-Rl2)/Rl2,
'dot_Mmt', mass transfer rate (Msun/yr),
'dot_M1', time derivative for the mass of star 1 (Msun/yr),
'dot_M2', time derivative for the mass of star 2 (Msun/yr),
'eff', mass transfer efficiency, computed as -dot_M2/dot_M1 (zero if dot_M1=0),
'dot_Medd', Eddington accretion rate (Msun/yr),
'L_acc', accretion luminosity when accreting to a point mass (ergs/s),
'Jorb', orbital angular momentum (g*cm^2/s)
'spin1', spin angular momentum of star 1 (g*cm^2/s),
'spin2', spin angular momentum of star 2 (g*cm^2/s),
'dot_J', time derivative of Jorb (g*cm^2/s^2),
'dot_Jgr', time derivative of Jorb due to gravitational waves (g*cm^2/s^2),
'dot_Jml', time derivative of Jorb due to mass loss (g*cm^2/s^2),
'dot_Jmb', time derivative of Jorb due to magnetic braking (g*cm^2/s^2),
'dot_Jls', time derivative of Jorb due to spin-orbit coupling (g*cm^2/s^2),
'rlo_iters', number of iterations for implicit calculation of mass transfer,
All this and more can be saved in binary_history.data during the run.
save LOGS1/profile1.data for model 1
save LOGS2/profile1.data for model 1
__________________________________________________________________________________________________________________________________________________
step lg_Tcntr Teff lg_LH lg_Lnuc Mass H_rich H_cntr N_cntr Y_surf X_avg eta_cntr zones retry
lg_dt_yr lg_Dcntr lg_R lg_L3a lg_Lneu lg_Mdot He_core He_cntr O_cntr Z_surf Y_avg gam_cntr iters bckup
age_yr lg_Pcntr lg_L lg_LZ lg_Psurf lg_Dsurf C_core C_cntr Ne_cntr Z_cntr Z_avg v_div_cs dt_limit
__________________________________________________________________________________________________________________________________________________
5 8.024590 3.441E+04 -8.888221 0.293805 0.400000 0.000000 0.000000 0.012389 0.980532 0.000000 1.718213 941 0
3.540366 4.647510 -1.045299 0.280588 -1.926298 -99.000000 0.400000 0.979729 0.000744 0.019468 0.980488 0.196749 6 0
1.4099E+04 20.566539 1.009370 -1.229450 5.124406 -6.999066 0.000000 0.000929 0.002123 2.027E-02 1.951E-02 -0.407E-10 max increase
__________________________________________________________________________________________________________________________________________________
step lg_Tcntr Teff lg_LH lg_Lnuc Mass H_rich H_cntr N_cntr Y_surf X_avg eta_cntr zones retry
lg_dt_yr lg_Dcntr lg_R lg_L3a lg_Lneu lg_Mdot He_core He_cntr O_cntr Z_surf Y_avg gam_cntr iters bckup
age_yr lg_Pcntr lg_L lg_LZ lg_Psurf lg_Dsurf C_core C_cntr Ne_cntr Si_cntr Z_avg v_div_cs dt_limit
__________________________________________________________________________________________________________________________________________________
2 5 7.474487 2.325E+04 -11.503953 -11.503953 0.799507 0.000000 0.000000 0.000000 0.980462 0.000000 202.064899 970 0
3.540366 7.019162 -1.988643 -29.858887 -2.710354 -99.000000 0.799507 0.000000 0.601497 0.019444 0.006412 35.254404 6 0
1.4099E+04 23.947243 -1.558301 -18.115684 7.836885 -4.252208 0.790012 0.372767 0.007869 0.002574 9.936E-01 0.293E-12 max increase
__________________________________________________________________________________________________________________________________________________
binary_step M1+M2 separ Porb e M2/M1 pm_i donor_i dot_Mmt eff Jorb dot_J dot_Jmb
lg_dt M1 R1 P1 dot_e vorb1 RL1 Rl_gap1 dot_M1 dot_Medd spin1 dot_Jgr dot_Jls
age_yr M2 R2 P2 Eorb vorb2 RL2 Rl_gap2 dot_M2 L_acc spin2 dot_Jml rlo_iters
__________________________________________________________________________________________________________________________________________________
bin 5 1.199507 0.409871 0.027765 0.000E+00 1.998767 0 1 0.000E+00 1.000000 1.130E+51 -4.049E+35 0.000E+00
3.540366 0.400000 0.090095 0.000000 0.000E+00 498.012330 0.131502 -3.149E-01 0.000E+00 1.000E+99 0.000E+00 -4.049E+35 0.000E+00
1.4099E+04 0.799507 0.010265 0.000000 -1.480E+48 249.159709 0.180323 -9.431E-01 0.000E+00 0.000E+00 0.000E+00 0.000E+00 1
stop because model_number >= max_model_number
10 10
__________________________________________________________________________________________________________________________________________________
step lg_Tcntr Teff lg_LH lg_Lnuc Mass H_rich H_cntr N_cntr Y_surf X_avg eta_cntr zones retry
lg_dt_yr lg_Dcntr lg_R lg_L3a lg_Lneu lg_Mdot He_core He_cntr O_cntr Z_surf Y_avg gam_cntr iters bckup
age_yr lg_Pcntr lg_L lg_LZ lg_Psurf lg_Dsurf C_core C_cntr Ne_cntr Z_cntr Z_avg v_div_cs dt_limit
__________________________________________________________________________________________________________________________________________________
10 8.030494 3.449E+04 -8.836364 0.529338 0.400000 0.000000 0.000000 0.012387 0.980532 0.000000 1.684673 941 0
3.936272 4.648437 -1.049135 0.515899 -1.905243 -99.000000 0.400000 0.979606 0.000744 0.019468 0.980467 0.194258 6 0
4.5089E+04 20.571752 1.005753 -0.986788 5.127909 -6.995673 0.000000 0.001051 0.002126 2.039E-02 1.953E-02 -0.403E-10 max increase
save LOGS1/profile2.data for model 10
save photos1/x010 for model 10
save LOGS2/profile2.data for model 10
save photos2/x010 for model 10
bin 10 1.199507 0.409584 0.027735 0.000E+00 1.998767 0 1 0.000E+00 1.000000 1.130E+51 -4.057E+35 0.000E+00
3.936272 0.400000 0.089303 0.000000 0.000E+00 498.187060 0.131410 -3.204E-01 0.000E+00 1.000E+99 0.000E+00 -4.057E+35 0.000E+00
4.5089E+04 0.799507 0.010270 0.000000 -1.481E+48 249.247129 0.180196 -9.430E-01 0.000E+00 0.000E+00 0.000E+00 0.000E+00 1
termination code: max_model_number
10 8.030494 3.449E+04 -8.836364 0.529338 0.400000 0.000000 0.000000 0.012387 0.980532 0.000000 1.684673 941 0
3.936272 4.648437 -1.049135 0.515899 -1.905243 -99.000000 0.400000 0.979606 0.000744 0.019468 0.980467 0.194258 6 0
4.5089E+04 20.571752 1.005753 -0.986788 5.127909 -6.995673 0.000000 0.001051 0.002126 2.039E-02 1.953E-02 -0.403E-10 max increase
2 10 7.474409 2.390E+04 -11.404806 -11.404806 0.799507 0.000000 0.000000 0.000000 0.980462 0.000000 202.101536 970 0
3.936272 7.019164 -1.988439 -29.867804 -2.711123 -99.000000 0.799507 0.000000 0.601497 0.019444 0.006412 35.260767 6 0
4.5089E+04 23.947245 -1.509920 -18.122845 7.877797 -4.259287 0.790012 0.372767 0.007869 0.002574 9.936E-01 0.550E-13 max increase
DATE: 2017-08-01
TIME: 17:23:34
Part 3b: Using run_binary_extras
MESA evolves a binary system by independently solving the structure of each component star and the orbital parameters, all using the same timestep. The timestep could be limited by changes in either star or by changes in the binary properties.
In the src/ subdirectory of a binary work directory there is a file
run_binary_extras.f. This provides routines analogous to those in
run_star_extras.f, but that apply to the binary evolution.
Importantly, there is still a run_star_extras.f with all of the same
routines we know and love.
Task 2: Determine the order of star and binary extras routines
Now that we have to think about the routines in both
run_star_extras.f and run_binary_extras.f, we should find out what
order they are called in. Add write statements in the binary routines
extras_binary_check_model and extras_binary_finish_step. Also add
write statements in the star routines extras_check_model and
extras_finish_step; for these routines, also print the value of
their argument named id. You may also find it helpful to print some
additional identifying information about the star (e.g. its mass).
Run MESA and use your terminal output to understand what happens when.
Answer
When I look at my terminal output, I see
check model for star 1 check model for star 2 check model for binary finish step for binary finish step for star 1 finish step for star 2
This tells us that check_model is called first for each of the stars
and then for the binary. Since the stars and the orbit are evolved
with the same timestep, if any of these asked for a redo/retry/backup,
MESA will have to recompute all three.
This also tells us that finish_step is called first for the binary
and then for each of the stars.