Evolve a population of stars
You can find this application in the demos folder of your Jupyter notebook environment.
- inlist
- inlist_create_zams
- inlist_hcoreburning
- multi_star_evolution.ipynb
In this tutorial, use Camber to generate some Modules for Experiments in Stellar Astrophysics (MESA) Stellar Tracks. Watch a population of stars evolve from the pre-main-sequence (pMS) to the zero-age-main-sequence (ZAMS), and then evolve them once more to the end of their hydrogen core-burning phase of evolution.
To execute the workload efficiently, this tutorial uses the create_scatter_job to execute many simulations in parallel.
In short, the flow is as follows:
- Initialize a set of inlist files, essentially a configuration file that defines the parameters to evolve one star.
- Execute one job for each of these inlists in parallel.
- Plot the output of each evolution in one chart.
Set up your MESA environment
First, import Camber and the libraries to do calculation and visualization:
import camber
import os
import pandas as pd
import numpy as np
import smplotlib
import matplotlib.pyplot as plt
%matplotlib inline
import mesa_reader as mrThen, define a function that creates an ensemble of stellar tracks, ranging in mass from minMass to maxMass with an interval size dm.
def create_table(minMass, maxMass, dm):
# Calculate the number of intervals
num_intervals = int((maxMass - minMass) / dm) + 1
# Generate the values based on the intervals
values = [minMass + i * dm for i in range(num_intervals)]
# Generate ZAMS file names
zams_file = [f"{round(val, 1)}Msun_zams.mod" for val in values]
hcore_file = [f"{round(val, 1)}Msun_hcore.mod" for val in values]
#outdir= [f"{val}Msun" for val in values]
# Create a DataFrame
df = pd.DataFrame({'mass': values, 'zams_file':zams_file,'hcore_file':hcore_file})
return dfminMass=1.0 ##
maxMass=4.5
dm=0.5
df = create_table(minMass,maxMass,dm)
dfGenerate inlist files
Now generate ZAMS stars from pMS collapse. To generate a ZAMS star, specify an inlist that includes the function create_pre_main_sequence_model = .true.
This inlist is called inlist_create_zams.
Note that the param_sets variable is going to be passed as an agrument to the MESA create_scatter_job method.
These params are the template for the jobs to execute in parallel.
jobs = camber.mesa.create_scatter_job(
param_sets=[
{
"INLIST": "'inlist_create_zams'",
"OUTPUT_MOD": '"' + str(row['zams_file']) + '"', # Add extra quotes around the value
"INIT_MASS": str(row['mass'])
}
for _, row in df.iterrows()
],
inlist_files = ['inlist','inlist_create_zams'],
engine_size='SMALL'
)This may take a few minutes on a small engine. To view the job status, you can return the jobs object:
jobsOnce all statuses are COMPLETED, you’ve generated your ZAMS files! Camber creates a final ZAMS model, appended with zams.mod, in each of the respective directories.
These ZAMS mod files are the requisite input files to evolve our stars from ZAMS to the next stage of evolution (in our case, H-Core burning).
Let’s move these files to our main directory:
!find $(pwd) -maxdepth 1 -type d -name "00*" -exec sh -c 'cp "$1"/*.mod $(pwd)' sh {} \;Now evolve each of our ZAMS stellar tracks along the H-core burning sequence.
To parameterize this, define another table of param_sets:
param_sets=[
{
"INLIST": "'inlist_hcoreburning'",
"INPUT_MOD": '"' + str(row['zams_file']) + '"', # Add extra quotes around the value
"INIT_MASS": str(row['mass']),
"OUTPUT_MOD": '"' + str(row['hcore_file']) + '"', # Add extra quotes around the value
}
for _, row in df.iterrows()]
param_setszams_mod_files = [file for file in os.listdir("./") if file.endswith("zams.mod")]
zams_mod_files=sorted(zams_mod_files)Run another scatter_job with the H-Core burning inlist files:
jobs = camber.mesa.create_scatter_job(
param_sets=[
{ "INLIST": "'inlist_hcoreburning'",
"INPUT_MOD":str(row['zams_file']),
"OUTPUT_MOD": str(row['hcore_file']),
"INIT_MASS": str(row['mass'])
}
for _, row in df.iterrows()
],
inlist_files = ['inlist','inlist_hcoreburning'],
model_files=zams_mod_files,
engine_size='SMALL'
)This too takes just a few minutes to complete, but you can start plotting as soon as your files begin to be generated. Check the status with:
for i in range(len(jobs)):
print(jobs[i].job_id, ": ", jobs[i].status)Once the statuses are RUNNING, you can start to plot results.
Plot results
fig = plt.figure(figsize=(4,6))
maxx=4.4
minx=3.5
miny=-0.5
maxy=4.5
dirs = sorted([file for file in os.listdir("./") if file.startswith("00")])
for d in dirs:
h=mr.MesaData(str(d)+'/LOGS/history.data')
plt.plot(h.log_Teff,h.log_L,linewidth=0.75,zorder=1,label=str(round(h.star_mass[0],4))+' $M_{\odot}$',color='black')
plt.annotate(str(round(h.star_mass[0],4))+' $M_{\odot}$', (max(h.log_Teff)+0.12,h.log_L[0]),fontsize=8)
plt.xlim(maxx,minx)
plt.ylim(miny,maxy)
plt.grid(alpha=0.25)
plt.xlabel('$\log$(Teff)')
plt.ylabel('$\log$(L)')
plt.show()
It’s always a good practice to clean-up after yourself. Go ahead and clear the jobs:
# Uncomment when you're sure you want to stop the jobs
# for job in jobs:
# job.clear()Next steps
What would the plot look like if you changed the value of minMass to 4 and the value of maxMass to 5?