Interpreting Output Files
This chapter demonstrates using Python to read and plot the summary and detail output files written by the GYRE frontends. Further information about these files is provided in the Output Files chapter.
PyGYRE
PyGYRE is a Python package maintained separately from GYRE, that provides a set of routines that greatly simplify the analysis of summary and detail files. Detailed information about PyGYRE can be found in the full documentation; here, we demonstrate how to use it to read and plot the output files from the Example Walkthrough section.
As a preliminary step, you’ll need to install PyGYRE from the Python Package Index (PyPI). This can be done using the pip command,
pip install pygyre
If PyGYRE is already installed, you can upgrade to a more-recent version via
pip install --upgrade pygyre
Analyzing a Summary File
To analyze the summary file written by gyre during the example walkthrough, change into your work directory and fire up your preferred interactive Python environment (e.g., Jupyter). Import PyGYRE and the other modules needed for plotting:
# Import modules
import pygyre as pg
import matplotlib.pyplot as plt
import numpy as np
(you may want to directly cut and paste this code). Next, read the summary file into the variable s:
# Read data from a gyre summary file
s = pg.read_output('summary.h5')
The pygyre.read_output()
function is able to read
files in both TXT and HDF formats, returning the data in an
astropy.table.Table
object. To inspect the data
on-screen, simply evaluate the table:
# Inspect the data
s
From this, you’ll see that there are three columns in the table,
containing the harmonic degree l
, radial order n_pg
and
frequency freq
of each mode found during the GYRE run.
Next, plot the frequencies against radial orders via
# Plot the data
plt.figure()
plt.plot(s['n_pg'], s['freq'].real)
plt.xlabel('n_pg')
plt.ylabel('Frequency (cyc/day)')
(the values in the freq
column are complex, and we plot the real
part). The plot should look something like Fig. 8.
Fig. 8 The frequency \(\nu\) of \(\ell=1\) and \(\ell=2\)
modes, plotted against their radial order \(\numpg\).
(Source
)
The straight line connecting the two curves occurs because we are plotting both the dipole and quadrupole modes together. To separate them, the table rows can be grouped by harmonic degree:
# Plot the data, grouped by harmonic degree
plt.figure()
sg = s.group_by('l')
plt.plot(sg.groups[0]['n_pg'], sg.groups[0]['freq'].real, label=r'l=1')
plt.plot(sg.groups[1]['n_pg'], sg.groups[1]['freq'].real, label=r'l=2')
plt.xlabel('n_pg')
plt.ylabel('Frequency (cyc/day)')
plt.legend()
The resulting plot, in Fig. 9, looks much better.
Fig. 9 The frequency nu of \(\ell=1\) and \(\ell=2\)
modes, grouped by \(\ell\) and plotted against their radial order \(\numpg\).
(Source
)
Analyzing a Detail File
Now let’s take a look at one of the detail files, for the mode with
\(\ell=1\) and \(\numpg=-7\). As with the summary file,
pygyre.read_output()
can be used to read the file
data into an astropy.table.Table
object:
# Read data from a GYRE detail file
d = pg.read_output('detail.l1.n-7.h5')
Inspecting the data using
# Inspect the data
d
shows there are 7 columns: the fractional radius x
, the radial
displacement eigenfunction xi_r
, the horizontal displacement
eigenfunction xi_h
, and 4 further columns storing structure
coefficients (see the Detail Files section for descriptions of
these data). Plot the two eigenfunctions using the code
# Plot displacement eigenfunctions
plt.figure()
plt.plot(d['x'], d['xi_r'].real, label='xi_r')
plt.plot(d['x'], d['xi_h'].real, label='xi_h')
plt.xlabel('x')
plt.legend()
Fig. 10 The radial (\(\txir\)) and horizontal (\(\txih\))
displacement eigenfunctions of the \(\ell=1\), \(n_{\rm
pg}=-7\) mode, plotted against the fractional radius \(x\).
(Source
)
The plot should look something like Fig. 10. From this figure , we see that the radial wavelengths of the eigenfunctions become very short around a fractional radius \(x \approx 0.125\). To figure out why this is, we can take a look at the star’s propagation diagram:
# Evaluate dimensionless characteristic frequencies
l = d.meta['l']
omega = d.meta['omega']
x = d['x']
V = d['V_2']*d['x']**2
As = d['As']
c_1 = d['c_1']
Gamma_1 = d['Gamma_1']
d['N2'] = d['As']/d['c_1']
d['Sl2'] = l*(l+1)*Gamma_1/(V*c_1)
# Plot the propagation diagram
plt.figure()
plt.plot(d['x'], d['N2'], label='N^2')
plt.plot(d['x'], d['Sl2'], label='S_l^2')
plt.axhline(omega.real**2, dashes=(4,2))
plt.xlabel('x')
plt.ylabel('omega^2')
plt.ylim(5e-2, 5e2)
plt.yscale('log')
Note how we access the mode harmonic degree l
and dimensionless
eigenfrequency omega
through the table metadata dict
d.meta
. The resulting plot (Fig. 11) reveals
that the Brunt-Väisälä frequency squared is large around \(x
\approx 0.125\); this feature is a consequence of the molecular weight
gradient zone outside the star’s convective core, and results in the
short radial wavelengths seen there in Fig. 10.
Fig. 11 Propagation diagram for the \(5\,\Msun\) model, plotting the
squares of the Brunt-Väisälä (\(N^{2}\)) and Lamb
(\(S_{\ell}^{2}\)) frequencies versus fractional radius
\(x\). The horizontal dashed line shows the frequency squared
\(\omega^{2}\) of the \(\ell=1\), \(n_{\rm pg}=-7\)
mode shown in Fig. 10. Regions where
\(\omega^{2}\) is smaller (greater) than both \(N^{2}\) and
\(S_{\ell}^{2}\) are gravity (acoustic) propagation regions;
other regions are evanescent. (Source
)