Lab II: Python + MESA-Web¶
Overview¶
In this lab, you will use Jupyter to read in, plot and analyze data created by the MESA-Web calculation you submitted in the Intermission.
Unpacking the Data¶
When your MESA-Web calculation completes, you will receive an email
with a link to a Zip archive. Download this Zip archive, saving it
into the astro_310 folder you created in Lab I. (The default location of this folder is in
/Users/username on Windows and Mac OS, and in /home/username
on Linux, where username is your user name). If for some reason
you never received the email, or if you failed to download the Zip
archive with 24 hours of notification, you can grab a replacement
here.
Next, unpack the Zip archive into the astro_310 folder. This
should produce a new folder with the name MESA-Web_Job_MMDDNNNNNN,
where MM and DD are the two-digit month and day when you
submitted the job, and NNNNNN is a six-digit unique calculation
identifier. Inside the folder you’ll find a number of files:
trimmed_history.datais the history file for the calculation, storing global data as a function of time.profileP.data(wherePis one or more digits) are the profile files for the calculation, storing spatial data for selected evolutionary stages.input.txtis a summary of the parameters you gave to MESA-Web when you submitted the calculation.MESA-Web_Job_MMDDNNNNNN.mp4(whereMMDDNNNNNNis the same as in the folder name) is an MP4-format movie summarizing the star’s evolution.
Installing the mesa-web Module¶
Your next step is to install the mesa_web Python module, which was
created by Rich to simplify the task of reading the history and
profile files produced by MESA-Web. Download the module from here, and save it in the astro_310 folder.
Exploring History Files¶
Getting Started¶
Create a new Jupyter notebook in the astro_310 folder, naming it
history-notebook. In this notebook, you’ll explore the contents of
the history file produced during MESA-Web run. 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 MMDDNNNNNN with the
# specific digits of your folder
hist_data = mw.read_history('MESA-Web_Job_MMDDNNNNNN/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 initial stellar mass, that apply to the whole MESA-Web calculation. History data are 1-dimensional arrays that record how specific properties of the star — e.g., its age, its luminosity, its radius — change during 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 in 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 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
i_pres = np.argmin(np.abs(log_LH))
print(i_pres)
# Print out values at this index
print('log(L_H) at present :', log_LH[i_pres])
print('log(L) at present :', log_L[i_pres])
print('log(Teff) at present :', log_Teff[i_pres])
print('log(R) at present :', hist_data['log_R'][i_pres])
(note, in the final line, 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[i_pres], log_L[i_pres], 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 depletion of hydrogen, and the enrichment of helium, as the Sun gets older:
# Plot central hydrogen abundance versus age (measured since the
# start of the calculation)
X_c = hist_data['center_h1']
Y_c = hist_data['center_he4']
age = hist_data['star_age']
plt.figure()
plt.plot(age/1E9, X_c, color='b', label='X')
plt.plot(age/1E9, Y_c, color='r', label='Y')
plt.xlabel('Age')
plt.legend()
We can also do the same plot against model number, which is an integer starting at 1 that indexes how many discrete steps MESA has taken during its calculation:
# 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(age/1E9, X_c, color='b', label='X')
plt.plot(age/1E9, Y_c, color='r', label='Y')
plt.xlabel('Model Number')
plt.legend()
Comparing these two plots reveals that MESA doesn’t step forward in time at a uniform rate — it takes large steps to traverse evolutionary stages where the stellar structure is evolving slowly, and small steps to handle stages where the structure is evolving rapidly. In the present case, small steps are required as the model evolves up the red-giant branch.
Exploring Profile Files¶
Getting Started¶
Now we’re going to have a look at the profile files produced by MESA-Web. These files provide snapshots of the full internal structure of the model star — i.e., variables such as radial coordinate \(r\), interior mass \(m\), temperature \(T\) — at selected stages during its evolution.
To determine which profile file corresponds to which stage requires a
little work. Let’s demonstrate the process for the present-day
Sun. First, find the model number corresponding to the index
i_pres determined previously:
# Print out the model number representative of the present-day Sun
# (the array mod_num was defined in the previous cell)
print(mod_num[i_pres])
Then, open up the file profiles.index (in the MESA-Web output
folder). Here’s a sample from that file:
139 models. lines hold model number, priority, and profile number.
1 2 1
50 1 2
100 1 3
150 1 4
200 1 5
244 2 6
250 1 7
297 2 8
300 1 9
350 1 10
The first column of this file lists model numbers for which profile files exist; and the third column gives the profile index number (don’t worry about the second column for now). Your task is to find the model number that’s closest to the model number just printed out, and then read off the corresponding profile index. So, if we’re looking for model number 298 (say), we would pick the line with model number 297 (highlighted above), and read off the profile index number as 8.
Reading Profile Data¶
Create a new Jupyter notebook in the astro_310 folder, naming it
profile-notebook. Import modules as before:
# 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
Then, customize the following code to read in the file
profileP.data from the MESA-Web folder (where P is the profile
index number determined above).
# Read history data. Be sure to replace the MMDDNNNNNN with the
# specific digits of your folder, and P with the profile index number
prof_data = mw.read_profile('MESA-Web_Job_MMDDNNNNNN/profileP.data')
# Inspect the prof_data variable
print(type(prof_data))
print(prof_data.keys())
As was the case with hist_data, we see that prof_data is a
dict containing many key/value pairs, and we can use help()
command to access the documentation describing these data:
# Print out documentation for the read_history function
help(mw.read_profile)
Let’s verify that we read in the correct profile, by checking the model number:
# Check the model number
print(prof_data['model_number'])
Plotting Profile Data¶
We’re now in a position to make some plots. First, a simple plot showing how the pressure, density and temperature vary as a function of fractional radius \(r/R\) throughout the star:
# Define constants
RSUN = 6.957E10 # Solar radius, in cm
MSUN = 1.989E33 # Solar mass, in g
# Plot pressure, density and temperature versus fractional radius
P = prof_data['pressure']
rho = 10**prof_data['logRho']
T = 10**prof_data['logT']
r = prof_data['radius']*RSUN # Convert radii to cm
R = prof_data['photosphere_r']*RSUN # Same
plt.figure()
plt.plot(r/R, P, color='r', label='P (dyn/cm^2)')
plt.plot(r/R, rho, color='g', label='rho (g/cm^3)')
plt.plot(r/R, T, color='b', label='T (K)')
plt.xlabel('r/R')
plt.yscale('log')
plt.legend()
Note how we make the vertical axis be logarithmic by using the
plt.yscale('log') command.
We’re not limited to having radius on the horizontal axis; we can use any variable we like. Here a plot of the hydrogen, helium, carbon and nitrogen mass fractions versus fractional mass \(m/M\):
# Plot abundances versus fractional mass
X = prof_data['h1']
Y = prof_data['he4']
Z_C12 = prof_data['c12']
Z_N14 = prof_data['n14']
m = prof_data['mass']*MSUN # Convert masses to g
M = prof_data['star_mass']*MSUN # Same
plt.figure()
plt.plot(m/M, X, color='r', label='X')
plt.plot(m/M, Y, color='g', label='Y')
plt.plot(m/M, Z_C12, color='b', label='Z_C12')
plt.plot(m/M, Z_N14, color='m', label='Z_N14')
plt.xlabel('m/M')
plt.yscale('log')
plt.legend()
Analyzing Profile Data¶
To finish the lab, let’s do a little analysis on the profile data. Specifically, let’s see if we can confirm that the model for the present-day Sun satisfies the hydrostatic equilibrium equation
First, let’s calculate the left-hand and right-hand sides of this
equation. To evaluate the pressure gradient, we can use the
np.gradient function, which uses finite differences to evaluate a
numerical gradient estimate:
# Define constants
G = 6.674E-8 # Gravitational constant, in cm^3/g/s^2
# Evaluate lhs and rhs of the hydrostatic equilibrium equation
lhs = np.gradient(P, r)
rhs = -G*m/r**2 * rho
Now, create a plot comparing the two:
# Plot the lhs and rhs of the hydrostatic equilibrium equation
plt.figure()
plt.plot(r/R, lhs, color='b', label='lhs')
plt.plot(r/R, rhs, color='r', label='rhs', dashes=(2,1))
plt.xlabel('r/R')
plt.legend()
The correspondence between the two is reasonable, but the pressure gradient data are very noisy. This is a consequence of using finite differences to evaluate the gradient. Internally, MESA uses a more sophisticated approach, assuring that hydrostatic balance is achieved to a high degree of accuracy.