1.3M_ms_high_Z
The test checks the evolution of metal-rich low-mass stars by evolving a 1.3 \({\rm M}_\odot\), metal-rich Z=0.04 model from the pre-main sequence to core hydrogen depletion.
This test case has two parts.
Part 1 (
inlist_zams) creates the pre-main-sequence model and stops near zams.Part 2 (
inlist_1.3M_ms_high_Z) evolves the model from zams to when the luminosity reaches log10(L/Lsun) = 0.7.
The final model, click on the image for a larger version, shows
The left plot shows the HR diagram. The yellow curve is a precalculated
HR track loaded from HR_OPAL.dat, while the green curve is the model calculated.
The right plot shows a profile of the metal mass fraction. The curve rises above the
background metallicity as hydrogen burns to helium and goes below the background metallicity
as core hydrogen depletes.
pgstar commands used for the plots above:
&pgstar
! MESA uses PGPLOT for live plotting and gives the user a tremendous
! amount of control of the presentation of the information.
! device
file_white_on_black_flag = .true.
file_device = 'vcps' ! postscript
file_extension = 'ps'
!file_device = 'png' ! png
!file_extension = 'png'
pgstar_interval = 10
Grid1_title = '1.3M_ms_high_Z'
Grid1_win_flag = .false.
Grid1_win_width = 15
Grid1_win_aspect_ratio = 0.4 ! aspect_ratio = height/width
Grid1_xleft = 0.05 ! fraction of full window width for margin on left
Grid1_xright = 0.95 ! fraction of full window width for margin on right
Grid1_ybot = 0.08 ! fraction of full window width for margin on bottom
Grid1_ytop = 0.92 ! fraction of full window width for margin on top
Grid1_num_cols = 5 ! divide plotting region into this many equal width cols
Grid1_num_rows = 1 ! divide plotting region into this many equal height rows
Grid1_num_plots = 2 ! <= 10
Grid1_plot_name(1) = 'HR'
Grid1_plot_row(1) = 1 ! number from 1 at top
Grid1_plot_rowspan(1) = 1 ! plot spans this number of rows
Grid1_plot_col(1) = 1 ! number from 1 at left
Grid1_plot_colspan(1) = 2 ! plot spans this number of columns
Grid1_plot_pad_left(1) = 0.025 ! fraction of full window width for padding on left
Grid1_plot_pad_right(1) = 0.05 ! fraction of full window width for padding on right
Grid1_plot_pad_top(1) = 0.05 ! fraction of full window height for padding at top
Grid1_plot_pad_bot(1) = 0.05 ! fraction of full window height for padding at bottom
Grid1_txt_scale_factor(1) = 1.0 ! multiply txt_scale for subplot by this
Grid1_plot_name(2) = 'Profile_Panels1'
Grid1_plot_row(2) = 1 ! number from 1 at top
Grid1_plot_rowspan(2) = 1 ! plot spans this number of rows
Grid1_plot_col(2) = 3 ! number from 1 at left
Grid1_plot_colspan(2) = 3 ! plot spans this number of columns
Grid1_plot_pad_left(2) = 0.05 ! fraction of full window width for padding on left
Grid1_plot_pad_right(2) = 0.025 ! fraction of full window width for padding on right
Grid1_plot_pad_top(2) = 0.05 ! fraction of full window height for padding at top
Grid1_plot_pad_bot(2) = 0.05 ! fraction of full window height for padding at bottom
Grid1_txt_scale_factor(2) = 1.0 ! multiply txt_scale for subplot by this
! file output
Grid1_file_flag = .true.
Grid1_file_dir = 'pgstar_out'
Grid1_file_prefix = 'grid1'
Grid1_file_interval = 1000000 ! output when mod(model_number,Grid1_file_interval)==0
Grid1_file_width = -1 ! negative means use same value as for window
Grid1_file_aspect_ratio = -1 ! negative means use same value as for window
! show HR diagram
! this plots the history of L,Teff over many timesteps
HR_win_flag = .false.
! set static plot bounds
HR_logT_min = 3.55
HR_logT_max = 3.85
HR_logL_min = 0.1
HR_logL_max = 1.0
! show OPAL results
HR_fname = 'HR_OPAL.dat' ! file name for extra HR data
! set window size (aspect_ratio = height/width)
HR_win_width = 6
HR_win_aspect_ratio = 1.0
Profile_Panels1_win_flag = .false.
Profile_Panels1_win_width = 6
Profile_Panels1_win_aspect_ratio = 0.75 ! aspect_ratio = height/width
Profile_Panels1_xleft = 0.15
Profile_Panels1_xright = 0.85
Profile_Panels1_ybot = 0.15
Profile_Panels1_ytop = 0.85
Profile_Panels1_txt_scale = 1.0
Profile_Panels1_title = 'Profile Panels1'
Profile_Panels1_xaxis_name = 'mass'
Profile_Panels1_xaxis_reversed = .false.
Profile_Panels1_xmin = -101d0 ! only used if /= -101d0
Profile_Panels1_xmax = -101d0 ! only used if /= -101d0
Profile_Panels1_xmargin = 0d0
Profile_Panels1_show_mix_regions_on_xaxis = .false.
Profile_Panels1_yaxis_name(:) = ''
Profile_Panels1_yaxis_reversed(:) = .false.
Profile_Panels1_yaxis_log(:) = .false. ! show log10 of abs value
Profile_Panels1_ymin(:) = -101d0 ! only used if /= -101d0
Profile_Panels1_ymax(:) = -101d0 ! only used if /= -101d0
Profile_Panels1_ymargin(:) = 0.1
Profile_Panels1_dymin(:) = -1
Profile_Panels1_other_yaxis_name(:) = ''
Profile_Panels1_other_yaxis_reversed(:) = .false.
Profile_Panels1_other_yaxis_log(:) = .false. ! show log10 of abs value
Profile_Panels1_other_ymin(:) = -101d0 ! only used if /= -101d0
Profile_Panels1_other_ymax(:) = -101d0 ! only used if /= -101d0
Profile_Panels1_other_ymargin(:) = 0.1
Profile_Panels1_other_dymin(:) = -1
Profile_Panels1_show_grid = .false.
! setup default plot
Profile_Panels1_num_panels = 1
Profile_Panels1_yaxis_name(1) = 'z_mass_fraction_metals'
Profile_Panels1_ymin(1) = 0.038 ! only used if /= -101d0
Profile_Panels1_ymax(1) = 0.042 ! only used if /= -101d0
! Profile_Panels1_other_yaxis_name(1) = 'entropy'
! Profile_Panels1_yaxis_name(2) = 'logRho'
! Profile_Panels1_other_yaxis_name(2) = 'logP'
! file output
Profile_Panels1_file_flag = .false.
Profile_Panels1_file_dir = 'pgstar_out'
Profile_Panels1_file_prefix = 'profile_panels1_'
Profile_Panels1_file_interval = 5 ! output when mod(model_number,Profile_Panels1_file_interval)==0
Profile_Panels1_file_width = -1 ! (inches) negative means use same value as for window
Profile_Panels1_file_aspect_ratio = -1 ! negative means use same value as for window
/ ! end of pgstar namelist
The figures can also be created in Python with matplotlib and reading in data with the mesa_reader package:
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
# pip install git+https://github.com/wmwolf/py_mesa_reader.git
# if the directory plt_out/ does not exits, make it
if not os.path.exists("plt_out"):
os.makedirs("plt_out")
# enable Latex for figure labels
plt.rcParams["text.usetex"] = True
# read in the history data from MESA
history = mr.MesaData("LOGS/history.data")
# load the last profile saved (largest model number)
load_dir = mr.MesaLogDir("./LOGS")
profile = load_dir.profile_data()
# load HR_OPAL.dat reference data
hr_opal = np.genfromtxt("HR_OPAL.dat", delimiter=" ", skip_header=1)
# Plot HR diagram
HR_logT_min = 3.55
HR_logT_max = 3.85
HR_logL_min = 0.1
HR_logL_max = 1.0
plt.figure()
plt.plot(hr_opal[:,0], hr_opal[:,1], label="OPAL.dat", color="tab:orange")
plt.plot(history.log_Teff, history.log_L, label="MESA", color="tab:blue")
plt.plot(history.log_Teff[-1], history.log_L[-1], "ro")
plt.xlabel(r"$\log_{10}(T_{\rm eff})$")
plt.ylabel(r"$\log_{10}(L/L_{\odot})$")
plt.title("HR Diagram")
plt.xlim(HR_logT_min, HR_logT_max)
plt.ylim(HR_logL_min, HR_logL_max)
plt.legend()
plt.gca().set_box_aspect(1)
plt.gca().invert_xaxis()
plt.savefig("plt_out/HR_diagram.svg", bbox_inches="tight", pad_inches=0)
# plt.show()
# Plot Profile of metal mass fraction
plt.figure()
plt.plot(profile.mass, profile.z_mass_fraction_metals)
plt.xlabel(r"$\log_{10}(M(r)/M_{\odot})$")
plt.ylabel(r"metal mass fraction")
plt.title("Metal mass fraction profile")
plt.xlim(0, profile.mass[0])
plt.ylim(0.038, 0.042)
plt.savefig("plt_out/z_mass_fraction_metals.svg", bbox_inches="tight", pad_inches=0)
# plt.show()
Last-Updated: 27May2021 (MESA ebecc10) by fxt
Last-Run: 26Jul2024 (MESA afd04d7d) by pmocz on C916PXT6XW in 159 seconds using 8 threads.