Embedded Files
I. Setting up the folder
I. Setting up the folder
In this Minilab 2 we will evolve a star from the end of core helium depletion (this was the stopping condition of your Minilab 1), i.e. when the central helium abundance is below 0.0005, to central carbon depletion, when the carbon central abundance drops below 0.001. This stage is quite complicated to model, as the star is now more evolved, reaching the tip of its Giant Branch, and the model is evolving at smaller and smaller timescales dictated by the fast burning cycles.
The first task of this Minilab 2 will be to compute a quantity of interest, the Kelvin-Helmholtz timescale of tKH , within MESA, and see how different choices of αMLT can affect this quantity. The calculation of a custom quantity can be done by modifying the run_star_extras.f90 routines in the source folder ./src. Additionally, we will modify the inlist_minilab2.
SOLUTIONS...? The complete inlists with fully coded solutions are provided in the folder MINILAB2_alpha_sol. Try to read first the provided hints and give it a go, but if you're stuck at any point you can always look at those.
Move to the Minilab 2 folder
Move to the Minilab 2 folder
Inside the Google Drive folder of the school, navigate into the folder of our session, MINILAB2_alpha. Download it into your working folder and inspect inlist_minilab2 with your favorite text editor.
You should see a bunch of well-commented options, all within the &star_jobs and &controls. These are the only ones that you will be touching for this Minilab 2. Pay special attention to the control that sets the αMLT parameter:
&controls
! The mixing length is this parameter times a local pressure scale height.
! Default is 2.0.
mixing_length_alpha = 2.0
...
/ ! end of controls namelist
As you can see, we are going to start with αMLT=2, a value that has been adopted historically from calibrations with our Sun which also happens to be the one that MESA picks by default if mixing_length_alpha is unspecified.
BEWARE: Many massive star models use αMLT=1.5, chosen due to agreement with (among other things) chemical yields from single-star models of SN-1987A. This choice directly impacts the model's stellar radius, as we will see in this lab.
Pick your mass
Pick your mass
Go to the Google spreadsheet of Minilab 2, Day2_Massive_lab2, and write your name close to your favorite mass, which we will refer to as M in M☉. Make sure to have something different than your neighbor(s)! Once your mass is picked, you can make use of the bank of models folder called PRELAB_models, which is in the Google Drive folder of the school at the same level as the MINILAB2_alpha folder. The starting model of interest for your simulation is going to be located into a path that follows this notation:
[M]Msun/alpha_2.0/after_core_he_burn_[M]Msun_alpha_2.0.mod.
Find your model, download it and copy it into your working folder:
# Substitute the picked mass M in place of [M]
$ cp ../PRELAB_models/[M]Msun/alpha_2.0/after_core_he_burn_[M]Msun_alpha_2.0.mod .
Don't forget to put your selected mass inside the inlist_minilab2, specifically in the controls namelist, as follows:
&controls
! Initial mass in Msun units.
! When loading a saved model, it needs to correspond to the model's mass.
initial_mass = [M]
...
/ ! end of controls namelist
HELPFUL note: You might notice that in the PRELAB_models folder you also have, per each M and αMLT, a folder called png. These are storing all the pgstar window pngs of the runs that evolve the starting model after_core_he_burn_[M]Msun_alpha_[alpha].mod all the way to carbon depletion, i.e. when the central carbon abundance drops under 0.001. In a sense, they contain the pngs of the solution models of this Minilab 2 for each parameter set!! Therefore, if you are curious (or you want to fill in more parameters in the Google spreadsheet of Minilab 2 for the collaborative task, Day2_Massive_lab2), you can give a look at those as well.
Load the starting model
Load the starting model
Similarly to what you did in Minilab 1, we are going to evolve our star from a starting model. In this case, the model is saved approximately towards the end of helium burning in the core. As you did before, you are going to need to load the saved model by modifying inlist_minilab2:
&star_job
...
! Instructions to use the specified model as starting point.
load_saved_model = .true.
load_model_filename = 'after_core_he_burn_[M]Msun_alpha_2.0.mod'
...
/ ! end of star_job namelist
We strongly suggest you also modify the name of the model saved after termination, save_model_filename, and the name of the output folder log_directory, to something that resembles the box below. We are going to be modifying αMLT soon, so you do not want to accidentally overwrite anything.
save_model_filename = 'after_core_c_burn_[M]Msun_alpha_2.0.mod'
...
log_directory = 'LOGS_[M]Msun_alpha_2.0'
...
Grid2_file_dir = 'png_[M]Msun_alpha_2.0'
The above will save a final model at the end of your simulation; the termination condition for your setup is set to be when the central carbon abundance drops under 0.001, and you can see it in the &controls namelist of inlist_to_end_core_c_burn.
Just as a check, try to compile and start running the model and see if the right LOGS_[M]Msun_alpha_2.0 and png_[M]Msun_alpha_2.0 appear. If everything is going correctly, you should only see a message complaining about not finding a quantity, i.e. failed_to_find_value t_KH. Kill the run after a few timesteps (CTRL+C), don't waste too much time on this, we will solve this issue in a minute.
# Give permissions to run
$ chmod a+x rn clean mk re rn1
# Clean and compile
$ ./clean && ./mk
# Run
$ ./rn
...
# Wait few timesteps and kill
(CTRL+C)
II. Modify the run_star_extras.f90
II. Modify the run_star_extras.f90
Now it's time to modify the ./src/run_star_extras.f90 to compute the Kelvin-Helmholtz timescale tKH. This is a global property of your star, since it depends on its mass M, the radius R and surface luminosity L:
As you saw this morning, in order to be able to save a wanted quantity that is a global property of the star, you can simply instruct MESA to save a custom, additional column in your history.data output file.
Find your spot
Find your spot
How do you instruct to produce an additional history column? The right spot is inside a routine called data_for_extra_history_columns. This routine, after declaring relevant input and output, is initializing the star_info structure in the pointer s. It returns two 1-dimensional arrays, called names and vals, in which you can store a string with the name of the custom variable and its value, respectively. You should be seeing something like this:
subroutine data_for_extra_history_columns(id, n, names, vals, ierr)
...
character (len=maxlen_history_column_name) :: names(n)
real(dp) :: vals(n)
...
! TO DO
! Kelvin-Helmholtz timescale computation, see Eq.(2) of
! your instructions.
names(1) = "t_KH"
vals(1) = 0.d0
end subroutine data_for_extra_history_columns
Whenever you see a ! TO DO line, that is your spot to work on. You see that the definitions of the names and vals arrays are such that you can store n names of variables, and per each n index you can store the value of the variable throughout the entire evolutionary history. The value of n is set by the routine called how_many_extra_history_columns, which you will modify soon.
Compute the Kelvin-Helmholtz timescale
Compute the Kelvin-Helmholtz timescale
You see that now the first element of the vals array is zero, vals(1) = 0d0. If you remember how to access the relevant information from the star_info structure with the s pointer, you can try and put the equation for tKH. Otherwise, you give a look at these hints here below, and / or at the solutions. Remember that you can find the parameter names for the mass M, radius R, and luminosity L inside the files contained in the directory: $MESA_DIR/star_data/public/. In this particular case, you want to have a look at the file $MESA_DIR/star_data/public/star_data_step_work.inc.
Hint 1: Global properties?
Your star is divided in a number nz of zones, which defines the spatial resolution of the simulation. Each zone has an index that can range from 1 to nz, and your surface corresponds to the index 1 and the center to nz. You want to use the values of the mass, radius and luminosity at the index of the surface, because those represent the integrated values from the center to the surface, thus the global (rather than the local!) properties.
Hint 2: The star pointer!
In MESA, all of the star's structure variables are contained in the star pointer, initialized in the line call star_ptr(id, s, ierr). This pointer is like a class in other programming languages, and the attributes of the pointer are called for some quantity at index i like s% quantity(i).
Hint 3: Mind the units!
In stellar astrophysics and in MESA we like to use the centimeter-gram-second units, therefore we saved our own useful constants to convert between energy units, or from seconds to years, and such. Those constants are readily accessible in run_star_extras.f90, if you know what their name is 🙃 You may want to use: secyer, which is the conversion factor between years and seconds (secyer = 365.25d0*24*60*60) and standard_cgrav, the gravitational constant G in c.g.s. units (cm3 s-1 g-2). These are loaded into run_star_extras via the declaration use const_def; for more definitions, see $MESA_DIR/const/public/const_def.f90.
How many columns?
How many columns?
Now that you implemented the tKH equation and you instructed the data_for_extra_history_columns routine to save the information on your history.data file, you still need to tell MESA that there is the need to fit this additional column into the output. Find the function called how_many_extra_history_columns and modify it accordingly:
integer function how_many_extra_history_columns(id)
...
! TO DO
! Add a column to the count
how_many_extra_history_columns = 0
end function how_many_extra_history_columns
Solution
integer function how_many_extra_history_columns(id)
...
! TO DO
! Add a column to the count
how_many_extra_history_columns = 1
end function how_many_extra_history_columns
subroutine data_for_extra_history_columns(id, n, names, vals, ierr)
...
! TO DO
! Kelvin-Helmholtz timescale computation, see Eq.(2) of
! your instructions.
names(1) = "t_KH"
vals(1) = 0.5*standard_cgrav*(s% m(1))*(s% m(1))/s% L(1)/s% r(1)/secyer
end subroutine data_for_extra_history_columns
III. The fiducial mixing length theory αMLT
III. The fiducial mixing length theory αMLT
Run the fiducial αMLT simulation
Run the fiducial αMLT simulation
You are all set to run your first simulation! We will first need to recompile MESA since, by playing with run_star_extras.f90, you modified the source fortran code. Let's clean and compile:
# Clean and compile
$ ./clean && ./mk
If the compilation is successful, you can go ahead and run the simulation with ./rn. Notice that we did not change any physics, especially the mixing_length_alpha = 2.0 value. We first want to see what happens when αMLT is posed to the fiducial value.
During the run, you can give a look at the pgstar window, which should look something like the picture here below. If you did everything correctly till this point, you should also be able to see the evolution of the Kelvin-Helmholtz timescale log tKH (yr) as a function of your star's age.
Look at the plot of log tKH as a function of model_number.
Notice how your star's radius log R is increasing in the Red Giant Branch, and how, consequently, the Kelvin-Helmholtz timescale tKH drops by few dex. Towards the end of the run, the behavior of tKH mirrors exactly that of the radius, as the star is basically sitting in the same spot of the HR diagram as it is burning carbon very fast in its core.
In the Profile panel all the way to the right in the upper corner, you can also notice the large convective envelope that your star is developing (see also the Kippenhahn diagram in the bottom left corner, with all the light blue present at the higher mass coordinates): we will think about it a bit more later.
NOTE: If your pgstar window is too big (or too small) for your screen, you can open inlist_pgstar_minilab2 and modify the following two options
&pgstar
...
! PGSTAR WINDOW DIMENSIONS
Grid2_win_width = 15
Grid2_win_aspect_ratio = 0.65 ! aspect_ratio = height/width
...
/ ! end of pgstar namelist
Example pgstar window for 12 M☉ model at the end of core Helium burning.
Crowdsourcing
Crowdsourcing
When the run is finished, the pgstar window will freeze on your screen, showing you the last model of the run. You should see something like:
termination code: xa_central_lower_limit
pause_before_terminate: hit RETURN to continue
This is caused by the line pause_before_terminate = .true. in &star_job. You may wait a bit before pressing your return button: stare at the final pgstar window.
Annotate, in the Google spreadsheet of Minilab 2, the values of the final surface luminosity log L (L☉), radius log R (R☉), effective temperature log Teff (K) and Kelvin-Helmholtz timescale tKH (yr). You can do this in three ways: if you have very good eyes, you can inspect the HR diagram to estimate log L (L☉) and log Teff (K), and the History panel plots to estimate log R (R☉) and tKH (yr). For a more accurate read, you can also read the summary text at the bottom of the pgstar window, in which you should find all you need (tKH will be in years). Otherwise, you can find these quantities in the last rows of your terminal output, see example down here (tKH is in years there too).
Great! We've calculated the thermal timescale globally, and run with a fiducial value of alpha_MLT! In the next Task, we will change 2 assumptions and explore the implications: (1) we will change our assumed alpha_MLT and compare the structure of the two different models, and (2) we will calculate the thermal timescale locally rather than globally and compare the two timescales within the same model!
Please navigate to Minilab 2 Task 2: The local thermal timescale !
Page updated
Google Sites
Report abuse