wd_nova_burst

wd_nova_burst#

This test case checks the evolution of a nova outburst for one cycle.

This test case has 2 parts. Click to see a larger version of a plot.

  • Part 1 (inlist_setup) loads 1.1M_lgTc_7.7.mod, a prebuilt carbon-oxygen white dwarf from the make_o_ne_wd test suite, sets the topical depth to 30, evolves for a few steps to an age of 1e7 year and terminates.

  • Part 2 (inlist_wd_nova_burst) continues the evolution with the accretion of a hydrogen-helium mixture at a rate of 1e-9 Msun/yr. Eventually hydrogen burning in the accerted envelope causes the luminosity to exceed 1e4 Lsun and an alert is written to the terminat via the run_star_extras.f90:

                                            num_bursts           1
        405   8.277879  4.365E+05   8.624940   8.625037   1.115473   0.000025   0.000000   0.000000   0.309197 329.049014   1497     36
-5.0945E+00   7.679057  -1.661720  -2.464781   7.731723  -9.000000   1.115449   0.000000   0.400729   0.019380  65.113954     10
 2.4766E+04   7.926520   4.191187   4.972103  -6.862099  -4.289137   1.115442   0.038204   0.474628   1.000000  0.327E-01  max increase
                               rel_E_err   -1.9543755810854924D-08
                       log_rel_run_E_err       -6.0755826271591920

and when the luminosity falls below 1e3 Lsun a nova cycle is considered finished and the run terminates:

        585   7.835616  3.884E+05   1.624521   1.624523   1.115450   0.000001   0.000000   0.000000   0.332277 329.045159   1617     48
-1.4608E+00   7.679051  -2.147858 -16.859525   0.457865  -9.000000   1.115449   0.000000   0.400729   0.019454  65.113548      5
 2.4768E+04   7.926494   3.015873  -3.748867 -67.588316  -3.897373   1.115442   0.038204   0.474628   1.000000 -0.133E-06        Lnuc_H
                               rel_E_err   -5.3661942520223188D-10
                       log_rel_run_E_err       -5.9590177286527863

 have finished burst
../_images/wd_nova_burst_grid.svg

The Python script used to create the figures above:

import matplotlib.pyplot as plt
import numpy as np
import os
import mesa_reader as mr

# if the directory plt_out/ does not exits, make it
if not os.path.exists("plt_out"):
    os.makedirs("plt_out")

# "MESA" styles for plotting
# note that this requires having a LaTeX installation locally.
# If you don't have that, you can comment out the "text.usetex" line.
plt.rcParams["figure.figsize"] = (3.38, 2.535)
plt.rcParams["lines.markersize"] = 4
plt.rcParams["lines.linewidth"] = 1.5
plt.rcParams["text.usetex"] = True
plt.rcParams["font.size"] = 10
plt.rcParams["font.family"] = "serif"
plt.rcParams["font.serif"] = "Computer Modern Roman"
plt.rcParams["axes.titlesize"] = "medium"
plt.rcParams["axes.labelsize"] = "medium"
plt.rcParams["legend.fontsize"] = 8
plt.rcParams["legend.frameon"] = False
plt.rcParams["figure.dpi"] = 300

plt.rcParams["xtick.direction"] = "in"
plt.rcParams["ytick.direction"] = "in"
plt.rcParams["xtick.top"] = True
plt.rcParams["ytick.right"] = True
plt.rcParams["xtick.minor.visible"] = True
plt.rcParams["ytick.minor.visible"] = True

plt.rcParams["savefig.bbox"] = "tight"
plt.rcParams["savefig.pad_inches"] = 0.1
plt.rcParams["savefig.dpi"] = 300
plt.rcParams["savefig.format"] = "svg"

plt.rcParams["axes.formatter.use_mathtext"] = True

# read in the history data from MESA
logs = mr.MesaLogDir("LOGS")
h = logs.history
p = logs.profile_data()

# Create a figure
fig = plt.figure(figsize=(10, 6))

# Create a GridSpec with 6 rows and 3 columns
gs_outer = fig.add_gridspec(1, 2, width_ratios=[1, 2])
gs_left = gs_outer[0].subgridspec(2, 1, hspace=0.3)
gs_right = gs_outer[1].subgridspec(3, 2, hspace=0, wspace=0.75)

# First column: Two plots of roughly equal height
ax1 = fig.add_subplot(gs_left[0])
ax2 = fig.add_subplot(gs_left[1])

# Second column: Three stacked axes sharing the same x-axis
ax3 = fig.add_subplot(gs_right[0, 0])
ax4 = fig.add_subplot(gs_right[1, 0])
ax5 = fig.add_subplot(gs_right[2, 0])

# Third column: Three stacked axes sharing a different x-axis
ax6 = fig.add_subplot(gs_right[0, 1])
ax7 = fig.add_subplot(gs_right[1, 1])
ax8 = fig.add_subplot(gs_right[2, 1])

# HR diagram
ax1.loglog(h.Teff, h.L)
ax1.set_xlabel(r"Effective Temperature [K]")
ax1.set_ylabel(r"Luminosity [$L_{\odot}$]")
ax1.invert_xaxis()
# add marker to last point
ax1.plot(h.Teff[-1], h.L[-1], "o", color="C3")

# T-Rho Profile
# load data from MESA_DIR and save to arrays in T-Rho space
plot_info_dir = os.path.join(os.environ["MESA_DIR"], "data/star_data/plot_info")

# degeneracy data
psi4_file = os.path.join(plot_info_dir, "psi4.data")
psi4_data = np.genfromtxt(psi4_file)
psi4_xs = 10 ** np.array(psi4_data.T[0])
psi4_ys = 10 ** np.array(psi4_data.T[1])

# hydrogen burn line
hydrogen_burn_file = os.path.join(plot_info_dir, "hydrogen_burn.data")
hydrogen_burn_data = np.genfromtxt(hydrogen_burn_file)
hydrogen_burn_xs = 10 ** np.array(hydrogen_burn_data.T[0])
hydrogen_burn_ys = 10 ** np.array(hydrogen_burn_data.T[1])

# helium burn line
helium_burn_file = os.path.join(plot_info_dir, "helium_burn.data")
helium_burn_data = np.genfromtxt(helium_burn_file)
helium_burn_xs = 10 ** np.array(helium_burn_data.T[0])
helium_burn_ys = 10 ** np.array(helium_burn_data.T[1])

# plot raw T-Rho data; sets limits for us
ax2.loglog(p.Rho, p.T)

# plot psi4 and hydrogen burn lines
# preserve limits... shouldn't have to do it like this, but
# autoscale commands always seem to fail
x_left, x_right = ax2.get_xlim()
y_bottom, y_top = ax2.get_ylim()
ax2.plot(psi4_xs, psi4_ys, color="lightgrey", ls=":", zorder=-5)
ax2.plot(hydrogen_burn_xs, hydrogen_burn_ys, ls="--", color="lightgrey", zorder=-5)
ax2.plot(helium_burn_xs, helium_burn_ys, ls="--", color="lightgrey", zorder=-5)
ax2.set_xlim(x_left, x_right)
ax2.set_ylim(y_bottom, y_top)
ax2.set_xlabel(r"Density [g/cm$^3$]")
ax2.set_ylabel(r"Temperature [K]")

# plot strong burning zones
high_burning = np.where(p.eps_nuc > 1e7, True, False)
mid_burning = np.where((p.eps_nuc > 1e3) & (p.eps_nuc < 1e7), True, False)
low_burning = np.where((p.eps_nuc > 1) & (p.eps_nuc < 1e3), True, False)
# temporarily set cap style to round
save_capstyle = plt.rcParams["lines.solid_capstyle"]
plt.rcParams["lines.solid_capstyle"] = "round"
ax2.plot(
    p.Rho[high_burning],
    p.T[high_burning],
    marker="o",
    ls="",
    ms=6,
    color="C3",
    zorder=-1,
)
ax2.plot(
    [],
    [],
    lw=6,
    color="C3",
    label=r"$\epsilon_{\mathrm{nuc}} > 10^7\,\mathrm{erg/g/s}$",
)
ax2.plot(
    p.Rho[mid_burning],
    p.T[mid_burning],
    marker="o",
    ls="",
    ms=4.5,
    color="C1",
    zorder=-2,
)
ax2.plot(
    [],
    [],
    color="C1",
    lw=4.5,
    label=r"$\epsilon_{\mathrm{nuc}} > 10^3\,\mathrm{erg/g/s}$",
)
ax2.plot(
    p.Rho[low_burning],
    p.T[low_burning],
    marker="o",
    ls="",
    ms=3,
    color="goldenrod",
    zorder=-3,
)
ax2.plot(
    [],
    [],
    color="goldenrod",
    lw=3,
    label=r"$\epsilon_{\mathrm{nuc}} > 1\,\mathrm{erg/g/s}$",
)
ax2.legend(loc="lower right")
# restore capstyle
plt.rcParams["lines.solid_capstyle"] = save_capstyle

# time series for last 5 years
window = 5

# top panel: photospheric and H-burning luminosities
mask = h.star_age > (max(h.star_age) - window)
ax3.set_title("Last 5 Years")
ax3.semilogy(h.star_age[mask], h.L[mask], label=r"$L$")
ax3.semilogy(h.star_age[mask], h.LH[mask], ls="--", label=r"$L_{\mathrm{H}}$")
ax3.set_ylabel(r"Luminosity [$L_{\odot}$]")
ax3.legend(loc="upper left")
ax3.set_xticklabels([])

# temperature and radius
ax4.semilogy(h.star_age[mask], h.Teff[mask], label=r"$T_{\mathrm{eff}}$")
ax4.semilogy([], [], label=r"$R$", ls="--", color="C1")
ax4.legend(loc="upper left")
ax4b = ax4.twinx()
ax4b.semilogy(h.star_age[mask], h.R[mask], ls="--", color="C1")
ax4.set_ylabel(r"Eff. Temp. [K]", color="C0")
ax4b.set_ylabel(r"Radius [$R_{\odot}$]", color="C1")
ax4.set_xticklabels([])

# mass loss and envelope mass
ax5.semilogy(
    h.star_age[mask],
    h.star_mass[mask] - h.he_core_mass[mask],
    label=r"$\Delta M_{\mathrm{H}}$",
)
ax5.semilogy([], [], ls="--", label=r"$|\dot{M}|$")
ax5.legend(loc="center left")
ax5.set_ylabel(r"Envelope Mass [$M_{\odot}$]", color="C0")
ax5b = ax5.twinx()
ax5b.semilogy(h.star_age[mask], h.abs_mdot[mask], ls="--", color="C1")
ax5b.set_ylabel(r"$\left\vert\dot{M}\right\vert$ [$M_{\odot}$/yr]", color="C1")

ax5.set_xlabel(r"Star Age [yr]")

# profile at the end of the simulation
# first up, common isotope mass fractions
xm = p.star_mass - p.mass
for iso in ["h1", "he4", "c12", "n14", "o16"]:
    element = "".join([i for i in iso if not i.isdigit()]).capitalize()
    mass_number = "".join([i for i in iso if i.isdigit()])
    tex_iso = r"${}^{" + mass_number + r"}$" + element
    ax6.loglog(xm, p.data(iso), label=tex_iso)

ax6.set_title("Final Profile")
ax6.set_ylim(4e-5, 1.5)
ax6.legend(loc="lower left")
ax6.set_ylabel(r"Mass Fraction")
ax6.set_xticklabels([])

# nuclear energy generation rates
ax7.loglog(xm, p.eps_nuc, label=r"$\epsilon_{\mathrm{nuc}}$")
ax7.loglog(xm, p.pp, ls="--", label=r"$\epsilon_{\mathrm{pp}}$")
ax7.loglog(xm, p.cno, ls=":", label=r"$\epsilon_{\mathrm{CNO}}$")
ax7.legend(loc="upper right")
ax7.set_ylabel(r"Specific Power [erg/g/s]")
ax7.set_xticklabels([])
ax7.set_ylim(1e-6, 1e10)

# basic thermodynamic quantities
ax8.loglog(xm, p.Rho, label="$\\rho$")
ax8.loglog([], [], ls="--", label="$T$")
ax8.set_ylim(2e-3, 2e6)
ax8.set_ylabel(r"Density [g/cm$^3$]", color="C0")
ax8b = ax8.twinx()
ax8b.loglog(xm, p.T, ls="--", color="C1")
ax8b.set_ylim(2e6, 1.25e8)
ax8b.set_ylabel(r"Temperature [K]", color="C1")
ax8.legend(loc="upper right")

for ax in [ax6, ax7, ax8]:
    ax.set_xlim(8e-4, 1e-10)
ax8.set_xlabel(r"Exterior Mass Coordinate [M$_{\odot}$]")


fig.savefig("plt_out/wd_nova_burst_grid.svg")

Last-Updated: 07Aug2024 (MESA 0b40398b) by wmwolf.

Last-Run: 22Oct2024 (MESA 9b2017ca) by pmocz on C916PXT6XW in 603 seconds using 8 threads.