MESA-Web Profile Files¶
Preliminaries¶
In the previous labs, you’ve explored the basics of Python and Jupyter notebooks, undertaken a MESA-Web calculation, and learned how to plot data from a MESA-Web history file. In this lab, you’ll use Jupyter to plot data from the profile files created during your calculation.
Exploring a Profile File¶
Using the Profiles Index¶
Profile files provide a snapshot of the full internal structure of the
model star — i.e., variables such as radial coordinate \(r\),
interior mass \(m\), temperature \(T\) — at a single
timestep. Because they can be quite large, MESA-Web only writes them
for selected timesteps. An index of the profile files, and the model
numbers they correspond to, is provided in the MESA-Web output folder
as the profiles.index file. Go ahead and take a look at this
index file for your calculation; you can do
this using Notepad (on Windows), TextEdit (on MacOS), or from the file
browser that Jupyter provides. You should see something like this:
11 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
368 3 11
The first column of this file lists model numbers for which corresponding profile files exist; and the third column gives the profile number (don’t worry about the second column for now).
Let’s figure out if there is a profile file for the present-day sun. In Lab III we found that the present-day Sun (defined by the point where \(L_{\rm H} = L_{\odot}\)) occurs in our calculation at model number 298. The index file shows that there isn’t a specific profile file for this model number; however, there is a file for model number 297 (which should have a very similar structure), and it has a profile number of 8.
Creating a Notebook¶
Fire up Jupyter in your web browser and create a new notebook in the
astro_310 folder, naming it profile-notebook. Then, import
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).
Then, customize the following code to read in the file
profileP.data from the MESA-Web folder (where P is the profile
number determined above).
# Read history data. Be sure to replace the NNNNNNNNNN with the
# specific digits of your folder, and P with the profile number
prof_data = mw.read_profile('MESA-Web_Job_NNNNNNNNNN/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_profile 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 interior mass \(m/M\):
# Plot abundances versus fractional interior 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.