MESA-Web History Files¶

Preliminaries¶

In the previous labs, you’ve explored the basics of Python and Jupyter notebooks, and undertaken a MESA-Web calculation. In this lab, you’ll use Jupyter to plot data from the history file created during your calculation.

How does MESA Work?¶

Before we get stuck in, let’s briefly review how the MESA stellar evolution code works. To solve the stellar structure equations (Handout 18), MESA discretizes them. This numerical technique involves replacing continuous quantities (e.g., the radial coordinate \(r\)) with a set of \(N\) discrete values (e.g., \(r=r_{1},r_{2},\ldots,r_{N}\)), representing the quantity sampled at discrete locations \(m=m_{1},m_{2},\ldots,m_{N}\) throughout the star (remember that the interior mass \(m\) is the independent variable in the stellar structure equations). By convention, MESA orders these locations outside-in; so, the first point is at the stellar surface (\(m_{1}=M\)), and subsequent points are ordered such that \(m_{i} < m_{i-1}\) (\(i=2,\ldots,N\)).

In addition to this spatial discretization of the structure equations, MESA uses a temporal discretization of the evolution equations (again, see Handout 18). Thus, it evaluates the structure of the star at a set of discrete times \(t=t_{1},t_{2},\ldots\) known as ‘timesteps’, or just ‘steps’. The time difference \({\rm d}t = t_{j} - t_{j-1}\) between adjacent steps is varied on-the-fly, in order to ensure that changes in the stellar structure are accurately calculated. Typically, \({\rm d}t\) is long when the star is changing slowly with time, and short when the star is changing rapidly with time.

What does MESA Output?¶

During a stellar evolution calculation, MESA writes out two types of data file:

Profile files contain snapshots of the stellar structure at a single timestep — i.e., the interior mass \(m_{i}\), radial coordinate \(r_{i}\), pressure \(P_{i}\), etc. (\(i=1,\ldots,N\)). There can be multiple profile files written during the calculation, representing the structure at different timesteps.
History files contain global properties of the star at each timestep of the calculation — i.e., the age \(t_{j}\), mass \(M_{j}\), effective temperature \(T_{{\rm eff},j}\) etc. (\(j=1,2,\ldots\)). There is only a single history file written during the calculation.

What does MESA-Web Output?¶

Although MESA allows a great deal of flexibility in choosing which specific data are written to profile and history files, MESA-Web simplifies things for the user by always writing a standard set of data. This means we always know exactly what to expect in profile files and history files produced by MESA-Web.

Exploring a History File¶

Installing the mesa-web Module¶

To help with reading and analyzing the profile and history files created by MESA-Web, your first step is to install the mesa_web Python module, which was created by Rich (so, report bugs to him!). Download the module from here, and save it in your astro_310 folder.

Creating a Notebook¶

Next, fire up Jupyter in your web browser and create a new notebook in the astro_310 folder, naming it history-notebook. In this notebook, you’ll explore the contents of the history file produced during your MESA-Web calculation. Start off by importing modules:

# Import the numpy module to provide numerical functionality
import numpy as np

# Import the matplotlib.pyplot module to provide plotting functionality
import matplotlib.pyplot as plt

# Tell matplotlib.pyplot to do inline plots
%matplotlib inline

# Import the mesa-web module to simplify reading MESA-Web files
import mesa_web as mw

(as usual, cut and paste this code into a notebook cell, and then execute the cell).

Reading History Data¶

Now, customize the following code to read in the file trimmed_history.data from the folder you unpacked:

# Read history data. Be sure to replace the NNNNNNNNNN with the
# specific digits of your folder

hist_data = mw.read_history('MESA-Web_Job_NNNNNNNNNN/trimmed_history.data')

# Inspect the hist_data variable

print(type(hist_data))
print(hist_data.keys())

From the last two commands, we see that hist_data is a dict containing many key/value pairs. To understand what’s actually contained in these pairs, we can use Python’s help() command to access the documentation for the mw.read_history function:

# Print out documentation for the read_history function

help(mw.read_history)

From reading through this documentation, we see that the dict contains a mix of header data and history data. Header data are single values, such as the date on which the MESA-Web calculation began. History data are 1-dimensional arrays that record properties of the star — e.g., its age, its luminosity, its radius — at each timestep of the calculation.

Plotting History Data¶

Let’s now extract data from the hist_data dict and create some plots exploring the evolution of the MESA-Web model. We’ll start with an evolutionary track in the Hertzsprung-Russell diagram:

# Extract data from hist_data, using dict indexing

log_Teff = hist_data['log_Teff'] # log(Teff/K)
log_L = hist_data['log_L']       # log(L/Lsun)

# Plot the HR diagram

plt.figure()

plt.plot(log_Teff, log_L)

plt.xlim(4.0, 3.5)
plt.ylim(-1,4)

plt.xlabel('log (T/K')
plt.ylabel('log (L/Lsun)')

plt.title('Hertzsprung-Russel Diagram for Sun')

Let’s now augment this plot by adding a marker for the present-day Sun. To do this, we must first choose a criterion for what constitutes the present-day Sun; and then find the index \(j\) of the timestep where the history data that best matches this criterion. Let’s set the criterion as \(L_{\rm H} = L_{\odot}\), where \(L_{\rm H}\) is the hydrogen burning luminosity (aside: have a think why \(L = L_{\odot}\) wouldn’t suffice). Then, we can set up the index thus:

# Extract log(LH) from hist_data

log_LH = hist_data['log_LH']

# Find the timestep index where log(LH) is closest to zero (i.e., LH closest
# to 1 Lsun), as representative of the present-day Sun. The np.abs()
# function returns the absolute value. The np.argmin() function
# returns the index of the smallest element (the final + 1 is because Python
# uses zero-based indexing).

j_pres = np.argmin(np.abs(log_LH)) + 1

print(j_pres)

# Print out values at this timestep (the -1 is again because of the zero-based indexing)

print('log(L_H) at present     :', log_LH[j_pres-1])
print('log(L) at present       :', log_L[j_pres-1])
print('log(Teff) at present    :', log_Teff[j_pres-1])
print('log(R) at present       :', hist_data['log_R'][j_pres-1])

Note

In the final line of the above code, note how we perform double indexing on hist_data: first to access an array within a dict, and then to access an element of the array).

Then, our HR diagram plot including the present-day Sun becomes this:

# Plot the HR diagram, marking the present-day Sun

plt.figure()

plt.plot(log_Teff, log_L)

plt.scatter(log_Teff[j_pres-1], log_L[j_pres-1], color='r', label='Present Day')

plt.xlim(4.0, 3.5)
plt.ylim(-1,4)

plt.xlabel('log (T/K')
plt.ylabel('log (L/Lsun)')

plt.title('Hertzsprung-Russel Diagram for Sun')

plt.legend()

We’re not just limited to HR diagrams; we can plot any one quantity against any other. Here’s a plot showing the central hydrogen and helium as a function of time. We can clearly see the depletion of hydrogen, and the enrichment of helium, as the Sun gets older!

# Plot central hydrogen & helium abundances versus time (measured
# since the start of the calculation)

X_c = hist_data['center_h1']
Y_c = hist_data['center_he4']
t = hist_data['star_age']

plt.figure()

plt.plot(t/1E9, X_c, color='b', label='X_c')
plt.plot(t/1E9, Y_c, color='r', label='Y_c')

plt.xlabel('t (Gyr)')

plt.legend()

We can also do the same plot against model number, which is an integer starting at 1 that labels MESA’s timesteps (this model number in fact equals the timestep index \(j\) we’ve mentioned previously).

# Plot central hydrogen abundance versus model number

X_c = hist_data['center_h1']
Y_c = hist_data['center_he4']
mod_num = hist_data['model_number']

plt.figure()

plt.plot(mod_num, X_c, color='b', label='X_c')
plt.plot(mod_num, Y_c, color='r', label='Y_c')

plt.xlabel('Model Number')

plt.legend()

Comparing these two plots (versus time, and versus model number) reminds us that MESA varies the time difference \({\rm d}t\) between steps. Typically, \({\rm d}t\) is long when the star is changing slowly with time, and short when the star is changing rapidly with time. We can see this most clearly by plotting time against model number:

# Plot time versus model number

t = hist_data['star_age']
mod_num = hist_data['model_number']

plt.figure()

plt.plot(mod_num, t/1E9, color='b')

plt.xlabel('Model Number')
plt.ylabel('t (Gyr)')

From this plot, we see many small steps are required to get through the pre-main sequence phase (\(1 \leq j \lesssim 250\)), while fewer, larger steps are required to traverse the main sequence (\(j \gtrsim 250\)).