controls

Contents

controls#

specifications for starting model#

NOTE: if you are loading a saved model, then the following initial values are NOT USED to modify the model. in particular, you cannot use these to change Y or Z of an existing model. if you want to do that, see star_job.defaults controls such as change_Y. however, these are reported in output as the initial values for the star.

initial_mass#

initial mass in Msun units. can be any value you’d like when you are creating a pre-main sequence model.

NOTE: this is not used when loading a saved model. however is reported in output as the initial mass of the star. don’t let that confuse you.

if you are loading a ZAMS model and the requested mass is in the range of prebuilt models, the code will interpolate in mass using the closest prebuilt models. if the requested mass is beyond the range of the prebuilt models, the code will load the closest one and then call “relax mass” to create a model to match the request. the prebuilt range is 0.08 Msun to 100 Msun, so the relax_mass method is only used for extreme cases. there are enough prebuilt models that the interpolation in mass seems to work fine for many applications.

initial_mass = 1

initial_z#

initial metallicity for create pre-ms and create initial model initial_z can be any value from 0 to 0.04

not used when loading a saved model. however is reported in output as the initial Z of the star.

however, if you are loading a zams model, then initial_z must match one of the prebuilt values. look in the 'data/star_data/zams_models' directory to see what prebuilt zams Z’s are available. at time of writing, only 0.02 was included in the standard version of star.

initial_z = 0.02d0

initial_y#

initial helium mass fraction for create pre-ms and create initial (< 0 means use default which is 0.24 + 2*initial_z)

not used when loading a saved model or a zams model. however is reported in output as the initial Y of the star.

NOTE: this is only used for create pre-main-sequence model and create initial model, and not when loading a zams model.

initial_y = -1

initial_he3#

initial mass fraction of he3. Must be smaller than initial_y. (< 0 means use default which is taken as solar scaled such that he4/he3 has the same value as the Sun)

not used when loading a saved model or a zams model.

initial_he3 = -1

controls for output#

terminal_interval#

write info to terminal when mod(model_number, terminal_interval) = 0. note: this replaces the obsolete control terminal_cnt.

terminal_interval = 1

write_header_frequency#

output the log header info to the terminal when mod(model_number, write_header_frequency*terminal_interval) = 0.

write_header_frequency = 10

extra_terminal_output_file#

if not empty, output terminal info to this file in addition to terminal. this does not capture all of the terminal output – just the common items. it is intended for use in situations where you cannot directly see the terminal output such as when running on a cluster. if you want to be able to monitor the progress for such cases, you can set extra_terminal_output_file = 'log' and then do tail -f log to view the terminal output as it is recorded in the file.

extra_terminal_output_file = ''

terminal_show_age_units#

terminal_show_timestep_units#

terminal_show_log_dt#

terminal_show_log_age#

options are ‘years’, ‘yrs’, ‘days’, ‘secs’, or ‘seconds’ this replaces the old controls terminal_show_age_in_years & terminal_show_age_in_days

terminal_show_age_units = 'years'
terminal_show_timestep_units = 'years'
terminal_show_log_dt = .true.
terminal_show_log_age = .false.

num_trace_history_values#

any valid name for a history data column, such as surf_v_rot for example if you have rapid rotation at the surface, you might want to try something like this:

num_trace_history_values = 7
trace_history_value_name(1) = 'surf_v_rot'
trace_history_value_name(2) = 'surf_omega_div_omega_crit'
trace_history_value_name(3) = 'log_rotational_mdot_boost'
trace_history_value_name(4) = 'log_total_angular_momentum'
trace_history_value_name(5) = 'center n14'
trace_history_value_name(6) = 'surface n14'
trace_history_value_name(7) = 'average n14'

value must be less than or equal to 10

num_trace_history_values = 0

trace_history_value_name(:)#

write values to terminal

trace_history_value_name(:) = ''

photo_directory#

directory for binary snapshots used in restarts

photo_directory = 'photos'

photo_interval#

save a photo file for possible restarting when mod(model_number, photo_interval) = 0. note: this replaces the obsolete control photostep.

photo_interval = 50

photo_digits#

use this many digits from the end of the model_number for the photo name

photo_digits = 3

log_directory#

directory for data files about the run

log_directory = 'LOGS'

do_history_file#

a history file is created if this is true

do_history_file = .true.

history_interval#

append an entry to the history.data file when mod(model_number, history_interval) = 0.

history_interval = 5

star_history_name#

name of history file

star_history_name = 'history.data'

star_history_header_name#

If not empty, then put star history header info in star_history_name file. In this case the history file has only data, making it easier to use with some plotting packages.

star_history_header_name = ''

star_history_dbl_format#

format for writing reals to star_history_name file

star_history_dbl_format = '(1pes40.16e3, 1x)'

star_history_int_format#

format for writing integer to star_history_name file

star_history_int_format = '(i40, 1x)'

star_history_txt_format#

format for writing characters to star_history_name file

star_history_txt_format = '(a40, 1x)'

write_profiles_flag#

profiles are written only if this is true

write_profiles_flag = .true.

profile_interval#

save a model profile info when mod(model_number, profile_interval) = 0.

profile_interval = 50

priority_profile_interval#

give saved profile a higher priority for retention when mod(model_number, priority_profile_interval) = 0.

priority_profile_interval = 1000

profiles_index_name#

name of the profile index file

profiles_index_name = 'profiles.index'

profiles_data_prefix#

prefix of the profile data

profile_data_prefix = 'profile'

profiles_data_suffix#

suffix of the profile data

profile_data_suffix = '.data'

profile_data_header_suffix#

If not empty, then put profile data header info here. In this case the profile data file has only data, making it easier to use with some plotting packages.

profile_data_header_suffix = ''

profile_dbl_format#

format for writing reals to profile file

profile_dbl_format = '(1pes40.16e3, 1x)'

profile_int_format#

format for writing integers to profile file

profile_int_format = '(i40, 1x)'

profile_txt_format#

format for writing characters to profile file

profile_txt_format = '(a40, 1x)'

max_num_profile_zones#

if nz > this, then only write a subsample of the zones. only used if > 1

max_num_profile_zones = -1

max_num_profile_models#

Maximum number of saved profiles. If there’s no limit on the number of profiles saved, you can fill up your disk – I’ve done it. So it’s a good idea to set this limit to a reasonable number such as 20 or 30. Once that many have been saved during a run, old ones will be discarded to make room for new ones. Profiles that were saved for key events are given priority and aren’t removed as long as there is a lower priority profile that can be discarded instead. Less than zero means no limit.

max_num_profile_models = 100

profile_model#

save profile when model_number equals this

profile_model = -1

profile_header_include_sys_details#

if this is true, profile header includes strings for compiler, build, etc.

profile_header_include_sys_details = .true.

write_model_with_profile#

if this is true, models are written at same time as profiles

write_model_with_profile = .false.

model_data_prefix#

prefix of the model data files

model_data_prefix = 'profile'

model_data_suffix#

suffix of the model data files

model_data_suffix = '.mod'

write_controls_info_with_profile#

if this is true, the values of the options in the controls inlist are written at same time as profiles

write_controls_info_with_profile = .false.

controls_data_prefix#

prefix of the control data files

controls_data_prefix = 'controls'

controls_data_suffix#

suffix of the control data files

controls_data_suffix = '.data'

mixing_D_limit_for_log#

if max D_mix in mixing region is less than this, don’t include the region in the log doesn’t apply to thermohaline or semiconvective regions

mixing_D_limit_for_log = 1d4

write_pulse_data_with_profile#

If .true., also save model data in a format compatible with pulsation codes.

write_pulse_data_with_profile = .false.

pulse_data_format#

Format of pulsation data. Current options (case insensitive) are:

'CAFEIN'

Format defined for CAFein (Valsecchi et al. 2013).

'FGONG'

Format originally defined for the GONG solar model project. A definition was given in 2005 for the CoRoT/ESTA project and GONG itself. MESA’s implementation largely follows this subsequent 2015 definition.

'OSC'

Format similar to FGONG originally produced for models from CESAM. Also defined for the CoRoT/ESTA project.

'GYRE'

The plain-text format defined for GYRE, which GYRE itself refers to as 'MESA' format.

'GSM'

The HDF5-based format defined for GYRE.

'SAIO'

Format for Saio’s pulsation code.

'GR1D'

Format for GR1D, defined in Sec. 3 of the GR1D documentation.

pulse_data_format = 'FGONG'

add_atmosphere_to_pulse_data#

if true, write atmosphere to pulse files. This is not valid when atm_option = 'table'

add_atmosphere_to_pulse_data = .false.

add_center_point_to_pulse_data#

if true, add point for r=0 to pulse files

add_center_point_to_pulse_data = .true.

keep_surface_point_for_pulse_data#

if true, add k=1 cell to pulse files

keep_surface_point_for_pulse_data = .false.

add_double_points_to_pulse_data#

add double points at discontinuities

add_double_points_to_pulse_data = .false.

interpolate_rho_for_pulse_data#

If true, then get rho_face by interpolating rho at cell center. If false, then calculate rho_face by dm/(4*pi*r^2*dr).

interpolate_rho_for_pulse_data = .true.

threshold_grad_mu_for_double_point#

threshold in grad_mu = dln(mu)/dln(P) for a double point to be written

threshold_grad_mu_for_double_point = 10d0

max_number_of_double_points#

maximum number of double points to be written (0 = no limit); when this limit is set, double points are chosen in order of decreasing \(|\nabla_\mu|\)

max_number_of_double_points = 0

fgong_header#

These are the four lines of arbitrary text that appear at the beginning of an FGONG file to describe its contents.

fgong_header(1) = 'FGONG file'
fgong_header(2) = 'Created by MESAstar'
fgong_header(3) = ''
fgong_header(4) = ''

fgong_ivers#

The version number for the FGONG file, which can only be 300 or 1300. 300 gives the old narrow format '(1P5E16.9,x)', which can produce numbers with no space separating them. 1300 gives the ‘wide’ FGONG format '(1P,5(X,E26.18E3))', as agreed on at the 5th Aarhus RGB workshop (University of Birmingham, UK, October 2015).

fgong_ivers = 1300 ! 300 or 1300 only

format_for_OSC_data#

float format for 'OSC' data format

format_for_OSC_data = '(1P5E19.12,x)'

gyre_data_schema#

data schema to use when passing model internally to GYRE, and/or writing files in GYRE/GSM format. For instance, to write v1.20 GYRE files use the value 120

gyre_data_schema = 101

max_num_gyre_points#

limit gyre output files to at most this number of points only used when > 1

max_num_gyre_points = -1

fgong_zero_A_inside_r#

when writing FGONG, if r < this and cell has mixing of some kind, force A = 0 Rsun units

fgong_zero_A_inside_r = 0d0

trace_mass_location#

location for trace_mass_radius, trace_mass_logT, etc. (Msun units)

trace_mass_location = 0

min_tau_for_max_abs_v_location#

controls choice of location in model for max_abs_v history info. can use this to exclude locations too close to surface. ignore if <= 0

min_tau_for_max_abs_v_location = 0

min_q_for_inner_mach1_location#

controls choice of location in model for innermost mach 1 history info. can use this to exclude locations too close to center.

min_q_for_inner_mach1_location = 0

max_q_for_outer_mach1_location#

controls choice of location in model for outermost mach 1 history info. can use this to exclude locations too close to surface.

max_q_for_outer_mach1_location = 1

burn_min1#

used for reporting where burning zone occur, for example in the pgstar TRho profiles. see star/public/star_data.inc for details. must be < burn_min2. In ergs/g/sec.

burn_min1 = 50

burn_min2#

used for reporting where burning zone occur, for example in the pgstar TRho profiles. see star/public/star_data.inc for details. In ergs/g/sec.

burn_min2 = 1000

width_for_limit_conv_vel#

look this number of cells on either side of boundary to see if any boundary k in that range has s% csound(k) < s% v(k) <= s% csound(k-1) i.e. transition from subsonic to supersonic as go inward if find any such transition then don’t allow increase in convection velocity. this implies no change from radiative to convective. the purpose of this is to prevent convective energy transport from moving energy from behind a shock to in front of the shock.

width_for_limit_conv_vel = 3

max_q_for_limit_conv_vel#

for q(k) <= this, don’t allow conv_vel to grow

max_q_for_limit_conv_vel = -1

max_r_in_cm_for_limit_conv_vel#

for r(k) <= this, don’t allow conv_vel to grow

max_r_in_cm_for_limit_conv_vel = -1

max_mass_in_gm_for_limit_conv_vel#

for m(k) <= this, don’t allow conv_vel to grow

max_mass_in_gm_for_limit_conv_vel = -1

center_avg_value_dq#

reported center values are averages over this fraction of star mass

center_avg_value_dq = 1d-8

surface_avg_abundance_dq#

reported surface abundances are averages over this fraction of star mass

surface_avg_abundance_dq = 1d-8

conv_core_gap_dq_limit#

skip non-convective gaps of less than this limit when reporting convective core size

conv_core_gap_dq_limit = 0d0

definition of core boundaries#

he_core_boundary_h1_fraction#

If fraction >= 0, boundary is outermost location where h1 mass fraction is <= this value, and he4 mass fraction >= min_boundary_fraction (see below). If fraction < 0, boundary is outermost location where he4 is the most abundant species.

he_core_boundary_h1_fraction = 0.1d0

co_core_boundary_he4_fraction#

If fraction >= 0, boundary is outermost location where he4 mass fraction is <= this value, and c12+o16 mass fraction >= min_boundary_fraction (see below). If fraction < 0, boundary is outermost location where c12+o16 is more abundant than any other species.

co_core_boundary_he4_fraction = 0.1d0

one_core_boundary_he4_c12_fraction#

If fraction >= 0, boundary is outermost location where combine he4+c12 mass fraction is <= this value, and o16+ne20 mass fraction >= min_boundary_fraction (see below). If fraction < 0, boundary is outermost location where o16+ne20 is more abundant than any other species.

one_core_boundary_he4_c12_fraction = 0.1d0

fe_core_boundary_si28_fraction#

For this case, “iron” includes any species with A > 46. If fraction >= 0, boundary is outermost location where si28 mass fraction is <= this value, and “iron” mass fraction >= min_boundary_fraction (see below). If fraction < 0, boundary is outermost location where “iron” is the most abundant species.

fe_core_boundary_si28_fraction = 0.1d0

neutron_rich_core_boundary_Ye_max#

Boundary is outermost location where Ye is <= this value.

neutron_rich_core_boundary_Ye_max = 0.48d0

min_boundary_fraction#

Value for deciding boundary regions.

min_boundary_fraction = 0.1d0

when to stop#

max_model_number#

The code will stop when it reaches this model number. Negative means no maximum.

max_model_number = -1

when_to_stop_rtol#

when_to_stop_atol#

Relative (rtol) and absolute (atol) error criteria for target stopping values. To compare how accurately the last step satisfied a stopping condition, MESA evaluates

error = abs(value - target_value)/ &
        (when_to_stop_atol + when_to_stop_rtol*max(abs(value), abs(target_value)))

and will redo with a smaller timestep if error is greater than 1. The default values 1d99 for both guarantee that error is tiny, so the run terminates as soon as a stopping condition is reached.

If you wish to use either rtol or atol, you should change the other to 0 (or your desired value). To see why, suppose we only set when_to_stop_rtol = 1d-3. Because when_to_stop_atol is still 1d99, error will still be very small and MESA won’t redo the last step.

when_to_stop_rtol = 1d99
when_to_stop_atol = 1d99

max_age#

max_age_in_days#

max_age_in_seconds#

Stop when the age of the star exceeds this value in years, days, or seconds. Only applies when > 0.

max_age = 1d36
max_age_in_days = -1
max_age_in_seconds = -1

num_adjusted_dt_steps_before_max_age#

This adjusts max_years_for_timestep so that hit max_age exactly, without needing possibly large change in timestep at end of run. only used if > 0

number of time steps to adjust to prior to hitting max age only used if > 0

num_adjusted_dt_steps_before_max_age = 0

dt_years_for_steps_before_max_age#

timestep in years

dt_years_for_steps_before_max_age = 1d6

reduction_factor_for_max_timestep#

per time step reduction limited to this

reduction_factor_for_max_timestep = 0.98d0

gamma_center_limit#

gamma is the plasma interaction parameter. Stop when the center value of gamma exceeds this limit.

gamma_center_limit = 1d99

eta_center_limit#

eta is the electron chemical potential in units of k*T. Stop when the center value of eta exceeds this limit.

eta_center_limit = 1d99

log_center_density_upper_limit#

log_center_density_lower_limit#

Stop when log10 of the center density is above/below the upper/lower limit.

log_center_density_upper_limit = 12d0
log_center_density_lower_limit = -1d99

log_center_temp_upper_limit#

log_center_temp_lower_limit#

Stop when log10 of the center temperature is above/below the upper/lower limit.

log_center_temp_upper_limit = 11d0
log_center_temp_lower_limit = -1d99

surface_accel_div_grav_limit = -1#

This is used when do not have a velocity variable. The acceleration ratio is abs(accel)/grav at surface, where accel is (rdot-rdot_old)/dt and grav is G*m/r^2. Stop if the ratio becomes larger than this limit. Ignored if <= 0.

surface_accel_div_grav_limit = -1

log_max_temp_upper_limit#

log_max_temp_lower_limit#

stop when log10 of the maximum temperature is above/below the upper/lower limit.

log_max_temp_upper_limit = 99
log_max_temp_lower_limit = -99

center_entropy_upper_limit#

center_entropy_lower_limit#

stop when the center entropy is above/below the upper/lower limit. in kerg per baryon, where kerg = (erg K^-1)

center_entropy_upper_limit = 1d99
center_entropy_lower_limit = -1d99

max_entropy_upper_limit#

max_entropy_lower_limit#

stop when the max entropy is above/below the upper/lower limit. in kerg per baryon, where kerg = (erg K^-1)

max_entropy_upper_limit = 1d99
max_entropy_lower_limit = -1d99

xa_central_lower_limit_species#

xa_central_lower_limit#

Lower limits on central mass fractions. Stop when central abundance drops below this limit. Can have up to num_xa_central_limits of these (see star_def.inc for value). xa_central_lower_limit_species contains an isotope name as defined in chem_def.f90. xa_central_lower_limit contains the lower limit value.

xa_central_lower_limit_species(1) = ''
xa_central_lower_limit(1) = 0

xa_central_upper_limit_species#

xa_central_upper_limit#

Upper limits on central mass fractions. Stop when central abundance rises above this limit. Can have up to num_xa_central_limits of these (see star_def.inc for value). E.g., to stop when center c12 abundance reaches 0.5, set

xa_central_upper_limit_species(1) = 'c12'
xa_central_upper_limit(1) = 0.5
xa_central_upper_limit_species(1) = ''
xa_central_upper_limit(1) = 0

xa_surface_lower_limit_species#

xa_surface_lower_limit#

Lower limits on surface mass fractions. Stop when surface abundance drops below this limit. Can have up to num_xa_surface_limits of these (see star_def for value) xa_surface_lower_limit_species contains an isotope name as defined in chem_def.f90 xa_surface_lower_limit contains the lower limit value

xa_surface_lower_limit_species(1) = ''
xa_surface_lower_limit(1) = 0

xa_surface_upper_limit_species#

xa_surface_upper_limit#

upper limits on surface mass fractions stop when surface abundance rises above this limit can have up to num_xa_surface_limits of these (see star_def for value) e.g., to stop when surface c12 abundance reaches 0.5, set

xa_surface_upper_limit_species(1) = 'c12'
xa_surface_upper_limit(1) = 0.5
xa_surface_upper_limit_species(1) = ''
xa_surface_upper_limit(1) = 0

xa_average_lower_limit_species#

xa_average_lower_limit#

lower limits on average mass fractions stop when average abundance drops below this limit can have up to num_xa_average_limits of these (see star_def for value)

xa_average_lower_limit_species(1) = ''
xa_average_lower_limit(1) = 0

xa_average_upper_limit_species#

xa_average_upper_limit#

upper limits on average mass fractions stop when average abundance rises above this limit can have up to num_xa_average_limits of these (see star_def for value)

xa_average_upper_limit_species(1) = ''
xa_average_upper_limit(1) = 0

HB_limit#

For detecting when the model has completed the horizontal branch. Only applies when center abundance by mass of h1 is < 1d-4. Stop when the center abundance by mass of he4 drops below this limit.

HB_limit = 0

star_mass_min_limit#

Stop when star mass in Msun units is < this. <= 0 means no limit.

star_mass_min_limit = 0

star_mass_max_limit#

Stop when star mass in Msun units is > this. <= 0 means no limit.

star_mass_max_limit = 0

remnant_mass_min_limit#

ejecta_mass_max_limit#

Stop when remnant mass in Msun units is < this. <= 0 means no limit. remnant_mass = star_mass - ejecta_mass ejecta_mass is outermost mass that all has v(k) >= v_escape(k)

remnant_mass_min_limit = 0
ejecta_mass_max_limit = 1d99

star_species_mass_min_limit#

star_species_mass_min_limit_iso#

star_species_mass_max_limit#

star_species_mass_max_limit_iso#

Stop when a particular species mass in units of Msun exceeds these limits.

star_species_mass_min_limit = 0
star_species_mass_min_limit_iso = ''
star_species_mass_max_limit = 0
star_species_mass_max_limit_iso = ''

star_H_mass_max_limit#

Stop when star hydrogen mass in Msun units is > this. <= 0 means no limit.

star_H_mass_min_limit replaced by star_species_mass_min_limit star_He_mass_min_limit replaced by star_species_mass_min_limit star_C_mass_min_limit replaced by star_species_mass_min_limit star_H_mass_max_limit replaced by star_species_mass_max_limit star_He_mass_max_limit replaced by star_species_mass_max_limit star_C_mass_max_limit replaced by star_species_mass_max_limit

envelope_mass_limit#

envelope_mass = star_mass - he_core_mass

Stop when envelope_mass drops below this limit, in Msun units.

envelope_mass_limit = 0

envelope_fraction_left_limit#

envelope_fraction_left = (star_mass - he_core_mass)/(initial_mass - he_core_mass)

Stop when envelope_fraction_left < this limit.

envelope_fraction_left_limit = 0

xmstar_min_limit#

xmstar = mstar - M_center

stop when xmstar in grams is < this. <= 0 means no limit.

xmstar_min_limit = 0

xmstar_max_limit#

xmstar = mstar - M_center

stop when xmstar in grams is > this. <= 0 means no limit.

xmstar_max_limit = 0

he_core_mass_limit#

stop when helium core reaches this mass, in Msun units

he_core_mass_limit = 1d99

co_core_mass_limit#

stop when CO core reaches this mass, in Msun units

co_core_mass_limit = 1d99

one_core_mass_limit#

stop when ONe core reaches this mass, in Msun units

one_core_mass_limit = 1d99

fe_core_mass_limit#

stop when iron core reaches this mass, in Msun units

fe_core_mass_limit = 1d99

neutron_rich_core_mass_limit#

stop when neutron rich core reaches this mass, in Msun units

neutron_rich_core_mass_limit = 1d99

he_layer_mass_lower_limit#

he layer mass is defined as he_core_mass - c_core_mass stop when c_core_mass > 0 and he layer mass < this limit (Msun units).

he_layer_mass_lower_limit = 0

abs_diff_lg_LH_lg_Ls_limit#

stop when abs(lg_LH - lg_Ls) <= abs_diff_LH_Lsurf_limit can be useful for deciding when pre-main sequence star has reached ZAMS set to negative value to disable

abs_diff_lg_LH_lg_Ls_limit = -1

Teff_upper_limit#

stop when Teff is greater than this limit.

Teff_upper_limit = 1d99

Teff_lower_limit#

stop when Teff is less than this limit.

Teff_lower_limit = -1d99

photosphere_r_upper_limit#

stop when photosphere_r is greater than this limit, in Rsun units

photosphere_r_upper_limit = 1d99

photosphere_r_lower_limit#

stop when photosphere_r is less than this limit, in Rsun units

photosphere_r_lower_limit = -1d99

photosphere_m_upper_limit#

stop when photosphere_m is greater than this limit, in Msun units

photosphere_m_upper_limit = 1d99

photosphere_m_lower_limit#

stop when photosphere_m is less than this limit, in Msun units

photosphere_m_lower_limit = -1d99

photosphere_m_sub_M_center_limit#

stop when photosphere_m is less than this limit above M_center, in Msun units

photosphere_m_sub_M_center_limit = -1d99

log_Teff_upper_limit#

stop when log10 of Teff is greater than this limit.

log_Teff_upper_limit = 1d99

log_Teff_lower_limit#

stop when log10 of Teff is less than this limit.

log_Teff_lower_limit = -1d99

log_Tsurf_upper_limit#

stop when log10 of T in outermost cell is greater than this limit.

log_Tsurf_upper_limit = 1d99

log_Tsurf_lower_limit#

stop when log10 of T in outermost cell is less than this limit.

log_Tsurf_lower_limit = -1d99

log_Rsurf_upper_limit#

stop when log10 of R/Rsun in outermost cell is greater than this limit.

log_Rsurf_upper_limit = 1d99

log_Rsurf_lower_limit#

stop when log10 of R/Rsun in outermost cell is less than this limit.

log_Rsurf_lower_limit = -1d99

log_L_upper_limit#

stop when log10(total luminosity in Lsun units) is greater than this limit. in order to skip pre-ms, this limit only applies when L_nuc > 0.01*L

log_L_upper_limit = 1d99

log_L_lower_limit#

stop when log10(total luminosity in Lsun units) is less than this limit.

log_L_lower_limit = -1d99

log_g_upper_limit#

stop when log10(gravity at surface) is greater than this limit.

log_g_upper_limit = 1d99

log_g_lower_limit#

stop when log10(gravity at surface) is less than this limit.

log_g_lower_limit = -1d99

log_Psurf_upper_limit#

stop when log10 of surface pressure is greater than this limit.

log_Psurf_upper_limit = 1d99

log_Psurf_lower_limit#

stop when log10 of surface pressure is less than this limit.

log_Psurf_lower_limit = -1d99

log_Dsurf_upper_limit#

stop when log10 of surface density is greater than this limit.

log_Dsurf_upper_limit = 1d99

log_Dsurf_lower_limit#

stop when log10 of surface density is less than this limit.

log_Dsurf_lower_limit = -1d99

power_nuc_burn_upper_limit#

stop when total power from all nuclear reactions (in Lsun units) is > this.

power_nuc_burn_upper_limit = 1d99

power_h_burn_upper_limit#

stop when total power from hydrogen-consuming reactions (in Lsun units) is > this.

power_h_burn_upper_limit = 1d99

power_he_burn_upper_limit#

stop when total power from reactions burning helium (in Lsun units) is > this.

power_he_burn_upper_limit = 1d99

power_z_burn_upper_limit#

stop when total power from reactions burning metals (in Lsun units) is > this

power_z_burn_upper_limit = 1d99

power_nuc_burn_lower_limit#

stop when total power from all nuclear reactions (in Lsun units) is < this.

power_nuc_burn_lower_limit = -1d99

power_h_burn_lower_limit#

stop when total power from hydrogen consuming reactions (in Lsun units) is < this.

power_h_burn_lower_limit = -1d99

power_he_burn_lower_limit#

stop when total power from reactions burning helium (in Lsun units) is < this.

power_he_burn_lower_limit = -1d99

power_z_burn_lower_limit#

stop when total power from reactions burning metals (in Lsun units) is < this.

power_z_burn_lower_limit = -1d99

max_abs_rel_run_E_err#

Stop if the abs value of cumulative_energy_error/total_energy exceeds this value. Ignore if < 0. Also ignore during relax operations.

max_abs_rel_run_E_err = -1

max_number_retries#

Stop if the number of retries exceeds this value. Ignore if < 0.

max_number_retries = -1
relax_max_number_retries = 300

min_timestep_limit#

stop if need timestep smaller than this limit, in seconds

min_timestep_limit = 1d-6

center_Ye_lower_limit#

stop if center_ye drops below this limit

center_Ye_lower_limit = -1

center_R_lower_limit#

stop if R_center drops below this limit (in cm)

center_R_lower_limit = -1

fe_core_infall_limit#

stop if at least fe_core_infall_mass of material has speed greater than this, at a location interior to fe_core_mass, in cm/s

fe_core_infall_limit = 3d7

fe_core_infall_mass#

Amount of mass to check if collapsing, the smaller this is the closer the velocity minima will be to fe_core_infall but there will be a greater chance of a transistent velocity spike causing the model to prematurely exit. In solar masses

fe_core_infall_mass = 0.1d0

non_fe_core_infall_limit#

stop if at least non_fe_core_infall_mass of material has speed greater than this, at a location interior to he_core_mass and exterior to fe_core_mass. in cm/s

non_fe_core_infall_limit = 1d99

non_fe_core_infall_mass#

Amount of mass to check if collapsing, the smaller this is the closer the velocity minima will be to non_fe_core_infall but there will be a greater chance of a transistent velocity spike causing the model to prematurely exit. In solar masses

non_fe_core_infall_mass = 0.1d0

non_fe_core_rebound_limit#

stop if max rebound speed (outward) at any location interior to he_core_mass. in cm/s

non_fe_core_rebound_limit = 1d99

v_div_csound_max_limit#

stop if any v/csound > this limit

v_div_csound_max_limit = 1d99

v_div_csound_surf_limit#

stop if v_surf/csound_surf > this limit

v_div_csound_surf_limit = 1d99

v_surf_div_v_kh_upper_limit#

stop if abs(v_surf/v_kh) > this limit, where v_kh = photosphere_r/kh_timescale

v_surf_div_v_kh_upper_limit = 1d99

v_surf_div_v_kh_lower_limit#

stop if abs(v_surf/v_kh) < this limit, where v_kh = photosphere_r/kh_timescale

v_surf_div_v_kh_lower_limit = -1d99

v_surf_div_v_esc_limit#

stop if v_surf/v_esc > this limit

v_surf_div_v_esc_limit = 1d99

v_surf_kms_limit#

stop if v_surf in km/s > this limit

v_surf_kms_limit = 1d99

Lnuc_div_L_zams_limit#

defines “near zams” – note: must also set stop_near_zams

Lnuc_div_L_zams_limit = 0.9d0

stop_near_zams#

if true, stop if Lnuc/L > Lnuc_div_L_zams_limit

stop_near_zams = .false.

stop_at_phase_PreMS#

stop_at_phase_ZAMS#

stop_at_phase_IAMS#

stop_at_phase_TAMS#

stop_at_phase_He_Burn#

stop_at_phase_ZACHeB#

stop_at_phase_TACHeB#

stop_at_phase_TP_AGB#

stop_at_phase_C_Burn#

stop_at_phase_Ne_Burn#

stop_at_phase_O_Burn#

stop_at_phase_Si_Burn#

stop_at_phase_WDCS#

if true, terminate model when phase of evolution reaches this point. Definitions for stop_at_phase_* can be found in $MESA_DIR/star/private/star_utils.f90 inside the subroutine set_phase_of_evolution

::

stop_at_phase_PreMS = .false. stop_at_phase_ZAMS = .false. stop_at_phase_IAMS = .false. stop_at_phase_TAMS = .false. stop_at_phase_He_Burn = .false. stop_at_phase_ZACHeB = .false. stop_at_phase_TACHeB = .false. stop_at_phase_TP_AGB = .false. stop_at_phase_C_Burn = .false. stop_at_phase_Ne_Burn = .false. stop_at_phase_O_Burn = .false. stop_at_phase_Si_Burn = .false. stop_at_phase_WDCS = .false.

Lnuc_div_L_upper_limit#

stop when Lnuc/L is greater than this limit. Here, Lnuc refers to the total thermal power from all burning, including photodisintegrations, power_nuc_burn in the history_columns.list file. L is the luminosity at the photosphere.

Lnuc_div_L_upper_limit = 1d99

Lnuc_div_L_lower_limit#

stop when Lnuc/L is less than this limit. Here, Lnuc refers to the total thermal power from all burning, including photodisintegrations, power_nuc_burn in the history_columns.list file. L is the luminosity at the photosphere.

Lnuc_div_L_lower_limit = -1d99

gamma1_limit#

stop if min gamma1 < this limit in a cell with q <= gamma1_limit_max_q

gamma1_limit = -1

gamma1_limit_max_q#

stop if gamma1 < this limit at any location with q <= gamma1_limit_max_q values near unity skip the outer envelope

gamma1_limit_max_q = 0.95d0

gamma1_limit_max_v_div_vesc#

stop if gamma1 < this limit at any location with v_div_vesc <= gamma1_limit_max_v_div_vesc

gamma1_limit_max_v_div_vesc = 0.95d0

Pgas_div_P_limit#

stop when Pgas/P <= Pgas_div_P_limit

Pgas_div_P_limit = 0

Pgas_div_P_limit_max_q#

stop if Pgas/P < Pgas_div_P_limit at any location with q <= Pgas_div_P_limit_max_q values near unity skip the outer envelope

Pgas_div_P_limit_max_q = 0.95d0

peak_burn_vconv_div_cs_limit#

limits ratio of convection velocity to sound speed at location of peak eps_nuc

peak_burn_vconv_div_cs_limit = 1d99

omega_div_omega_crit_limit#

stop if omega/omega_crit is > this anywhere in star ignore if < 0

omega_div_omega_crit_limit = -1

delta_nu_lower_limit#

stop when asteroseismology delta_nu in micro Hz is < this. <= 0 means no limit.

delta_nu_lower_limit = 0

delta_nu_upper_limit#

stop when asteroseismology delta_nu in micro Hz is > this. <= 0 means no limit.

delta_nu_upper_limit = 0

shock_mass_upper_limit#

stop when shock_mass is > this. <= 0 means no limit.

shock_mass_upper_limit = -1

mach1_mass_upper_limit#

stop when outer location of mach 1 is > this. <= 0 means no limit.

mach1_mass_upper_limit = -1

delta_Pg_lower_limit#

stop when delta_Pg in micro Hz is < this. <= 0 means no limit.

delta_Pg_lower_limit = 0

delta_Pg_upper_limit#

stop when delta_Pg in micro Hz is > this. <= 0 means no limit.

delta_Pg_upper_limit = 0

stop_when_reach_this_cumulative_extra_heating#

(ignore if <= 0)

stop_when_reach_this_cumulative_extra_heating = 0d0

mixing parameters#

mixing_length_alpha#

The mixing length is this parameter times a local pressure scale height. To increase R vs. L, decrease mixing_length_alpha.

mixing_length_alpha = 2

remove_small_D_limit#

If MLT diffusion coeff D (cm^2/sec) is less than this limit, then set D to zero and change the point to mixing_type == no_mixing.

remove_small_D_limit = 1d-6

use_Ledoux_criterion#

a location in the model is Schwarzschild stable when gradr < grada it is Ledoux stable when gradr < gradL, where gradL = grada + composition_gradient note that these are the same when composition_gradient = 0 so you can force the use of the Schwarzschild criterion by passing 0 for the composition_gradient argument to the mlt routine. that’s what happens if you set the control “use_Ledoux_criterion” to false.

overshooting and rotational mixing are dealt with separately

and are added after the MLT classifications are made.

use_Ledoux_criterion = .false.

num_cells_for_smooth_gradL_composition_term#

Number of cells on either side to use in weighted smoothing of gradL_composition_term. gradL_composition_term is set to the “raw” unsmoothed brunt_B and then optionally smoothed according num_cells_for_smooth_gradL_composition_term. In cases where the Ledoux criterion is used to evaluate the boundary for burning convective cores, you may need to set num_cells_for_smooth_gradL_composition_term = 0 to avoid smoothing the stabilizing composition jump into the convection zone and unphysically causing it to shrink. See section 3.2 in Moore, K., & Garaud, P. 2016, APJ, 817, 54

num_cells_for_smooth_gradL_composition_term = 3

threshold_for_smooth_gradL_composition_term#

Threshold for weighted smoothing of gradL_composition_term. Only apply smoothing (controlled by num_cells_for_smooth_gradL_composition_term) for contiguous regions where \(|\nabla_L|\) exceeds this threshold. Might be useful for preventing narrow composition jumps from being excessively broadened by smoothing

threshold_for_smooth_gradL_composition_term = 0

alpha_semiconvection#

Determines efficiency of semiconvective mixing. Semiconvection only applies if use_Ledoux_criterion is true.

alpha_semiconvection = 0

semiconvection_option#

  • 'Langer_85 mixing; gradT = gradr' : uses Langer scheme for mixing but sets gradT = gradr

  • 'Langer_85' : this calculates special gradT as well as doing mixing.

semiconvection_option = 'Langer_85 mixing; gradT = gradr'

thermohaline_coeff#

Determines efficiency of thermohaline mixing. was previously named thermo_haline_coeff. thermohaline mixing only applies if use_Ledoux_criterion is true.

thermohaline_coeff = 0

thermohaline_option#

determines which method to use for calculating thermohaline diffusion coef:

  • 'Kippenhahn' : use method of Kippenhahn, R., Ruschenplatt, G., & Thomas, H.-C. 1980, A&A, 91, 175.

  • 'Traxler_Garaud_Stellmach_11' : use method of Traxler, Garaud, & Stellmach, ApJ Letters, 728:L29 (2011).

  • 'Brown_Garaud_Stellmach_13' : use method of Brown, Garaud, & Stellmach, ApJ 768:34 (2013) Recommends thermohaline_coeff = 1, but it can nevertheless be changed.

thermohaline_option = 'Kippenhahn'

alt_scale_height_flag#

If false, then stick to the usual definition – P/(g*rho). If true, use min of the usual and sound speed * hydro time scale, sqrt(P/G)/rho. Note that the ‘TDC’ MLT_option does not respect the alt_scale_height option, and continues to use h = P / rho g even if that flag is set.

alt_scale_height_flag = .true.

mlt_use_rotation_correction#

When doing rotation, multiply grad_rad by ft_rot/ft_rot if this flag is true.

mlt_use_rotation_correction = .true.

mlt_Pturb_factor#

include MLT turbulent pressure at face k = mlt_Pturb_factor*0.5*(rho(k) + rho(k-1))*mlt_vc(k)**2/3 MLT turbulent pressure for cell k = avg of values at faces.

this replaces conv_dP_term_factor. also see extra_pressure vector and other_pressure routine

mlt_Pturb_factor = 0d0

MLT_option#

Options are:

  • ‘none’ : just give radiative values with no mixing.

  • ‘TDC’ : Time-dependent convection based on the Khufuss 1986 model. Reduces to Cox at long times.

  • ‘Cox’ : MLT as developed in Cox & Giuli 1968, Chapter 14.

  • ‘ML1’ : Bohm-Vitense 1958

  • ‘ML2’ : Bohm and Cassinelli 1971

  • ‘Mihalas’ : Mihalas 1978, Kurucz 1979

  • ‘Henyey’ : Henyey, Vardya, and Bodenheimer 1965

The ‘Cox’ and ‘TDC’ options assumes optically thick material. The other options are various ways of extending to include optically thin material. Note that TDC does not respect the alt_scale_height option, and continues to use h = P / rho g even if that flag is set.

We caution that combining different mixing models in a stellar evolution calculation might lead to physically inconsistent solutions, because the different models have been developed separately, and their underlying assumptions might not be compatible with each other. Examples include combining the newly implemented time-dependent local limit convection model TDC with an overshooting model, or combining TDC with other models for chemical composition gradients (predictive mixing or convective premixing), rotation, etc. (MESA VI)

MLT_option = 'TDC'

TDC options#

  • alpha_TDC_DAMP : The turbulent viscous damping parameter which determines the saturation of TDC. Increasing this decreases convection speeds.

  • alpha_TDC_DAMPR : The radiative damping parameter which determines the saturation of TDC. Increasing this decreases convection speeds.

  • alpha_TDC_PtdVdt : The prefactor on the term accounting for work done against turbulent pressure (P_turb * dV/dt). Physically this should be unity.

  • steps_before_use_TDC : TDC often struggles with models on the pre-main-sequence. Set this option to pick MLT_option=’Cox’ for the first several steps to get past the pre-MS. Note that if this option is positive then only either TDC or Cox will be used (depending on model number). THIS OVERRIDES MLT_option!

alpha_TDC_DAMP = 1d0
alpha_TDC_DAMPR = 0d0
alpha_TDC_PtdVdt = 0d0
steps_before_use_TDC = 0

Henyey_MLT_y_param#

Henyey_MLT_nu_param#

Values of the f1..f4 coefficients are taken from Table 1 of Ludwig et al. 1999, A&A, 346, 111 with the following exception: their value of f3 for Henyey convection is f4/8 when it should be 8*f4, i.e., f3=32*pi**2/3 and f4=4*pi**2/3. f3 and f4 are related to the henyey y parameter, so for the ‘Henyey’ case they are set based on the value of Henyey_y_param.

Henyey_MLT_y_param = 0.33333333d0
Henyey_MLT_nu_param = 8

make_gradr_sticky_in_solver_iters#

min_logT_for_make_gradr_sticky_in_solver_iters#

if true, then location that becomes radiative during solver iterations, stays radiative for rest of the solver iterations. to avoid flip-flopping between radiative and convective. also do this if max logT >= min_logT_for_make_gradr_sticky_in_solver_iters

make_gradr_sticky_in_solver_iters = .false.
min_logT_for_make_gradr_sticky_in_solver_iters = 1d99

no_MLT_below_shock#

if true, then no MLT below an outward going shock (just radiative).

no_MLT_below_shock = .false.

mlt_make_surface_no_mixing#

mlt_make_surface_no_mixing = .false.

T_mix_limit#

If there is any convection in surface zones with T < T_mix_limit, then extend the innermost such convective region outward all the way to the surface. For example,

  • T_mix_limit <= 0 means omit this operation.

  • T_mix_limit = 1d5 will effectively make the star convective down to the He++ region.

units in Kelvin

T_mix_limit = 0

mlt_gradT_fraction#

let f := mlt_gradT_fraction if f is >= 0 and <= 1, then gradT from mlt is replaced by f*grada_at_face(k) + (1-f)*gradr(k) see also the vector control adjust_mlt_gradT_fraction for fine grain control

mlt_gradT_fraction = -1

okay_to_reduce_gradT_excess#

gradT_excess = gradT_sub_grada = superadiabaticity.

Inefficient convection => large gradT excess and steep T gradient to enhance radiative transport. Reduce gradT excess by making gradT closer to adiabatic gradient. If true, code is allowed to adjust gradT to boost efficiency of energy transport See gradT_excess_f1, gradT_excess_f2, and gradT_excess_age_fraction below.

This is the treatment of convection, referred to as MLT++ in Section 7.2 of Paxton et al. (2013), that reduces the superadiabaticity in some radiation-dominated convective regions.

okay_to_reduce_gradT_excess = .false.

gradT_excess_f1#

gradT_excess_f2#

These are for calculation of efficiency boosted gradT.

gradT_excess_f1 = 1d-4
gradT_excess_f2 = 1d-3

gradT_excess_age_fraction#

These are for calculation of efficiency boosted gradT. Fraction of old to mix with new to get next.

gradT_excess_age_fraction = 0.9d0

gradT_excess_max_change#

These are for calculation of efficiency boosted gradT. Maximum change allowed in one timestep for gradT_excess_alpha. Ignored if negative.

gradT_excess_max_change = -1d0

gradT_excess_lambda1#

gradT_excess_beta1#

In some situations you might want to force alfa = 1. You can do that by setting gradT_excess_lambda1 < 0. The following are for the normal calculation of gradT_excess_alfa

gradT_excess_lambda1 = 1.0d0
gradT_excess_beta1 = 0.35d0

gradT_excess_lambda2#

gradT_excess_beta2#

The following are for the normal calculation of gradT_excess_alfa.

gradT_excess_lambda2 = 0.5d0
gradT_excess_beta2 = 0.25d0

gradT_excess_dlambda#

gradT_excess_dbeta#

The following are for the normal calculation of gradT_excess_alfa.

gradT_excess_dlambda = 0.1d0
gradT_excess_dbeta = 0.1d0

gradT_excess_max_center_h1#

No boost if center H1 > this limit.

gradT_excess_max_center_h1 = 1d0

gradT_excess_min_center_he4#

No boost if center He4 < this limit.

gradT_excess_min_center_he4 = 0d0

gradT_excess_max_logT#

No local boost if local logT > this limit.

gradT_excess_max_logT = 8

gradT_excess_min_log_tau_full_on#

gradT_excess_max_log_tau_full_off#

No local boost if local log_tau < gradT_excess_max_log_tau_full_off. Reduced local boost if local log_tau < gradT_excess_min_log_tau_full_on.

gradT_excess_min_log_tau_full_on = -99
gradT_excess_max_log_tau_full_off = -99

### max_logT_for_mlt

No mlt at cell if local logT > this limit.

max_logT_for_mlt = 99

use_superad_reduction#

superad_reduction_Gamma_limit#

superad_reduction_Gamma_limit_scale#

superad_reduction_Gamma_inv_scale#

superad_reduction_diff_grads_limit#

superad_reduction_limit#

Implicit alternative to okay_to_reduce_gradT_excess

use_superad_reduction = .false.
superad_reduction_Gamma_limit = 0.5d0
superad_reduction_Gamma_limit_scale = 5d0
superad_reduction_Gamma_inv_scale = 5d0
superad_reduction_diff_grads_limit = 1d-3
superad_reduction_limit = -1d0

overshooting#

There are two schemes implemented in MESA to treat overshooting: a step overshoot scheme and an exponential scheme.

Parameters for exponential diffusive overshoot are described in the paper by Falk Herwig, “The evolution of AGB stars with convective overshoot”, A&A, 360, 952-968 (2000).

Overshooting depends on the classification of the convective zone and can be different at the top and the bottom of the zone.

The overshooting controls are based on convection-zone and convection-boundary matching criteria. These criteria are overshoot_zone_type, overshoot_zone_loc and overshoot_bdy_loc. The overshooting parameter values corresponding to the first set of matching criteria will be used. Therefore, narrower criteria should precede more general ones (i.e have lower array indices).

These are arrays of size NUM_OVERSHOOT_PARAM_SETS which is defined in star_data/public/star_data_def.inc (currently 16)

overshoot_scheme(:) = '' ! ``exponential``, ``step``, ``other``
overshoot_zone_type(:) = '' !  ``burn_H``, ``burn_He``, ``burn_Z``, ``nonburn``, ``any``
overshoot_zone_loc(:) = '' ! ``core``, ``shell``, ``any``
overshoot_bdy_loc(:) = '' ! ``bottom``, ``top``, ``any``

Amount of overshooting from edge of convective zone These are arrays of size NUM_OVERSHOOT_PARAM_SETS which is defined in star_data/public/star_data_def.inc (currently 16)

overshoot_f(:) = 0d0
overshoot_f0(:) = 0d0

The switch from convective mixing to overshooting happens at a distance f0*Hp into the convection zone from the estimated location where grad_ad = grad_rad, where Hp is the pressure scale height at that location. A value <= 0 for f0 is a mistake – you are required to set f0 as well as f. take a look at the following from an email concerning this: Overshooting works by taking the diffusion mixing coefficient at the edge of the convection zone and extending it beyond the zone. But – and here’s the issue – at the exact edge of the zone the mixing coefficient goes to 0. So we don’t want that. Instead we want the value of the mixing coeff NEAR the edge, but not AT the edge. The “f0” parameter determines the exact meaning of “near” for this. It tells the code how far back into the zone to go in terms of scale height. The overshooting actually begins at the location determined by f0 back into the convection zone rather than at the edge where the diffusion coeff is ill-defined. So, for example, if you want overshooting of 0.2 scale heights beyond the normal edge, you might want to back up 0.05 scale heights to get the diffusion coeff from near the edge and then go out by 0.25 scale heights from there to reach 0.2 Hp beyond the old boundary. In the inlist this would mean setting the “f0” to 0.05 and the “f” to 0.25.

For step overshoot: overshooting extends a distance overshoot_f*Hp0 from r0 with constant diffusion coefficient D = overshoot_D0 + overshoot_Delta0*D_ob where D_ob is the convective diffusivity at the bottom (top) of the step overshoot region for outward (inward) overshooting. These are arrays of size NUM_OVERSHOOT_PARAM_SETS which is defined in star_data/public/star_data_def.inc (currently 16)

overshoot_D0(:) = 0d0
overshoot_Delta0(:) = 1d0

You can specify a range of star masses over which overshooting above H burning zones is gradually enabled. Do specified overshooting above H burning zone if star_mass >= this (Msun). These are arrays of size NUM_OVERSHOOT_PARAM_SETS which is defined in star_data/public/star_data_def.inc (currently 16)

overshoot_mass_full_on(:) = 0d0

You can specify a range of star masses over which overshooting above H burning zones is gradually enabled. No overshooting above H burning zone if star_mass <= this (Msun). These are arrays of size NUM_OVERSHOOT_PARAM_SETS which is defined in star_data/public/star_data_def.inc (currently 16)

overshoot_mass_full_off(:) = 0d0

overshoot_D_min#

Overshooting shuts off when the exponential decay has dropped the diffusion coefficient to this level.

overshoot_D_min = 1d2

overshoot_brunt_B_max#

Terminate overshoot region when encounter stabilizing composition gradient where (unsmoothed) brunt_B is greater than this limit. (<= 0 means ignore this limit) note: both brunt_B and gradL_composition_term come from unsmoothed_brunt_B and differ only in optional smoothing. (see num_cells_for_smooth_brunt_B and num_cells_for_smooth_gradL_composition_term).

overshoot_brunt_B_max = 0d0

min_overshoot_q#

Overshooting is only allowed at locations with mass m >= min_overshoot_q * mstar. E.g., if min_overshoot_q = 0.1, then only the outer 90% by mass can have overshooting. This provides a simple way of suppressing bogus center overshooting in which a small convective region at the core can produce excessively large overshooting because of a large pressure scale height at the center.

min_overshoot_q = 0d0

NOTE: In addition to giving these ‘f’ parameters non-zero values, you should also check the settings for mass_for_overshoot_full_on and mass_for_overshoot_full_off.

overshoot_alpha#

The value of Hp for overshooting is limited to the radial thickness of the convection zone divided by overshoot_alpha. only used when > 0. if <= 0, then use mixing_length_alpha instead.

overshoot_alpha = -1

limit_overshoot_Hp_using_size_of_convection_zone#

if false, allow large distance of overshoot for small convective zones.

limit_overshoot_Hp_using_size_of_convection_zone = .true.

burn_z_mix_region_logT#

burn_he_mix_region_logT#

burn_h_mix_region_logT#

max logT in convective region determines burn type for overshooting

burn_z_mix_region_logT = 8.7d0
burn_he_mix_region_logT = 7.7d0
burn_h_mix_region_logT = 6.7d0

max_Y_for_burn_z_mix_region#

max_X_for_burn_he_mix_region#

Even if a region reaches the above temperature to be considered as a z_burn region, only set it as such if the helium mass fraction in all points of the region is below max_Y_for_burn_z_mix_region. Similarly, max_X_for_burn_he_mix_region controls if a region is considered as a he_burn region in terms of the hydrogen mass fraction.

max_Y_for_burn_z_mix_region = 1d-4
max_X_for_burn_he_mix_region = 1d-4

Predictive mixing#

Predictive mixing is an approach for expanding convective boundaries until gradr = grada on the convective side of the boundary (as required by the criterion that the convective velocity and luminosity vanish at the boundary). It is discussed in detail in Paxton et al. 2018, ApJ, in press: “Modules for Experiments in Stellar Astrophysics (MESA): Convective boundaries, element diffusion, and massive star explosions”

Predictive mixing is controlled by specifying a set of parameters, which combines matching criteria (determining which boundaries to apply the predictive mixing to) together with values (determining how the predictive mixing should operate at those boundaries). Up to NUM_PREDICTIVE_PARAM_SETS of these parameter sets can be defined (see star_def.inc for value).

predictive_mix#

Set to .true. to enable this set of parameters

predictive_mix(1) = .false.

predictive_zone_type#

Matching criterion for the type of the convection zone. Possible values are burn_H (hydrogen burning), burn_He (helium burning), burn_Z (metal burning), nonburn (no burning) or any (which matches any type of zone).

predictive_zone_type(1) = ''

predictive_zone_loc#

Matching criterion for the location of the convection zone. Possible values are core (the core convection zone), shell (a convective shell), surf (the surface convection zone) or any (which matches any location).

predictive_zone_loc(1) = ''

predictive_bdy_loc#

Matching criterion for the location of the convective boundary. Possible values are top (the top of the convection zone), bottom (the bottom of the convection zone) or any (which matches any location).

predictive_bdy_loc(1) = ''

predictive_bdy_q_min#

Matching criterion for the minimum fractional mass coordinate of the convective boundary

predictive_bdy_q_min(1) = 0d0

predictive_bdy_q_max#

Matching criterion for the maximum fractional mass coordinate of the convective boundary

predictive_bdy_q_max(1) = 1d0

predictive_superad_thresh#

Threshold for minimum superadiabaticity in the predictive mixing scheme; boundary expansion stops when gradr/grada-1 drops below this threshold. Default value is usually good for main-sequence evolution; for core He-burning, set to 0.005, 0.01 or larger to prevent splitting of the core convection zone and/or core breathing pulses.

predictive_superad_thresh(1) = 0d0

predictive_avoid_reversal#

Species to monitor for reversals in abundance evolution. If this is set to the name of a species, then the predictive mixing scheme will try to avoid causing reversals in the abundance of that species (e.g., changing the abundance evolution from decreasing to increasing). Set to ‘he4’ during core He-burning to prevent splitting of the core convection zone and/or core breathing pulses.

predictive_avoid_reversal(1) = ''

predictive_limit_ingestion#

predictive_ingestion_factor#

Limit the rate of ingestion of a species, following the prescription given in equation (2) of Constantino, Campbell & Lattanzio (2017, MNRAS, 472, 4900). The control predictive_limit_ingestion specifies which species to limit, and the control predictive_ingestion_factor gives the multiplying factor. Setting this factor to 5/12 is the same as choosing alpha_i = 1 in their equation (2).

predictive_limit_ingestion(1) = ''
predictive_ingestion_factor(1) = 0d0

Convective premixing#

Convective premixing is a approach to handling mixing in convection zones that improves upon the predictive mixing scheme described above. Like predictive mixing, it expands convective boundaries until gradr = grada on the convective side of the boundary. Unlike predictive mixing, it directly modifies abundances in the stellar model, via a iterative series of mixings-to-homogeneity over a shifting window of cells. This iterative approach allows convective premixing to build ‘classical’ semiconvection regions, where the abundance gradient is tuned to yield convective neutrality.

NOTE: the history columns that give info on the convective and semi convective boundaries (i.e., mass_conv_core and mass_semiconv_core) do not work with CPM. Instead, one should look at the profiles to see where the boundaries are.

do_conv_premix#

Set to .true. to perform convective premixing. Note that this cannot be enabled at the same time as the predictive_mix control

do_conv_premix = .false.

do_premix_heating#

if true, calculate heating term associated with changes in internal energy due to any abundance changes from convective premixing, and include this term in the energy equation.

do_premix_heating = .true.

conv_premix_avoid_increase#

Attempt to avoid increases in the abundance of species being burned. Sometimes, the convective premixing scheme can cause the abundance of a species being burned (e.g., helium during core helium burning) to increase across a timestep. This typically arises when the scheme mixes a fresh supply of the species into the convection zone where the burning is occurring. Setting the conv_premix_avoid_increase control to .true. will tell the scheme to avoid such outcomes, if possible. In the case of core helium burning, this helps to reduces the incidences of core breathing pulses (although in some situations it doesn’t completely eliminate them).

conv_premix_avoid_increase = .false.

conv_premix_time_factor#

Scaling factor for deciding whether a convective boundary has enough time to advance during a timestep. Simple physical arguments suggest that a convective boundary requires a time delta_t ~ delta_r/v_conv to advance a distance delta_r. The convective premixing algorithm keeps a tally of how much time a boundary has spent advancing, and it prevents further advancing if this time would exceed conv_premix_time_factor*dt, where dt is the current timestep. Setting conv_premix_time_factor to a value <= 0 disables this check. STILL UNDER DEVELOPMENT, AND DISABLED BY DEFAULT

conv_premix_time_factor = 0.0

conv_premix_fix_pgas#

Flag to decide whether gas pressure is kept constant during premixing (.true.), or instead density is kept constant (.false.). In both cases, temperature is kept constant.

conv_premix_fix_pgas = .true.

conv_premix_dump_snapshots#

Flag to write out snapshots of the intermediate stages during the convective premixing iterations. Refer to the dump_snapshot_ routine in star/private/conv_premix.f90 to see what’s written out. Note that this can quickly fill up your disk!

conv_premix_dump_snapshots = .false.

Rayleigh-Taylor Instability#

derived from Paul Duffell’s code RT1D.

RTI_smooth_mass#

RTI_smooth_iterations#

RTI_smooth_fraction#

smoothing for dPdr_dRhodr_info done at start of step

RTI_smooth_mass = 0d0
RTI_smooth_iterations = 0
RTI_smooth_fraction = 1d0

alpha_RTI_diffusion_factor#

dudt_RTI_diffusion_factor#

dedt_RTI_diffusion_factor#

dlnddt_RTI_diffusion_factor#

composition_RTI_diffusion_factor#

max_M_RTI_factors_full_on#

min_M_RTI_factors_full_off#

alpha_RTI_diffusion_factor = 1d0
dudt_RTI_diffusion_factor = 1d0
dedt_RTI_diffusion_factor = 1d0
dlnddt_RTI_diffusion_factor = 1d0
composition_RTI_diffusion_factor = 1d0
max_M_RTI_factors_full_on = 1d99
min_M_RTI_factors_full_off = 1d99

alpha_RTI_src_max_q#

alpha_RTI_src_min_q#

option to set alpha_RTI source term to zero when cell q out of bounds. to turn off RTI near surface or center

alpha_RTI_src_max_q = 1d0
alpha_RTI_src_min_q = 0d0

alpha_RTI_src_min_v_div_cs#

option to set alpha_RTI source term to zero when v/cs < this min. e.g. to filter out false sources ahead of shock

alpha_RTI_src_min_v_div_cs = 1d0

Radial Stellar Pulsations (RSP)#

inspired by Radec Smolec’s Program

must set mass, Teff, L, X, and Z.

RSP_mass = -1
RSP_Teff = -1
RSP_L = -1
RSP_X = -1
RSP_Z = -1

Parameters of the convection model. Note that RSP_alfap, RSP_alfas, RSP_alfac, RSP_alfad and RSP_gammar are expressed in the units of standard values. Standard values are the ones for which static version of the Kuhfuss model reduces to standard MLT. See Table 1 in Smolec & Moskalik (2008) for standard values and the description of the parameters.

RSP_alfa = 1.5d0            ! mixing length; alfa = 0 gives a purely radiative model.
RSP_alfac = 1.0d0           ! convective flux; Lc ~ RSP_alfac
RSP_alfas = 1.0d0           ! turbulent source; Lc ~ 1/ALFAS; PII ~ RSP_alfas
RSP_alfad = 1.0d0           ! turbulent dissipation; damp ~ RSP_alfad
RSP_alfap = 0.0d0           ! turbulent pressure; Pt ~ alfap
RSP_alfat = 0.0d0           ! turbulent flux; Lt ~ RSP_alfat; overshooting.
RSP_alfam = 0.25d0          ! eddy viscosity; Chi & Eq ~ RSP_alfam
RSP_gammar = 0.0d0          ! radiative losses; dampR ~ RSP_gammar

time weighting for end-of-step vs start-of-step values in equations. 1 corresponds to fully implicit scheme - stable, but can have large numerical damping. 0.5 corresponds to trapezoidal rule integration - gives least numerical damping. do not use values less than 0.5. strongly recommend 0.5 for theta and thetat. don’t mess with any of these unless you know what you are doing or like to watch the code crash.

RSP_theta = 0.5d0           ! Pgas and Prad
RSP_thetat = 0.5d0          ! Pturb
RSP_thetae = 0.5d0          ! erad in terms using f_Edd
RSP_thetaq = 1.0d0          ! Pvsc
RSP_thetau = 1.0d0          ! Eq and Uq
RSP_wtr = 0.6667d0          ! Lr
RSP_wtc = 0.6667d0          ! Lc
RSP_wtt = 0.6667d0          ! Lt
RSP_gam = 1.0d0             ! Et src_snk

controls for building the initial model

RSP_nz = 150
RSP_nz_outer = 40
RSP_T_anchor = 11d3
RSP_T_inner = 2d6
RSP_testing = .false.
RSP_dq_1_factor = 2d0
RSP_max_outer_dm_tries = 100
RSP_max_inner_scale_tries = 100
RSP_T_anchor_tolerance = 1d-8

allowed relative difference between T at base of outer region and T_anchor if fail trying to create initial model, try increasing this to 1d-6 or more

RSP_T_inner_tolerance = 1d-8

allowed relative difference between T at inner boundary and T_inner if fail trying to create initial model, try increasing this to 1d-6 or more

RSP_relax_initial_model = .true.
RSP_relax_alfap_before_alfat = .true.
RSP_relax_max_tries = 1000
RSP_relax_dm_tolerance = 1d-6

Initial kick makes use of the scaled linear velocity eigenvector of a given mode or of the linear combination of the eigenvectors for the fundamental mode and first two radial overtones. The surface velocity is set to RSP_kick_vsurf_km_per_sec and the mode content is set by RSP_fraction_1st_overtone and RSP_fraction_2nd_overtone

RSP_kick_vsurf_km_per_sec = 0.1d0
RSP_fraction_1st_overtone = 0d0
RSP_fraction_2nd_overtone = 0d0

fraction from fundamental = 1d0 - (1st + 2nd) Note: This is important for models in which two or more modes are linearly unstable. Appropriate setting may help to arrive at the desired mode, since the final pulsation state may depend on initial conditions set by the three parameters above. Integration of the same model with different initial kicks is a way to study the nonlinear mode selection - for an example see Fig. 1 in Smolec & Moskalik (2010).

random initial velocity profile. added to any kick from eigenvector.

RSP_Avel = 0d0
RSP_Arnd = 0d0

period controls

RSP_target_steps_per_cycle = 600
RSP_min_max_R_for_periods = -1
RSP_min_deltaR_for_periods = -1
RSP_min_PERIOD_div_PERIODLIN = 0.5d0
RSP_mode_for_setting_PERIODLIN = 1
RSP_default_PERIODLIN = 34560

when to stop

RSP_max_num_periods = -1
RSP_GREKM_avg_abs_frac_new = 0.1d0
RSP_GREKM_avg_abs_limit = -1

timestep limiting

RSP_initial_dt_factor = 1d-2

start with smaller timestep to give time for initial model to adjust

RSP_v_div_cs_threshold_for_dt_limit = 0.8d0
RSP_max_dt_times_min_dr_div_cs = 2d0

i.e., make dt <= this times min sound crossing time dr/cs only considering cells with abs(v)/cs > threshold

RSP_max_dt_times_min_rad_diff_time = -1d0

make dt <= min time for radiative diffusion for RHD

RSP_max_dt = -1
RSP_report_limit_dt = .false.

artificial viscosity controls for the equations see: Appendix C in Stellingwerf (1975). In principle, for not too-non-adiabatic convective models artificial viscosity is not needed or should be very small. Hence a large cut-off parameter below (in purely radiative models the default value for cut-off was 0.01)

RSP_cq = 4.0d0
RSP_zsh = 0.1d0

zsh > 0 delays onset of artificial viscosity can eliminate most/all interior dissipation while still providing for extreme cases. using this parameter the dependence of limiting amplitude on cq is very weak. for Tscharnuter & Winkler form of artificial viscosity

RSP_Qvisc_linear = 0d0
RSP_Qvisc_quadratic = 0d0

as described in section 4.2 of mesa3, 2015. RSP_Qvisc_linear is analogous to shock_spread_linear RSP_Qvisc_quadratic is analogous to shock_spread_quadratic if switch to this form, set RSP_cq = 0 to shut off the Neumann & Richtmyer form. note that this form also uses RSP_zsh to delay onset of artificial viscosity

surface pressure. provides outer boundary condition for momentum equation.

RSP_use_Prad_for_Psurf = .false.
RSP_use_atm_grey_with_kap_for_Psurf = .false.
RSP_tau_surf_for_atm_grey_with_kap = 3d-3
RSP_fixed_Psurf = .true.
RSP_Psurf = 0d0
set_RSP_Psurf_to_multiple_of_initial_P1 = -1

set RSP_Psurf to this times initial surface cell pressure

RSP_surface_tau = 0.001d0

solver controls

RSP_tol_max_corr = 1d-8
RSP_tol_max_resid = 1d-6
RSP_max_iters_per_try = 100
RSP_max_retries_per_step = 50
RSP_report_undercorrections = .false.
RSP_nz_div_IBOTOM = 30d0
RSP_min_tau_for_turbulent_flux = 2d2

output data for work integrals during a particular period

RSP_work_period = -1
RSP_work_filename = 'work.data'

output data for 3d map. format same as for gnuplot pm3d

RSP_write_map = .false.
RSP_map_columns_filename = 'map_columns.list'

items listed in your map columns must also appear in your profile columns

RSP_map_filename = 'map.data'
RSP_map_first_period = -1
RSP_map_last_period = -1
RSP_map_zone_interval = 2
RSP_map_history_filename = 'map_history.data'

rsp hooks

use_other_RSP_linear_analysis = .false.
use_other_RSP_build_model = .false.

for special tests can set ALFA = 0 for pure radiative with no turbulence or convection. can set zero_gravity = .true. can set opacity to be constant times density.

RSP_kap_density_factor = -1

else set opacity to this times density

rsp misc

RSP_efl0 = 1.0d2
RSP_nmodes = 3
RSP_trace_RSP_build_model = .false.
use_RSP_new_start_scheme = .false.
RSP_Qvisc_linear_static = 0d0
RSP_relax_adjust_inner_mass_distribution = .true.
RSP_do_check_omega = .true.
RSP_report_check_omega_changes = .false.
RSP_hydro_only = .false.

this does not work with the standard build model, so requires use_other_RSP_build_model

mixing misc#

such as smoothing and editing of diffusion coefficients

mix_factor#

Mixing coefficients are multiplied by this factor. The mix_factor is applied in subroutine get_convection_sigmas in star/private/mix_info.f90 – the lagrangian diffusion coefficient sigma(k) at cell boundary k is set to mix_factor*D*(4*pi*r(k)^2*rho_face(k))^2. Note that the value of D is not changed – it is just used as a term in calculating sigma.

mix_factor = 1

max_conv_vel_div_csound#

convective velocities are limited to local sound speed times this factor.

max_conv_vel_div_csound = 1d99

max_conv_vel_div_csound_maxq#

only consider max_conv_vel_div_csound from center out to this location

max_conv_vel_div_csound_maxq = -1

max_v_for_convection#

disable convection for locations with v > than this limit. In km/s.

max_v_for_convection = 1d99

max_q_for_convection_with_hydro_on#

disable convection for locations with q > than this limit when either v_flag or u_flag are true.

max_q_for_convection_with_hydro_on = 1d99

max_v_div_cs_for_convection#

disable convection for locations with abs(v)/cs > this limit

max_v_div_cs_for_convection = 1d99

max_abs_du_div_cs_for_convection#

main purpose is to force radiative in shock face

max_abs_du_div_cs_for_convection = 0.03d0

prune_bad_cz_min_Hp_height#

Lower limit on radial extent of cz (<= 0 to disable). In units of average pressure scale height at top and bottom of region. Remove tiny convection zones unless have strong nuclear burning This allows emergence of very small cz at site of he core flash, for example. i.e., remove if size < prune_bad_cz_min_Hp_height .and. max_log_eps < prune_bad_cz_min_log_eps_nuc.

prune_bad_cz_min_Hp_height = 0

prune_bad_cz_min_log_eps_nuc#

Lower limit on max log eps nuc in cz. remove if size < prune_bad_cz_min_Hp_height .and. max_log_eps < prune_bad_cz_min_log_eps_nuc.

prune_bad_cz_min_log_eps_nuc = -99

redo_conv_for_dr_lt_mixing_length#

Check for small convection zones with total height less than mixing length and redo with reduced mixing_length_alpha to make mixing_length <= dr.

redo_conv_for_dr_lt_mixing_length = .false.

smooth_convective_bdy has been deleted.#

remove_mixing_glitches#

If true, then okay to remove gaps and singletons.

remove_mixing_glitches = .true.

glitches#

The following controls are for different kinds of “glitches” that can be removed.

okay_to_remove_mixing_singleton#

If true, remove singletons.

okay_to_remove_mixing_singleton = .true.

clip_D_limit#

Zero mixing diffusion coeffs that are smaller than this.

clip_D_limit = 0

min_convective_gap#

Close gap between convective regions if smaller than this (< 0 means skip this). Gap measured radially in units of pressure scale height.

min_convective_gap = -1

min_thermohaline_gap#

Close gap between thermohaline mixing regions if smaller than this (< 0 means skip this). Gap measured radially in units of pressure scale height.

min_thermohaline_gap = -1

min_thermohaline_dropout#

max_dropout_gradL_sub_grada#

If find radiative region embedded in thermohaline, and max(gradL - grada) in region is everywhere < max_dropout_gradL_sub_grada and region height is < min_thermohaline_dropout then convert the region to thermohaline. min_thermohaline_dropout <= 0 disables.

min_thermohaline_dropout = -1
max_dropout_gradL_sub_grada = 1d-3

min_semiconvection_gap#

Close gap between semiconvective mixing regions if smaller than this (< 0 means skip this). Gap measured radially in units of pressure scale height.

min_semiconvection_gap = -1

remove_embedded_semiconvection#

If have a semiconvection region bounded on each side by convection, convert it to be convective too.

remove_embedded_semiconvection = .false.

recalc_mix_info_after_evolve#

Re-evaluate mixing info after each evolve step. This is helpful if you want the profiles to reflect the mixing params after the step; otherwise, they give the mixing info from the start of the step (i.e., one step out-of-date)

recalc_mix_info_after_evolve = .false.

set_min_D_mix#

mass_lower_limit_for_min_D_mix#

mass_upper_limit_for_min_D_mix#

min_D_mix#

D_mix will be at least this large if set_min_D_mix is true. doesn’t apply for mass < lower limit or mass > upper limit.

set_min_D_mix = .false.
mass_lower_limit_for_min_D_mix = 0d0
mass_upper_limit_for_min_D_mix = 1d99
min_D_mix = 1d3

set_min_D_mix_in_H_He#

min_D_mix_in_H_He#

D_mix will be at least this large in regions where max mass fractions of H and He add to more that 0.5 if set_min_D_mix_in_H_He is true.

set_min_D_mix_in_H_He = .false.
min_D_mix_in_H_He = 1d3

set_min_D_mix_below_Tmax#

min_D_mix_below_Tmax#

D_mix will be at least this large for cells below location of max temperature if set_min_D_mix_below_Tmax is true.

set_min_D_mix_below_Tmax = .false.
min_D_mix_below_Tmax = 1d3

min_center_Ye_for_min_D_mix#

min_D_mix is only used when center_ye >= this i.e., when center_ye drops below this, min_D_mix = 0.

min_center_Ye_for_min_D_mix = 0.47d0

smooth_outer_xa_big#

smooth_outer_xa_small#

Soften composition jumps in outer layers. If smooth_outer_xa_big and smooth_outer_xa_small are bigger than 0, then starting from the outermost grid point, homogeneously mix a region of size smooth_outer_xa_small (in solar masses), and proceed inwards, linearly reducing the size of the homogeneously mixed region in such a way that it becomes zero. After going smooth_outer_xa_big solar masses in. In this way, the outer smooth_outer_xa_big solar masses are “cleaned” of composition jumps.

smooth_outer_xa_big = -1d0
smooth_outer_xa_small = -1d0

rotation controls#

In the following “am” stands for “angular momentum”.

the mesa implementation of rotation closely follows these papers:

  • Heger, Langer, & Woosley, ApJ, 528, 368. 2000

  • Heger, Woosley, & Spruit, ApJ, 626, 350. 2005

  • D_DSI = dynamical shear instability

  • D_SH = Solberg-Hoiland

  • D_SSI = secular shear instability

  • D_ES = Eddington-Sweet circulation

  • D_GSF = Goldreich-Schubert-Fricke

  • D_ST = Spruit-Tayler dynamo

skip_rotation_in_convection_zones#

if true, then set rotational diffusion coefficients to 0 in convective regions. This applies both for material mixing and diffusion of angular momentum.

skip_rotation_in_convection_zones = .false.

am_D_mix_factor#

Rotation and mixing of material. D_mix = diffusion coefficient for mixing of material. It is sum of non-rotational and rotational components. The rotational part is multiplied by this factor.

D_mix = D_mix_non_rotation + f*am_D_mix_factor*(
    D_DSI_factor * D_DSI +
    D_SH_factor  * D_SH +
    D_SSI_factor * D_SSI +
    D_ES_factor  * D_ES +
    D_GSF_factor * D_GSF +
    D_ST_factor  * D_ST)
f = 1  when logT <= D_mix_rotation_max_logT_full_on = full_on
  = 0  when logT >= D_mix_rotation_min_logT_full_on = full_off
  = (log(T)-full_on)/(full_off-full_on) else

note that for regions with brunt N^2 < 0, we set Richardson number to 1 which is > Ri_critical and therefore turns off DSI and SSI

according to Heger et al 2000 : 1/30d0 by default : 0

am_D_mix_factor = 0

am_nu_factor#

am_nu_non_rotation_factor#

diffusion of angular momentum

am_nu = diffusion coefficient for angular momentum

am_nu_non_rot = am_nu_factor*am_nu_non_rotation_factor*D_mix_non_rotation
am_nu_rot = am_nu_factor*(
    am_nu_visc_factor* D_visc +
    am_nu_DSI_factor * D_DSI +
    am_nu_SH_factor  * D_SH +
    am_nu_SSI_factor * D_SSI +
    am_nu_ES_factor  * D_ES +
    am_nu_GSF_factor * D_GSF +
    am_nu_ST_factor  * nu_ST)
am_nu = am_nu_non_rot + am_nu_rot

Note that for regions with brunt N^2 < 0, we set Richardson number to 1 which is > Ri_critical and therefore turns off DSI and SSI.

see also star_job controls for am_nu_rot_flag

am_nu_factor = 1
am_nu_non_rotation_factor = 1

am_nu_DSI_factor#

< 0 means use D_DSI_factor

am_nu_DSI_factor = -1

am_nu_SSI_factor#

< 0 means use D_SSI_factor

am_nu_SSI_factor = -1

am_nu_SH_factor#

< 0 means use D_SH_factor

am_nu_SH_factor = -1

am_nu_ES_factor#

< 0 means use D_ES_factor

am_nu_ES_factor = -1

am_nu_GSF_factor#

< 0 means use D_GSF_factor

am_nu_GSF_factor = -1

am_nu_ST_factor#

< 0 means use D_ST_factor

am_nu_ST_factor = -1

am_nu_visc_factor#

< 0 means use D_visc_factor. By default = 1 to mix angular momentum.

am_nu_visc_factor = 1

am_nu_omega_rot_factor#

am_nu_omega_non_rot_factor#

dj/dt = d/dm((4 pi r^2 rho)^2*(am_nu_omega*i_rot*domega/dm + am_nu_j*dj/dm)) am_nu_omega = am_nu_omega_non_rot_factor*am_nu_non_rot + am_nu_omega_rot_factor*am_nu_rot

am_nu_omega_rot_factor = 1
am_nu_omega_non_rot_factor = 1

am_nu_j_rot_factor#

am_nu_j_non_rot_factor#

dj/dt = d/dm((4 pi r^2 rho)^2*(am_nu_omega*i_rot*domega/dm + am_nu_j*dj/dm)) am_nu_j = am_nu_j_non_rot_factor*am_nu_non_rot + am_nu_j_rot_factor*am_nu_rot

am_nu_j_rot_factor = 0
am_nu_j_non_rot_factor = 0

set_uniform_am_nu_non_rot#

uniform_am_nu_non_rot#

You can specify a uniform value for am_nu_non_rot by setting this flag true. A large uniform am_nu will produce a uniform omega.

set_uniform_am_nu_non_rot = .false.
uniform_am_nu_non_rot = 1d20

set_min_am_nu_non_rot#

min_am_nu_non_rot#

You can also specify a minimum am_nu_non_rot. am_nu will be at least this large.

set_min_am_nu_non_rot = .false.
min_am_nu_non_rot = 1d8

min_center_Ye_for_min_am_nu_non_rot#

min_am_nu_non_rot is only used when center Ye >= this.

min_center_Ye_for_min_am_nu_non_rot = 0.47d0

Each rotationally induced diffusion coefficient has a factor that lets you control it. Value of 1 gives normal strength; value of 0 turns it off.

Note that for regions with brunt N^2 < 0, we set Richardson number to 1, which is > Ri_critical and therefore turns off DSI and SSI.

D_DSI_factor = 0
D_SH_factor = 0
D_SSI_factor = 0
D_ES_factor = 0
D_GSF_factor = 0
D_ST_factor = 0

D_visc_factor#

Kinematic shear viscosity. Should be = 0 because viscosity doesn’t mix chemical elements.

D_visc_factor = 0

am_gradmu_factor#

Sensitivity to composition gradients. In calculation of rotational induced mixing, grad_mu is multiplied by am_gradmu_factor. Value from from Heger et al 2000.

am_gradmu_factor = 0.05d0

Spatial smoothing is used in calculations of diffusion coefficients. These control the smoothing window widths (number of cells on each side).

smooth_D_DSI = 0
smooth_D_SH = 0
smooth_D_SSI = 0
smooth_D_ES = 0
smooth_D_GSF = 0
smooth_D_ST = 0
smooth_nu_ST = 0
smooth_D_omega = 0
smooth_am_nu_rot = 0

ST_angsmt#

ST_angsml#

Temporal smoothing of ST coefficients. See rotation_mix_info.f90 for details

ST_angsmt = 0.2d0
ST_angsml = 1d-3

simple_i_rot_flag#

If true, i_rot = (2/3)*r^2. If false, use slightly more complex expression that takes into account finite shell thickness. In practice, there doesn’t seem to be a significant difference.

simple_i_rot_flag = .false.

do_adjust_J_lost#

adjust_J_fraction#

adjust angular momentum With do_adjust_J_lost = .false., the angular momentum removed via winds from the star corresponds to that contained in the removed layers. However, since j_rot can increase steeply in the very outer layers, very small steps are required to obtain a convergent solution. To avoid this, the do_adjust_J_lost option adjusts the angular momentum content of layers below those removed, such that

actual_J_lost = &
   adjust_J_fraction*mass_lost*s% j_rot_surf + &
   (1d0 - adjust_J_fraction)*s% angular_momentum_removed

where s% angular_momentum_removed is the angular momentum contained in the removed layers of the star in that step. Note that s% angular_momentum_removed is set to actual_J_lost after this.

The region from which angular momentum is removed is chosen such that at its bottom q<min_q_for_adjust_J_lost, it contains at least min_J_div_delta_J times the angular momentum that needs to be accounted for. Angular momentum in these regions is adjusted in such a way that no artificial shear is produced at the inner boundary.

This can also be used to model mass loss mechanisms that remove more angular momentum than mass_lost*s% j_rot_surf, for instance magnetic braking or wind mass loss. In that case, you can use the use_other_j_for_adjust_J_lost option to specify a specific angular momentum of removed material different from j_rot_surf

In order to prevent the algorithm from digging to deep to adjust J, there is a timestep limit adjust_J_q_limit

do_adjust_J_lost = .true.
adjust_J_fraction = 1d0
min_q_for_adjust_J_lost = 0.995d0
min_J_div_delta_J = 3d0

premix_omega#

if premix_omega is true, then do 1/2 of the transport of angular momentum before updating the structure and 1/2 after. otherwise, do all of the transport after updating the structure. RECOMMENDED to turn it on when modelling an accreting star or when using do_adjust_J_lost.

premix_omega = .true.

angular_momentum_error_warn#

angular_momentum_error_stop#

if the relative change in total angular momentum exceeds these values, then a warning is given on the terminal output, or the simulation is stopped altogether. Not applied when using other_torque routines or for binaries.

angular_momentum_error_warn = 5d-6
angular_momentum_error_retry = 1d-6

recalc_mixing_info_each_substep#

if recalc_mixing_info_each_substep is true, then recalculate the omega mixing coefficients after each substep of the solve omega mix process.

recalc_mixing_info_each_substep = .false.

w_div_wcrit_max#

When fitted_fp_ft_i_rot = .true., limit fp and ft to their values at this w_div_wcrit

w_div_wcrit_max = 0.89d0

w_div_wcrit_max2#

When w_div_wc_flag is true, rather than a hard limit on w_div_wcrit we use two limiting values w_div_wcrit_max < w_div_wcrit_max2 to provide a smooth transition. Nothing is done when w_div_wc < w_div_wcrit_max, if w_div_wcrit_max < w_div_wc < w_div_wcrit_max2, we apply a sigmoid, and in the limit of j_rot -> infinity, the resulting w_div_wc will match w_div_wcrit_max2 (being on the top of the sigmoid)

w_div_wcrit_max2 = 0.90d0

D_mix_rotation_max_logT_full_on#

Use rotational components of D_mix for locations where logT <= this. For numerical stability, turn off rotational part of D_mix at very high T.

D_mix_rotation_max_logT_full_on = 9.4d0

D_mix_rotation_min_logT_full_off#

Drop rotational components of D_mix for locations where logT >= this. For numerical stability, turn off rotational part of D_mix at very high T.

D_mix_rotation_min_logT_full_off = 9.5d0

D_mix_rotation_min_tau_full_off#

Drop rotational components of D_mix for locations where tau <= this. For numerical stability, turn off rotational part of D_mix at very low tau.

D_mix_rotation_min_tau_full_off = 0d0

D_mix_rotation_max_tau_full_on#

Use rotational components of D_mix for locations where tau >= this. For numerical stability, turn off rotational part of D_mix at very low tau.

D_mix_rotation_min_tau_full_on = 0d0

D_omega_max_replacement_fraction#

D_omega_growth_rate#

D_omega_mixing_rate#

D_omega_mixing_across_convection_boundary (previously called D_omega_mixing_in_convection_regions)#

D_omega_mixing_rate = 1d0
D_omega_mixing_across_convection_boundary = .false.
max_q_for_D_omega_zero_in_convection_region = 0.8d0

nu_omega_max_replacement_fraction#

nu_omega_growth_rate#

nu_omega_mixing_rate#

nu_omega_mixing_across_convection_boundary#

nu_omega_mixing_rate = 1d0
nu_omega_mixing_across_convection_boundary = .false.
max_q_for_nu_omega_zero_in_convection_region = 0.8d0

atmosphere boundary conditions#

atm_option#

Controls how the surface temperature Tsurf and pressure Psurf are evaluated when setting up outer boundary conditions. We caution that the use of 'fixed_' atmosphere options might conflict with mlt_option = TDC.

  • 'T_tau': set Tsurf and Psurf by solving for the atmosphere structure given a T(tau) relation. The choice of relation is set by the atm_T_tau_relation control. See also the atm_T_tau_opacity, atm_T_tau_errtol, atm_T_tau_max_tries and atm_T_tau_max_steps controls.

  • 'table': set Tsurf and Psurf by interpolating in pre-calculated tables based on model atmospheres. The choice of table is set by the atm_table control. Requires tau_factor = 1, as surface of the model must always attach at the base of the tables, so there is no flexibility to move model surface inward or outward. Note that tau_base = tau_surf is the location at which the model attaches to the table BCs, and there is no particular location identified as the photosphere, so we fall back to the surface values of L, R, and m to calculate quantities such as Teff and log_g. This is consistent with the assumptions used for table construction: geometrically thin atmospheres with constant flux.

  • 'irradiated_grey': set Tsurf by solving for the atmosphere structure given the irradiated-grey T(tau) relation of Guillot, T, and Havel, M., A&A 527, A20 (2011). See also the atm_irradiated_opacity, atm_irradiated_errtol, atm_irradiated_T_eq, atm_irradiated_T_eq, atm_irradiated_kap_v, atm_irradiated_kap_v_div_kap_th, atm_irradiated_P_surf and atm_irradiated_max_tries controls.

  • 'fixed_Teff' :

    set Tsurf from Eddington T(tau) relation for current surface tau and Teff = atm_fixed_Teff. set Psurf = Radiation_Pressure(Tsurf)

  • 'fixed_Tsurf' :

    get value of Tsurf from control parameter atm_fixed_Tsurf. set Teff from Eddington T(tau) relation for given Tsurf and tau=2/3 set Psurf = Radiation_Pressure(Tsurf)

  • 'fixed_Psurf' :

    get value of Psurf from control parameter atm_fixed_Psurf. set Tsurf from L and R using L = 4*pi*R^2*boltz_sigma*T^4. set Teff using Eddington T(tau) relation for tau=2/3 and T=Tsurf.

  • 'fixed_Psurf_and_Tsurf' :

    get value of Psurf from control parameter atm_fixed_Psurf. get value of Tsurf from control parameter atm_fixed_Tsurf. see the conductive_flame test_suite for an example of this boundary condition implemented via the other_surface_PT hook.

atm_option = 'T_tau'

atm_off_table_option#

If have selected 'table' for atm_option, fallback to using this if the args are off the table. Possible choices are 'T_tau' or blank (in which case the code will halt when it encounters an off-table arg)

atm_off_table_option = 'T_tau'

atm_fixed_Teff#

Set this when using atm_option = 'fixed_Teff'

atm_fixed_Teff = 0

atm_fixed_Tsurf#

Set this when using atm_option = 'fixed_Tsurf'

atm_fixed_Tsurf = 0

atm_fixed_Psurf#

Set this when using atm_option = 'fixed_Psurf'

atm_fixed_Psurf = -1

Pextra_factor#

Parameter for extra pressure in surface boundary conditions. Pressure at optical depth tau is calculated as P = tau*g/kap*(1 + Pextra) Pextra takes into account nonzero radiation pressure at tau=0. The equation for Pextra includes Pextra_factor

Pextra = Pextra_factor*(kap/tau)*(L/M)/(6d0*pi*clight*cgrav)

For certain situations such super eddington L, you may need to increase Pextra to help convergence. e.g. try Pextra_factor = 2

Note that Pextra_factor is only applied when atm_option = 'T_tau' and atm_T_tau_opacity = 'fixed' or ‘iterated'.

Pextra_factor = 1

atm_T_tau_relation#

The T(tau) relation to use when atm_option = 'T_tau'

  • 'Eddington': use the grey Eddington T(tau) relation.

  • 'solar_Hopf': use a grey T(tau) relation with an approximate Hopf function tuned to solar data. See Paper II, Sec. A.5. Equivalent to the fit given by Sonoi et al. (2019, A&A, 621, 84) to the Vernazza et al. (1981) VAL-C model.

  • 'Krishna_Swamy': use the grey T(tau) relation described by K.S. Krishna-Swamy, ApJ 145, 174–194 (1966).

  • 'Trampedach_solar': use the analytic fit by Ball (2021, RNAAS 5, 7) to the Hopf function for the solar simulation by Trampedach et al. (2014, MNRAS 442, 805–820)

atm_T_tau_relation = 'Eddington'

atm_T_tau_opacity#

Controls how opacities are calculated throughout the atmosphere when atm_option = 'T_tau'

  • 'fixed': use a uniform opacity, fixed to the opacity of the outermost cell of the interior model

  • 'iterated': use a uniform opacity, iterated to be consistent with the final Tsurf and Psurf at the base of the atmosphere.

  • 'varying': use a varying opacity consistent with the local T and P throughout the atmosphere. Involves numerical integration of the hydrostatic equilibrium equation.

atm_T_tau_opacity = 'fixed'

atm_T_tau_errtol#

Error tolerance for iterations and integrations when atm_option = 'T_tau' and atm_T_tau_opacity = 'iterated' or 'varying'.

atm_T_tau_errtol = 1d-7

atm_T_tau_max_iters#

Maximum number of iterations for the opacity when atm_option = 'T_tau' and atm_T_tau_opacity = 'iterated'.

atm_T_tau_max_iters = 50

atm_T_tau_max_steps#

Maximum number of steps for integrating the hydrostatic equilibrium equation when atm_option = 'T_tau' and atm_T_tau_opacity = 'varying'.

atm_T_tau_max_steps = 500

atm_table#

Determines which table Tsurf and Psurf are interpolated in

  • 'tau_100', 'tau_10', 'tau_1', 'tau_1m1': use model atmosphere tables for Pgas and T at tau=100, 10, 1 or 0.1, respectively; solar Z only, as described in MESA Paper I (Paxton et al. 2011), Sec. 5.3. these tables are primarily for the evolution of low-mass stars, brown dwarfs, and giant planets. They are constructed from Castelli & Kurucz (2003) for Teff > 3000 K and the COND model atmospheres (Allard et al. 2001) for Teff < 3000 K. where no published results are available, the table has been filled in using integration of the Eddington T(tau) relation

  • 'photosphere': use model atmosphere tables for photosphere; range of Z’s, as described in MESA Paper I (Paxton et al. 2011), Sec. 5.3. the tables cover log(Z/Zsun) from -4 to 0.5 for a GN93 solar mixture, and span logg = -0.5 to 5.5 in steps of 0.5 dex and Teff = 2000 to 50 000 K in steps of 250 K. they are constructed, in order of precedence, using the PHOENIX model atmospheres (Hauschildt et al. 1999a,b) and the models by Castelli & Kurucz (2003). where neither is available, an entry is generated by integrating the Eddington T(tau) relation

  • 'WD_tau_25': hydrogen atmosphere tables for cool white dwarfs giving Pgas and T at log10(tau) = 1.4 (tau = 25.11886) Teff goes from 40,000 K down to 2,000K with step of 100 K Log10(g) goes from 9.5 down to 5.5 with step of 0.1. R.D. Rohrmann, L.G. Althaus, and S.O. Kepler, Lyman α wing absorption in cool white dwarf stars, Mon. Not. R. Astron. Soc. 411, 781–791 (2011)

  • 'DB_WD_tau_25': helium dominated (log(H/He)=-5.0) atmosphere tables for DB white dwarfs, provided by Odette Toloza and Detlev Koester. 5000K < Teff < 40000K

atm_table = 'tau_100'

atm_irradiated_opacity#

Controls how thermal opacities are calculated throughout the atmosphere when atm_option = 'irradiated_grey'

  • 'fixed': use a uniform opacity, fixed to the opacity of the outermost cell of the interior model.

  • 'iterated': use a uniform opacity, iterated to be consistent with the final Tsurf and Psurf at the base of the atmosphere.

atm_irradiated_opacity = 'fixed'

atm_irradiated_errtol#

Error tolerance for iterations when atm_option = 'irradiated_grey' and atm_irradiated_opacity = 'iterated'.

atm_irradiated_errtol = 1d-7

atm_irradiated_max_iters#

Maximum number of iterations for the opacity when atm_option = 'irradiated_grey' and atm_irradiated_opacity = 'iterated'.

atm_irradiated_max_iters = 50

atm_irradiated_T_eq#

Equilibrium temperature based on irradiation.

irrad_flux = Lstar/(4*pi*orbit**2)
  • Area of planet in plane perpendicular to irrad_flux = pi*Rplanet**2.

  • Stellar luminosity received by planet = irrad_flux*area.

  • This luminosity determines T_eq: T_eq**4 = irrad_flux/(4*sigma).

atm_irradiated_T_eq = 100

atm_irradiated_kap_v_div_kap_th#

atm_irradiated_kap_v#

The T(tau) relation when atm_option = 'irradiated_grey' depends on the ratio kap_v/kap_th where kap_v is the planet atmosphere opacity for stellar irradiation, and kap_th is the thermal opacity for the internally produced radiation. There are two ways to calculate the ratio:

  • if atm_irradiated_kap_v_div_kap_th > 0, then use it for kap_v/kap_th

  • if atm_irradiated_kap_v_div_kap_th == 0, then use atm_irradiated_kap_v to set kap_v, and then evaluate the ratio using kap_th

atm_irradiated_kap_v_div_kap_th = 0
atm_irradiated_kap_v = 4d-3

atm_irradiated_P_surf#

Surface pressure when atm_option = 'irradiated_grey'; default is 1 bar in cgs units.

atm_irradiated_P_surf = 1d6

use_compression_outer_BC#

gradient of compression vanishes at surface

see Grott, Chernigovski, Glatzel, 2005. d_dm(d_dm(r^2*v)) = 0 at surface by continuity, this is d_dm(d_dt(1/rho)) = 0 at surface finite volume form is (1/rho(1) - 1/rho_start(1)) = (1/rho(2) - 1/rho_start(2)) this BC determines the density for surface cell.

use_compression_outer_BC = .false.

use_momentum_outer_BC#

use P_surf from atm to set pressure gradient at surface in momentum equation calculate v(1) based on pressure difference P_surf - P(1)

use_momentum_outer_BC = .false.

use_zero_Pgas_outer_BC#

use Psurf = Radiation_Pressure(T_start(1))

use_zero_Pgas_outer_BC = .false.

use_fixed_vsurf_outer_BC#

fixed_vsurf#

v at outer boundary of model is set to be fixed_vsurf

use_fixed_vsurf_outer_BC = .false.
fixed_vsurf = 0

use_fixed_Psurf_outer_BC#

fixed_Psurf#

P at outer boundary of model is set to be fixed_Psurf

use_fixed_Psurf_outer_BC = .false.
fixed_Psurf = 0

Tsurf_factor#

used when use_momentum_outer_BC T_surf is set to Tsurf_factor*T_black_body(L_surf,R_surf)

Tsurf_factor = 1

irradiation_flux#

column_depth_for_irradiation#

irradiation_flux = 0
column_depth_for_irradiation = -1

atm_build_tau_outer#

atm_build_dlogtau#

atm_build_errtol#

Parameters controlling atmosphere structure building. MESA can evaluate the spatial structure of the atmosphere for the following atm_option choice:

  • 'T_tau'

The atmosphere structure data are appended to the interior model when add_atmosphere_to_pulsation = .true.. They do not affect the surface boundary conditions applied to the interior model.

atm_build_tau_outer specifies the outermost optical depth to include in the atmosphere; atm_build_dlntau specifies the spacing of atmosphere points in (base-10) logarithmic optical depth; and atm_build_errtol specifies the error tolerance for evaluating the structure.

atm_build_tau_outer = 1d-3
atm_build_dlogtau = 0.01
atm_build_errtol = 1d-8

use_T_tau_gradr_factor#

If .true., modify the radiative gradient so that the temperature profile of the optically thin layers follow the T(τ) relation chosen by atm_T_tau_relation.

use_T_tau_gradr_factor = .false.

starspots#

do_starspots#

If .true., switch on impedance of the surface flux due to magnetic pressure from starspots, parameterized in the style of an atmospheric boundary modification. First described by Somers et al. (2015; ApJ). Detailed discussion of this functionality can be found in MESA V.

do_starspots = .false.

fspot#

Filling Factor of starspots. Valid values between 0.0 and 1.0 (no spots to 100% coverage)

fspot = 0d0

xspot#

Temperature contrast between the spotted and unspotted regions. Valid values are between 1.0 (no contribution from magnetic pressure) and 0.5 (half of the total pressure is due to magnetic pressure)

xspot = 1d0

mass gain or loss#

mass_change#

Rate of accretion (Msun/year). Negative for mass loss. This only applies when the wind_scheme = ''.

mass_change = 0d0

mdot_omega_power#

Enhanced mass loss due to rotation as in Heger, Langer, and Woosley, 2000, ApJ, 528:368-396.

Mdot = Mdot_no_rotation/(1 - Osurf/Osurf_crit)^mdot_omega_power

where

Osurf = angular velocity at surface Osurf_crit^2 = (1 - Gamma_edd)*G*M/R^3 Gamma_edd = kappa*L/(4 pi c G M), Eddington factor

Typical value for mdot_omega_power = 0.43.

Set to 0 to disable this feature.

mdot_omega_power = 0.43d0

max_rotational_mdot_boost#

This limits the rotational boost.

max_rotational_mdot_boost = 1d4

max_mdot_jump_for_rotation#

Don’t increase prev mdot by more that this.

max_mdot_jump_for_rotation = 2

lim_trace_rotational_mdot_boost#

Output to terminal if boost > this.

lim_trace_rotational_mdot_boost = 1d99

rotational_mdot_boost_fac#

Increase mdot.

rotational_mdot_boost_fac = 1d5

rotational_mdot_kh_fac#

Kelvin-helmholtz boost.

rotational_mdot_kh_fac = 0.3d0

surf_avg_tau_min#

Use mass avg starting from this optical depth.

surf_avg_tau_min = 1

surf_avg_tau#

Use mass avg down to this optical depth.

surf_avg_tau = 100

hot_wind_scheme#

cool_wind_RGB_scheme#

cool_wind_AGB_scheme#

This section replaces the old “RGB_wind_scheme” and “AGB_wind_scheme” with temperature-dependent hot_wind and cool_wind. You can still use the RGB and AGB wind scheme as before, the functionality remains.

Now you can also select a hot wind scheme that takes effect above some temperature, set by hot_wind_full_on_T. Similarly, the cool wind scheme has temperature controls that set the temperature below which they are relevant (cool_wind_full_on_T).

As before, an empty string ‘’ means no wind.

The wind “eta” values, which are constant scaling factors, have all renamed *_wind_eta -> *_scaling_factor.

Here is an example of how to translate an existing inlist from the old style to the new:

Before

After

RGB_wind_scheme = ‘Reimers’ Reimers_wind_eta = 0.1 AGB_wind_scheme = ‘Blocker’ Blocker_wind_eta = 0.5 RGB_to_AGB_wind_switch = 1d-4

cool_wind_RGB_scheme = ‘Reimers’ Reimers_scaling_factor = 0.1 cool_wind_AGB_scheme = ‘Blocker’ Blocker_scaling_factor = 0.5 RGB_to_AGB_wind_switch = 1d-4

! only use the cool_wind_scheme cool_wind_full_on_T = 1d10 !K hot_wind_full_on_T = 1.1d10 !K hot_wind_scheme = ‘’

suggested hot and cool wind schemes follow but any valid wind option will work for either hot or cool.

Empty string means no wind

Suggested hot wind option:

  • ‘Vink’

  • ‘Bjorklund’

Suggested cool wind options:

  • ‘Reimers’

  • ‘Blöcker’

  • ‘de Jager’

  • ‘van Loon’

  • ‘Nieuwenhuijzen’

For now the ‘Dutch’ scheme can be used in either capacity.

NOTE: for schemes that scale with metallicitity, we use Zbase rather than Z (as long as Zbase > 0). This is because wind mass loss rate is primarily determined by iron opacity, which is unlikely to change during the evolution.

hot_wind_scheme = ''
cool_wind_RGB_scheme = ''
cool_wind_AGB_scheme = ''

cool_wind_full_on_T#

hot_wind_full_on_T#

NOTE: hot_wind_full_on_T was previously called ‘cool_wind_full_off_T’

use only cool wind schemes for T_phot < cool_wind_full_on_T use only hot wind schemes for T_phot > hot_wind_full_on_T if cool_wind_full_on_T /= hot_wind_full_on_T then ramp between these limits requires hot_wind_full_on_T > cool_wind_full_on_T

cool_wind_full_on_T = 0.8d4
hot_wind_full_on_T = 1.2d4

RGB_to_AGB_wind_switch#

If center hydrogen abundance is < 0.01 and center helium abundance by mass is less than RGB_to_AGB_wind_switch, then system will use AGB_wind_scheme rather than RGB_wind_scheme.

RGB_to_AGB_wind_switch = 1d-4

The code will automatically choose between an RGB wind and an AGB wind. The following names for the different schemes are recognized:

  • ‘Reimers’

  • ‘Blocker’

  • ‘de Jager’

  • ‘van Loon’

  • ‘Nieuwenhuijzen’

  • ‘Vink’

  • ‘Dutch’

  • ‘Bjorklund’

  • ‘other’ — for wind options implemented using other_wind hook

Reimers_scaling_factor#

Reimers mass loss for red giants.

D. Reimers “Problems in Stellar Atmospheres and Envelopes” Baschek, Kegel, Traving (eds), Springer, Berlin, 1975, p. 229.

Parameter for mass loss by Reimers wind prescription. Reimers mdot is eta*4d-13*L*R/M (Msun/year), with L, R, and M in solar units. Typical value is 0.5.

Reimers_scaling_factor = 0

Blocker_scaling_factor = 0#

Blocker’s mass loss for AGB stars.

T. Blocker “Stellar evolution of low and intermediate-mass stars” A&A 297, 727-738 (1995).

Parameter for mass loss by Blocker’s wind prescription. Blocker mdot is eta*4.83d-9*M**-2.1*L**2.7*4d-13*L*R/M (Msun/year), with L, R, and M in solar units. Typical value is 0.1d0.

Blocker_scaling_factor = 0

de_Jager_scaling_factor#

de Jager mass loss for various applications. de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259. Parameter for mass loss by de Jager wind prescription.

de_Jager_scaling_factor = 0d0

van_Loon_scaling_factor#

see van Loon et al. 2005, A&A, 438, 273 “An empirical formula for the mass-loss rates of dust-enshrouded red supergiants and oxygen-rich Asymptotic Giant Branch stars”

van_Loon_scaling_factor = 0d0

Nieuwenhuijzen_scaling_factor#

See Nieuwenhuijzen, H.; de Jager, C. 1990, A&A, 231, 134.

Nieuwenhuijzen_scaling_factor = 0d0

Vink_scaling_factor#

Vink, J.S., de Koter, A., & Lamers, H.J.G.L.M., 2001, A&A, 369, 574. “Mass-loss predictions for O and B stars as a function of metallicity”

Vink_scaling_factor = 0d0

Bjorklund_scaling_factor#

Björklund, R., Sundqvist, J.O., Puls, J., & Najarro, F., 2021, A&A, 648, A36. “New predictions for radiation-driven, steady-state mass-loss and wind-momentum from hot, massive stars II. A grid of O-type stars in the Galaxy and the Magellanic Clouds”

This wind scheme does not feature a bistability jump.

Bjorklund_scaling_factor = 0d0

Dutch_scaling_factor#

The “Dutch” wind scheme for massive stars combines results from several papers, all with authors mostly from the Netherlands.

The particular combination we use is based on Glebbeek, E., et al, A&A 497, 255-264 (2009) [more Dutch authors!]

For Teff > 1e4 and surface H > 0.4 by mass, use Vink et al 2001 Vink, J.S., de Koter, A., & Lamers, H.J.G.L.M., 2001, A&A, 369, 574.

For Teff > 1e4 and surface H < 0.4 by mass, use Nugis & Lamers 2000 Nugis, T.,& Lamers, H.J.G.L.M., 2000, A&A, 360, 227 Some folks use 0.8 for non-rotating models (Maeder & Meynet, 2001).

Dutch_scaling_factor = 0d0

Dutch_wind_lowT_scheme#

For Teff < 1e4

Use de Jager if Dutch_wind_lowT_scheme = 'de Jager' de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259.

Use van Loon if Dutch_wind_lowT_scheme = 'van Loon' van Loon et al. 2005, A&A, 438, 273.

Use Nieuwenhuijzen if Dutch_wind_lowT_scheme = 'Nieuwenhuijzen' Nieuwenhuijzen, H.; de Jager, C. 1990, A&A, 231, 134

Dutch_wind_lowT_scheme = 'de Jager'

Kudritzki_scaling_factor#

Radiation driven winds of hot stars. See Kudritzki et al, Astron. Astrophys. 219, 205-218 (1989). this is now implemented using the other_wind hook. see other_physics_hooks test case.

Grafener_scaling_factor#

Grafener, G. & Hamann, W.-R. 2008, A&A 482, 945 contributed to mesa by Nilou Afsari this is now implemented using the other_wind hook. see other_physics_hooks test case.

Stern51_scaling_factor#

this is now implemented using the other_wind hook. see other_physics_hooks test case.

use_accreted_material_j#

Angular momentum of accreted material.

use_accreted_material_j = .false.

If false, then accreted material is given j so that it is rotating at the same angular velocity as the surface. If true, then accreted material is given j = accreted_material_j.

accreted_material_j = 0

no_wind_if_no_rotation#

Use this to delay start of wind until after have started rotation.

no_wind_if_no_rotation = .false.

min_wind#

Min wind in Msun/year > 0; ignore this limit if it is <= 0. e.g., might have low level wind even when normal scheme doesn’t call for any.

min_wind = 0d0

max_wind#

Max wind in Msun/year > 0; ignore this limit if it is <= 0.

max_wind = 0d0

For critical rotation mass loss Redo step as needed to find mdot that brings model to just below critical. if max_mdot_redo_cnt > 0, and surf_omega_div_omega_crit > surf_omega_div_omega_crit_limit, then recompute the step while increasing mdot, until surf_omega_div_omega_crit < surf_omega_div_omega_crit_limit. Once an upper limit for mdot is found, the solution for mdot is further refined by bisection until it is computed to a tolerance of surf_omega_div_omega_crit_tol. During iterations, mdot is adjusted alternately by multiplication by mdot_revise_factor, and by adjusting it by implicit_mdot_boost*mdot_initial, where mdot_initial is the value of mdot at the first iteration. This is done to deal with mass accreting stars, where mdot might need to change sign for the star to remain below critical. This is a direct replacement for surf_w_div_w_crit_limit and surf_w_div_w_crit_tol

max_mdot_redo_cnt = 0
min_years_dt_for_redo_mdot = 0
surf_omega_div_omega_crit_limit = 0.99d0
surf_omega_div_omega_crit_tol = 0.05d0
mdot_revise_factor = 1.1d0
implicit_mdot_boost = 0.1d0

implicit wind computation.#

max_tries_for_implicit_wind#

The implicit method will modify the mass transfer rate and redo the step until it either finds a solution, or the number of tries hits max_tries_for_implicit_wind. If max_tries_for_implicit_wind = 0, the wind computation is explicit, meaning that the value of mdot is set using values at the start of the step. This only applies when mdot < 0.

max_tries_for_implicit_wind = 0

iwind_tolerance#

Tolerance for which a solution is considered valid. A solution is valid if

abs(explicit_mdot - implicit_mdot) <
       abs(implicit_mdot)*iwind_tolerance

where

explicit_mdot = mstar_dot at start of step
implicit_mdot = mstar_dot at end of step
iwind_tolerance = 1d-3

iwind_lambda#

If do not satisfy tolerance, redo with a different mdot as follows:

mstar_dot = explicit_mdot + &

iwind_lambda*(implicit_mdot - explicit_mdot)

iwind_lambda = 1d0

super_eddington_scaling_factor#

For super eddington wind we use Ledd averaged by mass to optical depth tau = surf_avg_tau.

super_eddington_scaling_factor = 0

super_eddington_wind_Ledd_factor#

Parameter for mass loss driven by super Eddington luminosity. Divide L by this factor when computing super Eddington wind, e.g., if this is 2, then only get wind when L/2 > Ledd.

super_eddington_wind_Ledd_factor = 1

wind_boost_full_off_L_div_Ledd#

Boost off for L/Ledd <= this (set large to disable this). This alternative form is used when super_eddington_scaling_factor == 0.

wind_boost_full_off_L_div_Ledd = 1.5d0

wind_boost_full_on_L_div_Ledd#

Do max boost for L/Ledd >= this. This alternative form is used when super_eddington_scaling_factor == 0.

wind_boost_full_on_L_div_Ledd = 5

super_eddington_wind_max_boost#

Multiply wind mdot by up to this amount. This alternative form is used when super_eddington_scaling_factor == 0.

super_eddington_wind_max_boost = 1

trace_super_eddington_wind_boost#

Send super eddington wind information to terminal.

trace_super_eddington_wind_boost = .false.

mass_change_full_on_dt#

mass_change_full_off_dt#

These params provide the option to turn off mass change when have very small timesteps. Between mass_change_full_on_dt and mass_change_full_off_dt mass change is gradually reduced. Units in seconds.

mass_change_full_on_dt = 1d-99
mass_change_full_off_dt = 1d-99

trace_dt_control_mass_change#

trace_dt_control_mass_change = .false.

max_star_mass_for_gain#

Automatic stops for mass loss/gain in Msun units (negative means ignore this parameter). Turn off mass gain when star mass reaches this limit.

max_star_mass_for_gain = -1

min_star_mass_for_loss#

Automatic stops for mass loss/gain in Msun units (negative means ignore this parameter). Turn off mass loss when star mass reaches this limit.

min_star_mass_for_loss = -1

max_T_center_for_any_mass_loss#

No mass loss for T center > this.

max_T_center_for_any_mass_loss = 2d9

max_T_center_for_full_mass_loss#

No reduction in mass loss for T center <= this. This must be <= max_T_center_for_full_mass_loss. Reduce mass loss rate to 0 as T center climbs from max_for_full to max_for_any. The idea behind this is that during final stages of burning, there is so little time left in the life of the star, that any mass loss to winds will be negligible, but the inclusion of that insignificant mass loss can actually make convergence more difficult, so you are better off without it.

max_T_center_for_full_mass_loss = 1d9

wind_H_envelope_limit#

Winds automatically shut off when star_mass - he_core_mass mass is less than this limit. The value of he_core_boundary_h1_fraction defines he_core_mass. Mass in Msun units. Previously called wind_envelope_limit.

wind_H_envelope_limit = -1

wind_H_He_envelope_limit#

Winds automatically shut off when star_mass - co_core_mass is less than this limit. The value of co_core_boundary_he4_fraction defines co_core_mass. Mass in Msun units.

wind_H_He_envelope_limit = -1

wind_He_layer_limit#

Winds automatically shut off when he_core_mass - co_core_mass is less than this limit. Mass in Msun units.

wind_He_layer_limit = -1

rlo_scaling_factor#

Amplitude of mass loss. “rlo” wind scheme provides a simple radius-determined-wind with exponential increase.

rlo_scaling_factor = 0

rlo_wind_min_L#

Only on when L > this limit. (Lsun)

rlo_wind_min_L = 1d-6

rlo_wind_max_Teff#

Only on when Teff < this limit.

rlo_wind_max_Teff = 1d99

rlo_wind_roche_lobe_radius#

Only on when R > this (Rsun).

rlo_wind_roche_lobe_radius = 0.40d0

rlo_wind_base_mdot#

Base rate of mass loss when R = roche lobe radius (Msun/year).

rlo_wind_base_mdot = 1d-3

rlo_wind_scale_height#

Determines exponential growth rate of mass loss (Rsun).

rlo_wind_scale_height = 1d-1

roche_lobe_xfer_full_on#

Full accretion when R/RL <= this. Limit accretion when Roche lobe is nearing full (only with rlo_scaling_factor > 0).

roche_lobe_xfer_full_on = 0.5d0

roche_lobe_xfer_full_off#

No accretion when R/RL >= this.

roche_lobe_xfer_full_off = 1.0d0

controls for adjust_mass#

max_logT_for_k_below_const_q#

max_q_for_k_below_const_q#

min_q_for_k_below_const_q#

Move k_below_const_q inward from surface until q(k) <= max_q. Then continue moving inward until reach logT(k) >= max_logT or q(k) <= min_q.

max_logT_for_k_below_const_q = 5
max_q_for_k_below_const_q = 1.0d0
min_q_for_k_below_const_q = 0.999d0

max_logT_for_k_const_mass#

max_q_for_k_const_mass#

min_q_for_k_const_mass#

Move k_const_mass inward from k_below_const_q+1 until q(k) <= max_q. Then continue moving inward until reach logT(k) >= max_logT or q(k) <= min_q.

max_logT_for_k_const_mass = 6
max_q_for_k_const_mass = 1.0d0
min_q_for_k_const_mass = 0.995d0

composition controls#

accrete_same_as_surface#

If true, composition of accreted material is identical to the current surface composition. If false, then the composition is determined by accrete_given_mass_fractions.

The actual mass accretion rate can be set up using the mass_change option.

accrete_same_as_surface = .true.

accrete_given_mass_fractions#

If true, use accrete_given_mass_fractions, num_accretion_species, accretion_species_id and accretion_species_xa to determine composition of accreted material – they must add to 1.0.

If false, then the composition is determined using accretion_h1, accretion_h2, accretion_he3, accretion_he4 and accretion_zfracs.

The actual mass accretion rate can be set up using the mass_change option.

Note that this control is ignored if accrete_same_as_surface is true.

accrete_given_mass_fractions = .false.

num_accretion_species#

accretion_species_id#

accretion_species_xa#

If accrete_same_as_surface is false and accrete_given_mass_fractions is true, then composition of accreted material is determined by these options. The actual mass accretion rate can be set up using the mass_change option.

num_accretion_species can be up to s% max_num_accretion_species, see star/public/star_def.inc for the value of this parameter.

For each of num_accretion_species, the id for the isotope needs to be specified by accretion_species_id as defined in chem/public/chem_def.f90.

Mass fractions for each isotope are defined by accretion_species_xa

num_accretion_species = 0
accretion_species_id(1) = ''
accretion_species_xa(1) = 0

accretion_h1#

accretion_h2#

accretion_he3#

accretion_he4#

If accrete_same_as_surface is false and accrete_given_mass_fractions is false, then the mass fractions of h1, h2, he3 and h4 are determined by these options. Mass fractions for metals are set with the accretion_zfracs control. The actual mass accretion rate can be set up using the mass_change option.

If no h2 in current net, then it is automatically added to h1.

accretion_h1 = 0
accretion_h2 = 0
accretion_he3 = 0
accretion_he4 = 0

accretion_zfracs =#

One of the following identifiers for different Z fractions from chem_def.

  • AG89_zfracs = 1, Anders & Grevesse 1989

  • GN93_zfracs = 2, Grevesse & Noels 1993

  • GS98_zfracs = 3, Grevesse & Sauval 1998

  • L03_zfracs = 4, Lodders 2003

  • AGS05_zfracs = 5, Asplund, Grevesse & Sauval 2005

or set accretion_zfracs = 0 to use the following list of z fractions

accretion_zfracs = -1

accretion_dump_missing_metals_into_heaviest#

this controls the treatment metals that are not included in the current net. if this flag is true, then the mass fractions of missing metals are added to the mass fraction of the most massive metal included in the net. if this flag is false, then the mass fractions of the metals in the net are renormalized to make up for the total mass fraction of missing metals.

accretion_dump_missing_metals_into_heaviest = .true.

Special list of z fractions. If you use these, they must add to 1.0.

z_fraction_li = 0
z_fraction_be = 0
z_fraction_b = 0
z_fraction_c = 0
z_fraction_n = 0
z_fraction_o = 0
z_fraction_f = 0
z_fraction_ne = 0
z_fraction_na = 0
z_fraction_mg = 0
z_fraction_al = 0
z_fraction_si = 0
z_fraction_p = 0
z_fraction_s = 0
z_fraction_cl = 0
z_fraction_ar = 0
z_fraction_k = 0
z_fraction_ca = 0
z_fraction_sc = 0
z_fraction_ti = 0
z_fraction_v = 0
z_fraction_cr = 0
z_fraction_mn = 0
z_fraction_fe = 0
z_fraction_co = 0
z_fraction_ni = 0
z_fraction_cu = 0
z_fraction_zn = 0

lgT_lo_for_set_new_abundances#

lgT_hi_for_set_new_abundances#

Composition controls for set_new_abundances.

lgT_lo_for_set_new_abundances = 5.2d0
lgT_hi_for_set_new_abundances = 5.5d0

mesh adjustment#

max_allowed_nz#

Maximum number of grid points allowed. Array allowed to grow arbitrarily large if max_allowed_nz <= 0 If nz > max_allowed_nz, MESA will prematurely exit and call this error.

max_allowed_nz = 8000

remesh_max_allowed_logT#

Turn off remesh if any cell has logT > this.

remesh_max_allowed_logT = 1d99

mesh_max_allowed_ratio#

Must be >= 2.5. Max ratio for mass of adjacent cells. If have ratio exceeding this, split the larger cell.

mesh_max_allowed_ratio = 2.5d0

max_delta_x_for_merge#

Don’t merge neighboring cells if any abundance differs by more than this.

max_delta_x_for_merge = 0.1d0

mesh_delta_coeff#

A larger value increases the max allowed deltas and decreases the number of grid points. and a smaller does the opposite.

analogous to time_delta_coeff for better time resolution.

E.g., you’ll roughly double the number of grid points if you cut mesh_delta_coeff in half. Don’t expect it to exactly double the number however since other parameters in addition to gradients also influence the details of the grid spacing.

this factor also applies to max_dq, max_center_cell_dq, and max_surface_cell_dq.

mesh_delta_coeff = 1.0d0

mesh_Pgas_div_P_exponent#

Multiply mesh_delta_coeff by (Pgas/Ptotal) to this power.

mesh_Pgas_div_P_exponent = 0

mesh_delta_coeff_for_highT#

Use different mesh_delta_coeff at higher temperatures.

mesh_delta_coeff_for_highT = 3.0d0

logT_max_for_standard_mesh_delta_coeff#

Use mesh_delta_coeff for center logT <= this. This value should be less than logT_min_for_highT_mesh_delta_coeff.

logT_max_for_standard_mesh_delta_coeff = 9.0d0

logT_min_for_highT_mesh_delta_coeff#

Use mesh_delta_coeff_for_highT for center logT >= this. Linearly interpolate in logT for intermediate center temperatures.

logT_min_for_highT_mesh_delta_coeff = 9.5d0

max_dq#

Max size for cell as fraction of total mass.

max_dq = 1d-2

min_dq#

Min size for cell as fraction of total mass.

min_dq = 1d-14

Min size for cell to be split.

min_dq_for_split = 1d-14

min_dq_for_xa#

Min size for splitting because of composition gradient. only for non-convective regions if have set min_dq_for_xa_convective > 0.

min_dq_for_xa = 1d-14

min_dq_for_xa_convective#

Min size for splitting because of composition gradient in convective regions. if <= 0, then use min_dq_for_xa instead of this.

min_dq_for_xa_convective = 1d-6

Min size for cell to be split because of jump in logT.

min_dq_for_logT = 1d-14

mesh_min_dlnR#

Limit on difference in lnR across cell for mesh refinement. Do not make this smaller than about 1d-14 or will fail with numerical problems.

mesh_min_dlnR = 1d-9

merge_if_dlnR_too_small#

If true, mesh adjustment will force merge if difference in lnR across cell is too small.

merge_if_dlnR_too_small = .false.

mesh_min_dr_div_dRstar#

Limit on relative radial extent for mesh refinement. dRstar = s% r(1) - s% R_center Don’t split if dr/dRstar would drop below this limit.

mesh_min_dr_div_dRstar = -1

merge_if_dr_div_dRstar_too_small#

If true, mesh adjustment will force merge if dr_div_dRstar too small.

merge_if_dr_div_dRstar_too_small = .true.

mesh_min_dr_div_cs#

Limit (in seconds) on sound crossing time for mesh refinement. Don’t split if sound crossing time would drop below this limit.

mesh_min_dr_div_cs = -1

merge_if_dr_div_cs_too_small#

If true, mesh adjustment will force merge if dr_div_cs too small.

merge_if_dr_div_cs_too_small = .true.

max_center_cell_dq#

Largest allowed dq at center.

max_center_cell_dq = 1d-7

max_surface_cell_dq#

Largest allowed dq at surface.

max_surface_cell_dq = 1d-12

max_num_subcells#

Limits number of new cells from 1 old one.

max_num_subcells = 2

max_num_merge_cells#

Limits number of old cells to merge into 1 new one.

max_num_merge_cells = 2

mesh_adjust_get_T_from_E#

If true, then use internal energy conservation to set new temperature. If false, just use average temperature based on reconstruction polynomials.

mesh_adjust_get_T_from_E = .true.

mesh_ok_to_merge#

mesh_max_k_old_for_split#

mesh_min_k_old_for_split#

mesh_ok_to_merge = .true.
mesh_max_k_old_for_split = 999999999
mesh_min_k_old_for_split = 0

E_function_weight#

internal energy gradient, E_function = E_function_weight*max(E_function_param,log10(energy)).

E_function_weight = 0
E_function_param = 16d0

P_function_weight#

Pressure gradient, P_function = P_function_weight*log10(P).

P_function_weight = 40

T_function1_weight#

Temperature gradient, T_function1 = T_function1_weight*log10(T). NOTE: The T gradient mesh controls below seems to be necessary to allow burning that starts off center to be able to reach the center. You can see this in the pre_zahb test_suite case if you try running it without the T function. The center temperature will fail to rise.

T_function1_weight = 110

T_function2_weight#

T_function2_param#

T_function2 = T_function2_weight*log10(T / (T + T_function2_param))

Largest change in T_function2 happens around T = T_function2_param. Default value puts this in the envelope ionization region.

T_function2_weight = 0
T_function2_param = 2d4

R_function_weight#

R_function_param#

log radius gradient

R_function = R_function_weight*log10(1 + (r/Rsun)/R_function_param)
R_function_weight = 0
R_function_param = 1d-4

R_function2_weight#

R_function2_param1#

R_function2_param2#

R_function2 = R_function2_weight*min(R_function2_param1,max(R_function2_param2,r/Rstar))

where Rstar = radius of outer edge of model.

R_function2_weight = 0
R_function2_param1 = 0.4d0
R_function2_param2 = 0

R_function3_weight#

radius gradient

R_function3 = R_function3_weight*(r/Rstar)
R_function3_weight = 0

M_function_weight#

M_function_param#

log mass gradient

M_function = M_function_weight*log10(1 + (m/Msun)/M_function_param)
M_function_weight = 0
M_function_param = 1d-6

gradT_function_weight#

gradT gradient, gradT_function = gradT_function_weight*gradT

gradT_function_weight = 0

log_tau_function_weight#

log_tau gradient (optical depth)

log_tau_function = log_tau_function_weight*log10(tau)
log_tau_function_weight = 0

log_kap_function_weight#

log_kap gradient (optical depth)

log_kap_function = log_kap_function_weight*log10(kap)
log_kap_function_weight = 0

omega_function_weight#

omega gradient (rotation omega in rad/sec)

omega_function = omega_function_weight*log10(omega)
omega_function_weight = 0

gam_function_weight#

gam_function_param1#

gam_function_param2#

For extra resolution around liquid/solid transition.

gam = plasma interaction parameter
gam_function = gam_function_weight*tanh((gam - gam_function_param1)/gam_function_param2)
gam_function_weight = 0
gam_function_param1 = 170
gam_function_param2 = 20

xa_function_species#

xa_function_weight#

Mass fraction gradients.

xa_function = xa_function_weight*log10(xa + xa_function_param),

Up to num_xa_function of these - see star_def for value of num_xa_function. 0 length string means skip, otherwise name of nuclide as defined in chem_def. weight <= 0 means skip.

xa_function_species(:) = ''
xa_function_weight(:) = 0
xa_function_species(1) = 'he4'
xa_function_weight(1) = 30
xa_function_param(1) = 1d-2

xa_mesh_delta_coeff#

Useful if you want to increase mesh_delta_coeff during advanced burning. If xa_function_species(j) has the largest atomic number in current set of species, then multiply mesh_delta_coeff by xa_mesh_delta_coeff(j).

xa_mesh_delta_coeff(:) = 1

mesh_delta_coeff_factor_smooth_iters#

Some smoothing is useful when using local changes to mesh_delta_coeff.

mesh_delta_coeff_factor_smooth_iters = 3

“Indirect” mesh controls work by increasing sensitivity in selected regions. They work in the same way as mesh_delta_coeff – values less than 1.0 mean smaller allowed jumps in mesh functions and hence smaller grid points and higher resolution. But whereas mesh_delta_coeff applies uniformly to all cells, the “extra” coefficients can vary in value from one cell to the next.

Note that you can set your own local changes by means of the hook other_mesh_delta_coeff_factor.

mesh_logX_species#

mesh_logX_min_for_extra#

Increase resolution at points with large abs(dlogX/dlogP); logX = log10(X mass fraction).

mesh_logX_species(1) = ''
mesh_logX_min_for_extra(1) = -6

mesh_dlogX_dlogP_extra(1)#

mesh_dlogX_dlogP_full_on(1)#

mesh_dlogX_dlogP_full_off(1)#

Only increase resolution if logX >= mesh_logX_min_for_extra. Make mesh_dlogX_dlogP_extra < 1 for smaller allowed change in logP and hence higher resolution. Full effect if abs(dlogX/dlogP) >= mesh_dlogX_dlogP_full_on. No effect if abs(dlogX/dlogP)) <= mesh_dlogX_dlogP_full_off. Up to num_mesh_logX of these (see star_def for value of num_mesh_logX).

mesh_dlogX_dlogP_extra(1) = 1
mesh_dlogX_dlogP_full_on(1) = 2
mesh_dlogX_dlogP_full_off(1) = 1

Multiply mesh_delta_coeff near convection zone boundary (czb) by the following factors. Value < 1 gives increased resolution.

Increase resolution at points with large abs(dlog_eps/dlogP) for nuclear power eps (ergs/g/sec). At any particular location, only use eps nuc category with max local value e.g., only use mesh_dlog_pp_dlogP_extra at points where pp is the max burn source.

mesh_dlog_eps_min_for_extra#

Only increase resolution if log_eps >= mesh_dlog_eps_min_for_extra.

mesh_dlog_eps_min_for_extra = -2

mesh_dlog_eps_dlogP_full_on#

Full effect if abs(dlog_eps/dlogP) >= mesh_dlog_eps_dlogP_full_on.

mesh_dlog_eps_dlogP_full_on = 4

mesh_dlog_eps_dlogP_full_off#

No effect if abs(dlog_eps/dlogP)) <= mesh_dlog_eps_dlogP_full_off.

mesh_dlog_eps_dlogP_full_off = 1

Multiply the allowed change between adjacent cells by the following factors; (small factor => smaller allowed change => more cells).

pp and cno burning

mesh_dlog_pp_dlogP_extra = 0.25d0
mesh_dlog_cno_dlogP_extra = 0.25d0

triple alpha, c, n, and o burning

mesh_dlog_3alf_dlogP_extra = 0.25d0
mesh_dlog_burn_c_dlogP_extra = 0.25d0
mesh_dlog_burn_n_dlogP_extra = 0.25d0
mesh_dlog_burn_o_dlogP_extra = 0.25d0

ne, na, and mg burning

mesh_dlog_burn_ne_dlogP_extra = 0.25d0
mesh_dlog_burn_na_dlogP_extra = 0.25d0
mesh_dlog_burn_mg_dlogP_extra = 0.25d0

c12+c12. c12+o16, and o16+o16 burning

mesh_dlog_cc_dlogP_extra = 0.25d0
mesh_dlog_co_dlogP_extra = 0.25d0
mesh_dlog_oo_dlogP_extra = 0.25d0

si to iron along alpha chain burning

mesh_dlog_burn_si_dlogP_extra = 0.25d0
mesh_dlog_burn_s_dlogP_extra = 0.25d0
mesh_dlog_burn_ar_dlogP_extra = 0.25d0
mesh_dlog_burn_ca_dlogP_extra = 0.25d0
mesh_dlog_burn_ti_dlogP_extra = 0.25d0
mesh_dlog_burn_cr_dlogP_extra = 0.25d0
mesh_dlog_burn_fe_dlogP_extra = 0.25d0

photodisintegration burning

mesh_dlog_pnhe4_dlogP_extra = 0.25d0
mesh_dlog_other_dlogP_extra = 0.25d0
mesh_dlog_photo_dlogP_extra = 1

convective_bdy_weight#

convective_bdy_dq_limit#

convective_bdy_min_dt_yrs#

Mesh function to enhance resolution near convective boundaries

convective_bdy_weight = 0
convective_bdy_dq_limit = 3d-5
convective_bdy_min_dt_yrs = 1d-3

max_rel_delta_IE_for_mesh_total_energy_balance#

remeshing can adjust internal energy of cell by this fraction in order to maintain total internal + potential + kinetic energy.

max_rel_delta_IE_for_mesh_total_energy_balance = 0.05d0

trace_mesh_adjust_error_in_conservation#

If true, report relative errors for total PE, KE, and IE. (potential, kinetic, internal).

trace_mesh_adjust_error_in_conservation = .false.

okay_to_remesh#

If false, then no remeshing.

okay_to_remesh = .true.

restore_mesh_on_retry#

If true, then after a retry the remeshing is undone for the step, and the following try is performed with the same mesh used in the previous step. This can help with the retry by reducing the changes.

restore_mesh_on_retry = .false.

num_steps_to_hold_mesh_after_retry#

When restore_mesh_on_retry=true, then after a retry remeshing is not done for this number of steps.

num_steps_to_hold_mesh_after_retry = 0

remesh_dt_limit#

No remesh if dt < remesh_dt_limit, in seconds.

remesh_dt_limit = -1

use_split_merge_amr#

use_split_merge_amr = .false.

split_merge_amr_logtau_zoning#

split_merge_amr_log_zoning#

split_merge_amr_hybrid_zoning#

split_merge_amr_flipped_hybrid_zoning#

if split_merge_amr_logtau_zoning, target is even grid spacing in logtau else if split_merge_amr_log_zoning, target is even grid spacing in logr else if split_merge_amr_hybrid_zoning, target is even r spacing for core, even logr outside else if split_merge_amr_flipped_hybrid_zoning, target is even logr spacing for core, even r outside else target is even grid spacing in r

split_merge_amr_logtau_zoning = .false.
split_merge_amr_log_zoning = .true.
split_merge_amr_hybrid_zoning = .false.
split_merge_amr_flipped_hybrid_zoning = .false.

split_merge_amr_nz_baseline#

split_merge_amr_nz_r_core#

split_merge_amr_nz_r_core_fraction#

split_merge_amr_mesh_delta_coeff#

split_merge_amr_nz_baseline = 1000
split_merge_amr_nz_r_core = 0d0 ! ignore if <= 0
split_merge_amr_nz_r_core_fraction = 0d0 ! ignore if <= 0; else r_core = r_center + fraction*(r(1) - r_center)
split_merge_amr_mesh_delta_coeff = 1d0 ! like mesh_delta_coeff, but for amr

split_merge_amr_MaxLong#

split_merge_amr_MaxShort#

split cell if ratio of actual/desired size is > split_merge_amr_MaxLong; ignore if <= 0 merge cell if ratio of desired/actual size is < split_merge_amr_MaxShort; ignore if <= 0 be careful to avoid inconsistent limits such as when a required split triggers a required merge. to be safe, make sure product of limits > 2.

split_merge_amr_MaxLong = 1.5d0
split_merge_amr_MaxShort = 1.5d0

merge_amr_max_abs_du_div_cs#

merge_amr_max_abs_du_div_cs = 0.1d0

merge_amr_du_div_cs_limit_only_for_compression#

If true, then merge_amr_max_abs_du_div_cs limit will only be considered for cells that would undergo compression

merge_amr_du_div_cs_limit_only_for_compression = .false.

merge_amr_inhibit_at_jumps#

merge_amr_inhibit_at_jumps = .false.

merge_amr_ignore_surface_cells#

merge_amr_k_for_ignore_surface_cells#

Merging surface cells can cause problems. If merge_amr_ignore_surface_cells is true, then the outermost merge_amr_k_for_ignore_surface_cells cells are ignored for merge.

merge_amr_ignore_surface_cells = .true.
merge_amr_k_for_ignore_surface_cells = 2

split_merge_amr_avoid_repeated_remesh#

If true, then after a cell has been merged or split, the resulting cell will not be considered in further remeshing for this step.

split_merge_amr_avoid_repeated_remesh = .false.

split_merge_amr_dq_min#

split_merge_amr_dq_max#

split_merge_amr_dq_min = 1d-14
split_merge_amr_dq_max = 1d0

split_merge_amr_r_core_cm#

split_merge_amr_r_core_cm = 1d8

split_merge_amr_max_iters#

split_merge_amr_max_iters = 100

split_merge_amr_okay_to_split_1#

split_merge_amr_okay_to_split_nz#

split_merge_amr_okay_to_split_1 = .true.
split_merge_amr_okay_to_split_nz = .true.

equal_split_density_amr#

equal_split_density_amr = .false.

trace_split_merge_amr#

trace_split_merge_amr = .false.

nuclear reaction controls#

default_net_name#

Name of base reaction network. Each net corresponds to a file in $MESA_DIR/data/net_data/nets. Look in that directory to see your network options, or learn how to create your own net.

default_net_name = 'basic.net'

screening_mode#

  • empty string means no screening

  • ' extended' : extends the Graboske (1973) method using results from Alastuey and Jancovici (1978), along with plasma parameters from Itoh et al (1979) for strong screening.

  • 'salpeter' : weak screening only. following Salpeter (1954), with equations (4-215) and (4-221) of Clayton (1968).

  • 'chugunov' : based on code from Sam Jones Implements screening from Chugunov et al (2007) for weak and strong screening MESA versions <=11435 used extended as the default value

screening_mode = 'chugunov'

net_logTcut_lo#

strong rates are zero logT < logTcut_lo use default from net if this is <= 0

net_logTcut_lo = -1

net_logTcut_lim#

strong rates cutoff smoothly for logT < logTcut_lim use default from net if this is <= 0

net_logTcut_lim = -1

max_abar_for_burning#

if abar > this, suppress all burning e.g., if want an “inert” core heavy elements, set this to 55 or, if want to turn off the net, set this to -1

max_abar_for_burning = 199

dxdt_nuc_factor#

Control for abundance changes by burning. Changes dxdt_nuc (rate of change of abundances) without changing the rates or eps_nuc (rate of energy generation).

dxdt_nuc_factor = 1

weak_rate_factor#

all weak rates are multiplied by this factor

weak_rate_factor = 1

reaction_neuQs_factor#

all neutrino Q factors are multiplied by this factor

reaction_neuQs_factor = 1

nonlocal_NiCo_kap_gamma#

nonlocal_NiCo_kap_gamma = 0

nonlocal_NiCo_decay_heat#

if true, do non-local deposition of gamma-ray energy from Ni56 and Co56 decays. only for approx nets including co56. intended for use with stripped envelope supernovae.

nonlocal_NiCo_decay_heat = .false.

dtau_gamma_NiCo_decay_heat#

dtau_gamma_NiCo_decay_heat = 1d0

max_logT_for_net#

max_logT_for_net = 10.2d0

element diffusion#

gravitational settling and chemical diffusion.

show_diffusion_info#

terminal output for diffusion

show_diffusion_info = .false.

show_diffusion_substep_info#

terminal output for diffusion

show_diffusion_substep_info = .false.

show_diffusion_timing#

show time for each call on diffusion

show_diffusion_timing = .false.

do_element_diffusion#

determines whether or not we do element diffusion

do_element_diffusion = .false.

diffusion_dt_limit#

no element diffusion if dt < this limit (in seconds)

diffusion_dt_limit = 3.15d7

diffusion_use_paquette#

if true, use atomic diffusion coefficients according to Paquette et al. (1986). if false, use Stanton & Murillo PhRvE, 93, 043203 (2016) for diffusion coefficients. (Paquette coefficients still used for electron-ion because Stanton & Murillo did not do calculations for attractive potentials.)

diffusion_use_paquette = .false.

diffusion_use_caplan#

if true, use atomic diffusion coefficients according to Caplan, Bauer, & Freeman MNRAS, 513, L52 (2022) at strong coupling (Gamma > 10), relevant for white dwarf interiors.

diffusion_use_caplan = .true.

diffusion_use_iben_macdonald#

if true, use diffusion coefficients similar to Iben & MacDonald (1985). if false, use Stanton & Murillo (2016) for diffusion coefficients. this was previously called diffusion_use_pure_coulomb.

diffusion_use_iben_macdonald = .false.

diffusion_use_cgs_solver#

if false, solve the system of equations described by Thoul et al. (1994) if true, solve the unmodified Burgers equations in cgs units

diffusion_use_cgs_solver = .true.

cgs_thermal_diffusion_eta_full_on#

cgs_thermal_diffusion_eta_full_off#

When diffusion_use_cgs_solver = .true. for eta < cgs_thermal_diffusion_eta_full_on, includes the heat flow vector terms in the Burgers equations. Then smoothly turns off use of these terms so that they are not included for eta > cgs_thermal_diffusion_eta_full_off, since these terms are problematic when distribution function become non-Maxwellian.

cgs_thermal_diffusion_eta_full_on = 0d0
cgs_thermal_diffusion_eta_full_off = 2d0

do_WD_sedimentation_heating#

min_xa_for_WD_sedimentation_heating#

if true, include heating associated with sedimentation when element diffusion is on. Only elements with mass fraction > min_xa_for_WD_sedimentation_heating will be included in this calculation. For best results, set diffusion_use_full_net = .true. This will affect white dwarf cooling times. See also eps_WD_sedimentation_factor

do_WD_sedimentation_heating = .false.
min_xa_for_WD_sedimentation_heating = 1d-5

do_diffusion_heating#

if true, calculate heating term associated with changes in internal energy due to any abundance changes from element diffusion, and include this term in the energy equation. To avoid double-counting, this control can only be used if do_WD_sedimentation_heating = .false.

do_diffusion_heating = .true.

diffusion_min_dq_at_surface#

treat at least this much at surface as a single cell for purposes of diffusion

diffusion_min_dq_at_surface = 1d-9

diffusion_min_T_at_surface#

treat cells cells at surface with T < this as a single cell for purposes of diffusion default should be large enough to ensure hydrogen ionization

diffusion_min_T_at_surface = 1d4

diffusion_min_dq_ratio_at_surface#

combine cells at surface until have total mass >= this factor times the next cell below them this helps with surface boundary condition for diffusion by putting large cell at surface

diffusion_min_dq_ratio_at_surface = 10

diffusion_dt_div_timescale#

dt is at most this fraction of timescale. Each stellar evolution step can be divided into many substeps for diffusion. The substep timescale is set by rates of flow in and out for each species in each cell. The substep size, dt, is initially set to timescale*diffusion_dt_div_timescale.

diffusion_dt_div_timescale = 1

diffusion_min_num_substeps#

Max substep dt is total time divided by this.

diffusion_min_num_substeps = 1

diffusion_max_iters_per_substep#

If the substep requires too many iterations, the substep time is decreased for a retry.

diffusion_max_iters_per_substep = 10

diffusion_max_retries_per_substep#

If the substep requires too many retries, diffusion fails and forces a retry for the star.

diffusion_max_retries_per_substep = 10

diffusion_tol_correction_max#

diffusion_tol_correction_norm#

Tolerances for solver iterations. Corrections smaller will be treated as converged. Corrections larger will cause another solver iteration.

diffusion_tol_correction_max = 1d-1
diffusion_tol_correction_norm = 1d-3

diffusion_min_X_hard_limit#

tolerance for negative mass fraction errors errors larger will cause retry; errors smaller will be corrected. Tightening this control may help with “failed in fixup” errors when diffusion_use_isolve = .true.

diffusion_min_X_hard_limit = -1d-3

diffusion_X_total_atol#

diffusion_X_total_rtol#

tolerances for errors in total species conservation errors larger will cause retry; errors smaller will be corrected.

diffusion_X_total_atol = 1d-9
diffusion_X_total_rtol = 1d-6

diffusion_upwind_abs_v_limit#

switch to upwind for i at face k if abs(v(i,k)) > this limit mainly for use with radiative levitation where get very much higher velocities

diffusion_upwind_abs_v_limit = 1d99

diffusion_v_max#

Max velocity (cm/sec). We can get extremely large velocities in the extreme outer envelope that cause problems numerically without really effecting the results, so we allow a max for the velocities that should help the numerics without changing the results. Note: change diffusion_v_max to at least 1d-2 when using radiative levitation.

diffusion_v_max = 1d-3

D_mix_ignore_diffusion#

Diffusion is turned off in core and surface convection zones, since it is overwhelmed by other mixing there. D_mix_ignore_diffusion roughly defines the mixing coefficient below which diffusion is included again. The code finds the location where D_mix falls to this value, backs up some, and turns on diffusion from there onward.

D_mix_ignore_diffusion = 1d5

diffusion_gamma_full_off#

diffusion_gamma_full_on#

gamma_full_on <= gamma_full_off Shut off diffusion for large gamma (i.e. for gamma >= gamma_full_off). Gradually decrease diffusion as gamma increases from full_on to full_off. Allow normal diffusion for gamma <= gamma_full_on. Default is diffusion off when get well into liquid regime.

diffusion_gamma_full_off = 1d99
diffusion_gamma_full_on = 1d99

diffusion_T_full_on#

diffusion_T_full_off#

T_full_on >= T_full_off Shut off diffusion for small T (i.e., for T <= T_full_off) Gradually decrease diffusion as T decreases from T_full_on to T_full_off. Allow normal diffusion for T >= T_full_on.

diffusion_T_full_on = 1d3
diffusion_T_full_off = 1d3

diffusion_calculates_ionization#

If diffusion_calculates_ionization is false, MESA uses typical charges for a set of representative species as defined in diffusion_class_typical_charge and diffusion_class_representative for all points rather than calculating the ionization from the local conditions.

diffusion_calculates_ionization = .true.

diffusion_nsmooth_typical_charge#

smoothing over charge

diffusion_nsmooth_typical_charge = 10

diffusion_SIG_factor#

diffusion_GT_factor#

factors for playing with SIG and GT terms for concentration diffusion and advection

diffusion_SIG_factor = 1d0
diffusion_GT_factor = 1d0

diffusion_AD_dm_full_on#

diffusion_AD_dm_full_off#

diffusion_AD_boost_factor#

artificial concentration diffusion near surface (mainly for radiative levitation) Msun units for full_on and full_off boost only used if > 0

diffusion_AD_dm_full_on = -1
diffusion_AD_dm_full_off = -1
diffusion_AD_boost_factor = 0

diffusion_Vlimit_dm_full_on#

diffusion_Vlimit_dm_full_off#

in Msun units artificial velocity limitation near surface (mainly for radiative levitation)

diffusion_Vlimit_dm_full_on = -1
diffusion_Vlimit_dm_full_off = -1

diffusion_Vlimit#

In units of local cell crossing velocity (only used if > 0). When full on, limit abs(v) <= Vlimit*dr/dt, cell size dr, substep time dt.

diffusion_Vlimit = 0

D_mix_zero_region_bottom_q#

D_mix_zero_region_top_q#

dq_D_mix_zero_at_H_He_crossover#

dq_D_mix_zero_at_H_C_crossover#

D_mix_zero_region_bottom_q = 1d99
D_mix_zero_region_top_q = -1d99
dq_D_mix_zero_at_H_He_crossover = -1d0
dq_D_mix_zero_at_H_C_crossover = -1d0

diffusion_min_T_for_radaccel#

diffusion_max_T_for_radaccel#

If T between these limits, then include radiative levitation at that location. Calculation of radiative levitation is costly, so only use it where necessary. Note: change diffusion_v_max to at least 1d-2 when using radiative levitation.

Note that radiative levitation requires OP calculations of g_rad for each class, and only 17 elements are supported (H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Cr, Mn, Fe, Ni). If you want to include radiative levitation, your options are: + Define diffusion classes such that all class representatives are among the 17 elements listed above. + Use a net with only elements from the 17 above, and set diffusion_use_full_net = .true.

diffusion_min_T_for_radaccel = 0
diffusion_max_T_for_radaccel = 0

diffusion_min_Z_for_radaccel#

diffusion_max_Z_for_radaccel#

If Z between these limits, then include radiative levitation for that element. Calculation of radiative levitation is costly, so only use it where necessary. e.g., limit to Fe and Ni by min_Z = 26 and max_Z = 28

diffusion_min_Z_for_radaccel = 0
diffusion_max_Z_for_radaccel = 1000

diffusion_screening_for_radaccel#

Include screening for radiative levitation.

diffusion_screening_for_radaccel = .true.

diffusion_use_full_net#

If true, don’t lump elements into classes for diffusion. Instead, each isotope in the network is treated as its own separate class. This can cause significant slowdowns for large nets, so it is off by default. This works for nets with up to 100 isotopes; larger nets require lumping into classes.

diffusion_use_full_net = .false.

diffusion_num_classes#

Number of representative classes of species for diffusion calculations. (maximum of 100)

diffusion_num_classes = 5

diffusion_class_representative(:)#

isotope names for diffusion representatives

diffusion_class_representative(1) = 'h1'
diffusion_class_representative(2) = 'he3'
diffusion_class_representative(3) = 'he4'
diffusion_class_representative(4) = 'o16'
diffusion_class_representative(5) = 'fe56'

diffusion_class_A_max(:)#

atomic number A. in ascending order. species goes into 1st class with A_max >= species A

diffusion_class_A_max(1) = 2
diffusion_class_A_max(2) = 3
diffusion_class_A_max(3) = 4
diffusion_class_A_max(4) = 16
diffusion_class_A_max(5) = 10000

diffusion_class_typical_charge(:)#

Typical charges for use if diffusion_calculates_ionization is false Use charge 21 for Fe in the sun, from Thoul, Bahcall, and Loeb (1994), ApJ, 421, 828.

diffusion_class_typical_charge(1) = 1
diffusion_class_typical_charge(2) = 2
diffusion_class_typical_charge(3) = 2
diffusion_class_typical_charge(4) = 8
diffusion_class_typical_charge(5) = 21

diffusion_class_factor(:)#

Arbitrarily enhance or inhibit diffusion effects by class.

diffusion_class_factor(:) = 1d0

parameters for ionization solver#

diffusion_use_isolve#

Activate iterative solver.

diffusion_use_isolve = .false.

diffusion_rtol_for_isolve#

diffusion_atol_for_isolve#

Relative and absolute error parameters for iterative solver.

diffusion_rtol_for_isolve = 1d-4
diffusion_atol_for_isolve = 1d-5

diffusion_maxsteps_for_isolve#

Maximum number of steps to take in iterative solver.

diffusion_maxsteps_for_isolve = 1000

diffusion_isolve_solver#

Which ode solver to use for iterative.

Options include:

  • 'ros2_solver'

  • 'rose2_solver'

  • 'ros3p_solver'

  • 'ros3pl_solver'

  • 'rodas3_solver'

  • 'rodas4_solver'

  • 'rodasp_solver'

diffusion_isolve_solver = 'ros2_solver'

diffusion_dump_call_number#

debugging info of diffusion at call number

diffusion_dump_call_number = -1

WD phase separation#

do_phase_separation#

Phase separation upon crystallization in WD cores using the implementation of Bauer (2023).

do_phase_separation = .false.

phase_separation_option#

Choice of appropriate option for the WD core mixture:

phase_separation_option = 'CO'

do_phase_separation_heating#

if true, calculate heating term associated with changes in internal energy due to any abundance changes from phase separation, and include this term in the energy equation.

do_phase_separation_heating = .true.

phase_separation_mixing_use_brunt#

if true, the phase separation mixing recalculates relevant EOS quantities and evaluates the Ledoux criterion, including the brunt B term that depends on the composition gradient. This can be somewhat expensive, so this option can be set to false to instead just mix outward until there is no more negative mu gradient. These will produce similar final chemical profiles, but setting this option to true is the only way to properly evaluate the physical criterion.

phase_separation_mixing_use_brunt = .true.

eos controls#

fix_d_eos_dxa_partials#

The star solver uses the partial derivatives of lnPgas and lnE with respect to composition. When the EOS fails to provide these, replace them with a finite-difference approximation.

fix_d_eos_dxa_partials = .true.

opacity controls#

more opacity controls can be found in star_job.defaults

use_simple_es_for_kap#

for experiments with simple electron scattering if true, opacity = 0.2*(1 + X)

use_simple_es_for_kap = .false.

use_starting_composition_for_kap#

for experiments with partials of opacity with respect to composition if true, calls on kap during solver iterations use the starting composition

use_starting_composition_for_kap = .false.

opacity_min#

limit minimum opacities to this value (ignore this if value is < 0)

opacity_min = -1

opacity_max#

limit opacities to this value (ignore this if value is < 0)

opacity_max = -1

opacity_factor#

opacities are multiplied by this value

opacity_factor = 1

min_logT_for_opacity_factor_off#

min_logT_for_opacity_factor_on and#

max_logT_for_opacity_factor_on#

max_logT_for_opacity_factor_off#

temperature controls for where the opacity_factor is applied if, for example, you only want the opacity factor to apply in the iron bump region you can give a logT range such as

min_logT_for_opacity_factor_off = 5.2
min_logT_for_opacity_factor_on = 5.3
max_logT_for_opacity_factor_on = 5.7
max_logT_for_opacity_factor_off = 5.8

ignore these if < 0.

min_logT_for_opacity_factor_off = -1
min_logT_for_opacity_factor_on = -1
max_logT_for_opacity_factor_on = -1
max_logT_for_opacity_factor_off = -1

if you need cell-by-cell control of opacity factor, set the vector “extra_opacity_factor” using the routine “other_opacity_factor

OP mono opacities#

The OP_mono opacities use data and code from the OP website as modified by Haili Hu. Since the tar.xz file is large (462 MB), it is not included in the standard mesa download.

You can get OP4STARS_1.3.tar.xz from https://zenodo.org/records/4390522

Put it any place you want on your disk.

tar -xf OP4STARS_1.3.tar.xz

Set the inlist controls for the “mono” directory with the data files. For example, in my case it looks like the following, but you can put the directory anywhere you like – it doesn’t need to be in the mesa/data directory. And the cache file doesn’t need to be in the mono directory.

op_mono_data_path = '/Users/bpaxton/OP4STARS_1.3/mono'
op_mono_data_cache_filename = '/Users/bpaxton/OP4STARS_1.3/mono/op_mono_cache.bin'

If you use these opacities, you should cite Seaton (2005).

op_mono_data_path#

if this path is set to the empty string, ‘’, then it defaults to the environment variable $(MESA_OP_MONO_DATA_PATH)

op_mono_data_path = ''      ! '' defaults to $MESA_OP_MONO_DATA_PATH

op_mono_data_cache_filename#

if this is set to the empty string, ‘’, then it defaults to the environment variable $(MESA_OP_MONO_DATA_CACHE_FILENAME)

op_mono_data_cache_filename = ''      ! '' defaults to $MESA_OP_MONO_DATA_CACHE_FILENAME

op_mono_method#

Compute the Rosseland mean opacity and radiative accelerations from the OP mono data by brute force ('hu') or use the faster method by allowing for a small tolerance on the mixture used for the computations of these variables ('mombarg').

op_mono_method = 'hu'

emesh_data_for_op_mono_path#

path to the OP_mono_master_grid_MESA_emesh.txt file containing the data need for when op_mono_method = 'mombarg'. If this is set to the empty string, ‘’, then it defaults to the environment variable $(MESA_OP_MONO_MASTER_GRID)

You can get OP_mono_master_grid_MESA_emesh.txt from https://doi.org/10.5281/zenodo.6850861

You can either download the uncompressed .txt file, or download the compressed .xz file and then unpack it in place with unxz -v OP_mono_master_grid_MESA_emesh.txt.xz

emesh_data_for_op_mono_path = ''      ! '' defaults to $MESA_OP_MONO_MASTER_GRID

high_logT_op_mono_full_off#

high_logT_op_mono_full_on#

low_logT_op_mono_full_off#

low_logT_op_mono_full_on#

Blending controls for turning op_mono opacities off above (high_logT) and below (low_logT) specified temperature ranges. When not using op_mono, the code will use standard opacity tables. For example, you might only use high T limits so that op_mono is only used in the envelope, or you might set both low and high T limits so that op_mono is used around the Fe peak logT but not for other locations in the star. These controls should satisfy the following inequalities:

high_logT_op_mono_full_off >= high_logT_op_mono_full_on
high_logT_op_mono_full_on >= low_logT_op_mono_full_on
low_logT_op_mono_full_on >= low_logT_op_mono_full_off
op_mono opacities fully on if

log10T <= high_logT_op_mono_full_on and log10T >= low_logT_op_mono_full_on

op_mono opacities full off if

log10T >= high_logT_op_mono_full_off or log10T <= low_logT_op_mono_full_off

Note: OP mono ignored if either high_logT control < 0

high_logT_op_mono_full_off = -99
high_logT_op_mono_full_on = -99
low_logT_op_mono_full_off = -99
low_logT_op_mono_full_on = -99

op_mono_min_X_to_include#

skip iso if mass fraction < this

op_mono_min_X_to_include = 1d-20

use_op_mono_alt_get_kap#

if true, call the op_mono_alt_get_kap routine instead of op_mono_get_kap. see mesa/kap/public/kap_lib.f for details about these routines.

use_op_mono_alt_get_kap = .false.

min_kap_for_dPrad_dm_eqn#

min_kap_for_dPrad_dm_eqn = 1d-4

asteroseismology controls#

get_delta_nu_from_scaled_solar#

If get_delta_nu_from_scaled_solar is .false., the large separation delta_nu is the inverse of the sound crossing time from one side of the star to the other, through the center. This is sometimes called the “asymptotic” large separation.

Otherwise, delta_nu is calculated from the asteroseismic scaling relations (see Ulrich 1986, Brown et al. 1991 and Kjeldsen & Bedding 1995) using solar reference values nu_max_sun, delta_nu_sun and astero_Teff_sun.

nu_max is always computed from the scaling relation.

get_delta_nu_from_scaled_solar = .false.

nu_max_sun#

delta_nu_sun#

astero_Teff_sun#

Solar reference values used in the asteroseismic scaling relations for delta_nu (if get_delta_nu_from_scaled_solar is .false.) and nu_max (always).

The default nu_max_sun is the Sun-as-as-star value reported by Lund et al. (2017), which is consistent with but conceptually different from the result of 3073.59 ± 0.18 μHz by Kiefer et al. (2018).

The default delta_nu_sun is also taken from Lund et al. (2017).

The default astero_Teff_sun is the value adopted in IAU 2015 Resolution B3. This should not be confused with the constant Teffsun, which is always equal to the IAU value and not controlled by a parameter. The “asteroseismic” value can be changed in case one needs to reproduce previous calculations using the scaling relations.

nu_max_sun = 3078d0 ! μHz
delta_nu_sun = 134.91d0 ! μHz
astero_Teff_sun = 5772d0 ! kelvin

delta_Pg_traditional#

delta_Pg_mode_freq#

For calculating the asymptotic g-mode period spacing delta_Pg

if delta_Pg_traditional = .true., then calculate the asymptotic gravity mode period space by directly integrating N2 over the entire stellar model.

if delta_Pg_traditional = .false. use the method of Bildsten et al. (2012): //ui.adsabs.harvard.edu/abs/2012ApJ…744L…6B/abstract

delta_Pg_mode_freq is only used when delta_Pg_traditional = .false. it is provided in units of uHz. if <=0, use nu_max from scaled solar value

delta_Pg_traditional = .true.
delta_Pg_mode_freq = 0d0

Brunt controls#

calculate_Brunt_B#

calculate_Brunt_N2#

Only calculate if this is true.

calculate_Brunt_B = .true.
calculate_Brunt_N2 = .true.

brunt_N2_coefficient#

Standard N2 is multiplied by this value.

brunt_N2_coefficient = 1

num_cells_for_smooth_brunt_B#

Number of cells on either side to use in weighted smoothing of brunt_B.

num_cells_for_smooth_brunt_B = 2

threshold_for_smooth_brunt_B#

Threshold for weighted smoothing of brunt_B. Only apply smoothing (controlled by num_cells_for_smooth_brunt_B) for contiguous regions where \(|B|\) exceeds this threshold. Might be useful for preventing narrow peaks from being excessively broadened by smoothing

threshold_for_smooth_brunt_B = 0d0

min_magnitude_brunt_B#

If set brunt_B to 0 if absolute value is < this.

min_magnitude_brunt_B = -1d99

structure equations#

energy_eqn_option#

Available options are 'dedt' or 'eps_grav'. See below for descriptions of each form and form-specific options.

energy_eqn_option = 'dedt'

dedt form#

This form of the energy equation is used when energy_eqn_option = 'dedt'.

It is a conservative equation for the local specific total energy introduced in MESA V, Section 3. See in particular Eq. (8) and surrounding discussion.

Because this equation is written in a conservative form, it should always do an excellent job of numerical energy conservation. The error in numerical energy conservation (quantified by rel_E_err) reflects the energy equation residuals (i.e., the extent to which the energy equation was not satisfied). When using this form of the equation, models should generally have a small (\(\lesssim 10^{-8}\)) value of rel_E_err, roughly independent of space and time resolution. A small value of rel_E_err and its cumulative counterpart rel_run_E_err only demonstrates that the equation residuals are small and is not evidence that a model is converged or reliable. Convergence studies targeting the physical quantities of interest remain essential. A large value of rel_run_E_err (\(\gtrsim 10^{-2}\)) should be a cause for concern and should be investigated further (see also warn_when_large_rel_run_E_err and max_abs_rel_run_E_err )

no_dedt_form_during_relax#

dedt_eqn_r_scale#

no_dedt_form_during_relax = .true.
dedt_eqn_r_scale = 1d8

eps_grav form#

This form of the energy equation is used when energy_eqn_option = 'eps_grav'.

The quantity eps_grav is defined as \(\epsilon_{\rm grav} = -\frac{De}{Dt} - P \frac{DV_\rho}{Dt}\), where \(e\) is the specific internal energy, \(P\) is the pressure, \(V_\rho \equiv 1/\rho\) is the specific volume, and \(D/Dt\) is the Lagrangian time derivative. See MESA IV, Section 8 for more discussion.

This quantity is then re-written into the following convenient-to-evaluate form (see MESA IV, eq. 63):

\(\epsilon_{\rm grav} = -T c_P \left[(1 - \nabla_{\rm ad} \chi_T)\frac{D\ln T}{Dt} - \nabla_{\rm ad} \chi_\rho \frac{D\ln \rho}{Dt}\right] + \epsilon_{\rm grav, X}\).

The final term reflects the change in internal energy due to changes in composition (at fixed density and temperature) and is referred to in MESA as eps_grav_composition_term. It is defined as

\(\epsilon_{\rm grav, X} \equiv -\sum_i \left(\frac{\partial e}{\partial X_i}\right)_{\rho,T, \{X\ne X_i\}} \frac{DX_i}{Dt}\),

and MESA evaluates this term using the finite difference

\(\epsilon_{\rm grav, X} = -\frac{1}{\delta t}\left[e(\rho, T, X) - e(\rho, T, X_{\rm start})\right]\),

where \(\delta t\) is the timestep and \(X_{\rm start}\) is the start-of-step mass fractions. (Other quantities take their end-of-step values.)

Note: In a phase transition, eps_grav includes the latent heat.

As with the dedt form, the error in numerical energy conservation (quantified by rel_E_err) reflects the energy equation residuals (i.e., the extent to which the energy equation was not satisfied). However, because this version of the energy equation is not written in a conservative form, it also includes error associated with the time discretization. An additional source of error enters because the equation of state provided by the eos module does not precisely satisfy the mathematical and thermodynamic identities that are used in rewriting the total and partial derivatives present in the equation. This inconsistency is usually worst at the conditions where MESA blends different component EOSes. It is important to understand that time discretization error and eos inconsistencies also affect models using the dedt form of the energy equation but manifest in different ways (e.g., as entropy generation). Under degenerate conditions, it is often preferable to incur energy errors rather than entropy errors, and the eps_grav form should generally be preferred in such circumstances.

In practice, the error sources usually exhibit the ordering (time discretization) > (eos inconsistency) >> (equation residuals). With increasing time resolution, the time discretization error can be driven down to the floor imposed by the eos inconsistency error. The value of this floor depends on the physical conditions, but may be \(\sim 10^{-4}\), well above the level of the residuals (\(\lesssim 10^{-8}\)). The control use_time_centered_eps_grav provides a time-centered implementation of eps_grav that can often reach this floor at larger timesteps. Convergence studies targeting the physical quantities of interest remain essential.

include_composition_in_eps_grav#

If true, evaluate eps_grav_composition_term and include this quantity in eps_grav.

When this flag is true, the composition derivatives of eps_grav are also included in the Jacobian.

If this flag is set to false, the eps_grav form will not conserve energy in situations with changing composition.

include_composition_in_eps_grav = .true.

use_time_centered_eps_grav#

If true, use a time-centered version of eps_grav and eps_grav_composition_term. (Disabled during relax.)

use_time_centered_eps_grav = .true.

Gamma_lnS_eps_grav_full_off#

Gamma_lnS_eps_grav_full_on#

Automatic switch to lnS form for eps_grav (\(\epsilon_{\rm grav} = -T\frac{Ds}{Dt}\)) in regions with high Gamma (plasma interaction parameter). These controls only apply when using the PC EOS. This is necessary to get the latent heat associated with the crystallization phase transition. The composition term \(-\sum_i (\partial e/\partial Y_i)_{s,\rho} dY_i\) is never included when the lnS form is used, independent of the setting of the control include_composition_in_eps_grav.

Gamma_lnS_eps_grav_full_on = 150d0
Gamma_lnS_eps_grav_full_off = 120d0

eps_grav_factor#

multiply eps_grav by this factor

eps_grav_factor = 1d0

velocity_q_upper_bound#

Local override for global v_flag. If local q > this bound, local v_flag is set false, else local v_flag is set to global v_flag. this lets you force v = 0 in outer envelope.

velocity_q_upper_bound = 1d99

velocity_tau_lower_bound#

Local override for global v_flag. If local tau < this bound, local v_flag is set false, else local v_flag is set to global v_flag. this lets you force v = 0 in outer envelope.

velocity_tau_lower_bound = -1d99

velocity_logT_lower_bound#

Local override for global v_flag. If local logT < this bound, local v_flag is set false, else local v_flag is set to global v_flag. this lets you force v = 0 in outer envelope.

velocity_logT_lower_bound = -1d99

max_dt_yrs_for_velocity_logT_lower_bound#

Only apply velocity_logT_lower_bound when timestep < this limit.

max_dt_yrs_for_velocity_logT_lower_bound = 1d99

use_gravity_rotation_correction#

With rotation, multiply gravity by fp_rot if this flag is true. See the 2nd MESA paper (2013), equation 22. previously called “use_dP_dm_rotation_correction”.

use_gravity_rotation_correction = .true.

non_nuc_neu_factor#

Multiplies power from non-nuclear reaction neutrinos. i.e., thermal neutrinos such as computed by mesa/neu.

non_nuc_neu_factor = 1

eps_nuc_factor#

Multiplies eps_nuc without changing rates or dxdt_nuc. Thus controls energy production without modifying the amount of change in abundances.

eps_nuc_factor = 1

eps_WD_sedimentation_factor#

This controls energy production from sedimentation of Ne22 (and possibly other neutron-rich elements in WD interiors).

eps_WD_sedimentation_factor = 1

max_abs_eps_nuc#

Limit magnitude of eps_nuc to this.

max_abs_eps_nuc = 1d99

fe56ec_fake_factor#

min_T_for_fe56ec_fake_factor#

Multiplier on ni56 electron capture rate to take isotopes in hardwired networks to more neutron rich isotopes.

fe56ec_fake_factor = 1d-7
min_T_for_fe56ec_fake_factor = 3d9

eps_mdot_factor#

multiply eps_mdot by this factor

eps_mdot_factor = 1d0

max_num_surf_revisions#

Max number of forced reconverges for changes in surf_lnS.

max_num_surf_revisions = 1

max_abs_rel_change_surf_lnS#

Force solver reconverge if surf_lnS changed more than this.

max_abs_rel_change_surf_lnS = 5d-4

extra_power_source#

erg/g/sec applied uniformly throughout the model This can be used to push a pre-ms model up the track to lower center temperatures. Can be used simultaneously with inject_extra_ergs_sec and inject_uniform_extra_heat

extra_power_source = 0

inject_uniform_extra_heat#

extra heat in erg g^-1 s^-1 Added to cells in range min_q_for_uniform_extra_heat to max. Can be used simultaneously with inject_extra_ergs_sec and extra_power_source.

inject_uniform_extra_heat = 0

min_q_for_uniform_extra_heat#

sets bottom of region for inject_uniform_extra_heat

min_q_for_uniform_extra_heat = 0

max_q_for_uniform_extra_heat#

sets top of region for inject_uniform_extra_heat

max_q_for_uniform_extra_heat = 1

inject_extra_ergs_sec#

added to mass equal to grams_for_inject_extra_core_ergs_sec can be used simultaneously with extra_power_source and inject_uniform_extra_heat

inject_extra_ergs_sec = 0

base_of_inject_extra_ergs_sec#

(units: Msun) sets bottom of region for inject_extra_ergs_sec note: actual base is at max of this and the center of the model

base_of_inject_extra_ergs_sec = 0

total_mass_for_inject_extra_ergs_sec#

(units: Msun) sets size of region for inject_extra_ergs_sec

total_mass_for_inject_extra_ergs_sec = 0

start_time_for_inject_extra_ergs_sec#

(units: sec) start time for injecting extra ergs/s

start_time_for_inject_extra_ergs_sec = -1d99

duration_for_inject_extra_ergs_sec#

(units: sec) length of time for injecting extra ergs/s set to negative value to keep injecting indefinitely or until reach target

duration_for_inject_extra_ergs_sec = -1

inject_until_reach_model_with_total_energy#

(units: ergs) target for model total energy usually want to set duration_for_inject_extra_ergs_sec = -1 for this option. continue injecting until total energy of model reaches min of inject_until_reach_model_with_total_energy, and initial total energy

inject_until_reach_model_with_total_energy = 1d99

steps_before_use_velocity_time_centering#

include_P_in_velocity_time_centering#

include_L_in_velocity_time_centering#

use_P_d_1_div_rho_form_of_work_when_time_centering_velocity#

P_theta_for_velocity_time_centering#

L_theta_for_velocity_time_centering#

for time weighting in energy and momentum equations to give intrinsic conservation of total energy (conservation -> perfect as residuals -> 0), and to minimize numerical damping. as discussed in the 3rd mesa instrument paper (2015). time centering applies to velocities, pressures, and luminosities. not for u_flag. steps_before_use_velocity_time_centering < 0 means no time centering. P_theta and L_theta = 0.5 for time centered, = 1.0 for fully implicit.

steps_before_use_velocity_time_centering = -1
include_P_in_velocity_time_centering = .false.
P_theta_for_velocity_time_centering = 0.5d0
include_L_in_velocity_time_centering = .false.
L_theta_for_velocity_time_centering = 0.5d0
use_P_d_1_div_rho_form_of_work_when_time_centering_velocity = .false.

use_dPrad_dm_form_of_T_gradient_eqn#

use_flux_limiting_with_dPrad_dm_form#

use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn#

These are for alternatives ways to determine the T gradient. The standard form of the equation is

dT/dm = dP/dm * T/P * grad_T, grad_T = dlnT/dlnP from MLT.

use hydrostatic value for dP/dm in this. this is because of limitations of MLT for calculating grad_T. (MLT assumes hydrostatic equilibrium) see comment in K&W chpt 9.1.

The alternatives forms are for dynamic situations where the use of hydrostatic dP/dm is inappropriate. In order of priority,

if (use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn)
equation is gradT(k)*(lnPeos(k-1) - lnPeos(k)) = lnT(k-1) - lnT(k)
recall that gradT(k) is expected dlnT/dlnPeos at face k from MLT.

if (use_dPrad_dm_form_of_T_gradient_eqn)
if (gradT < gradr) then
use L_rad = L*gradT/gradr (see, e.g., Cox&Giuli 14.109)
else
use L_rad = L
if (use_flux_limiting_with_dPrad_dm_form)
flxR = 4 * pi * r^2 * abs(dT^4/dm) / kap
flxLambda = (6 + 3*flxR) / (6 + 3*flxR + flxR^2) (see Levermore & Pomraning 1981)
L_rad = L_rad / flxLambda

With the resulting L_rad, determine the expected dT/dm by

d_Prad/dm = -kap*L_rad/(clight*area^2) – see, e.g., K&W (5.12)

use_dPrad_dm_form_of_T_gradient_eqn = .false.
use_flux_limiting_with_dPrad_dm_form = .false.
use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false.

Hydrodynamic drag#

drag_coefficient#

min_q_for_drag#

only when v_flag = .true.. Adjusts both v and energy transfer from kinetic to thermal. only for v(k) when q(k) > min_q_for_drag. kill off fraction of v = drag_coefficient (i.e. set to 1 to keep v near 0) useful for preventing the development radial pulsations during advanced burning in massive stars and AGB stars.

Under certain circumstances we will not have drag in the surface k=1 zone To force the drag term to be on in the outer zone you must enable one of the following surface boundary conditions: use_momentum_outer_BC, use_zero_Pgas_outer_BC, or use_fixed_Psurf_outer_BC

If use_drag_energy = .true. adds the work done by the drag force to the energy balance. Since the drag force is proportional to the velocity which is subject to numerical spikes and oscillations, you may not want to use the work done by it in the model If set to .false. but drag_coefficient is nonzero, the drag force will be applied in the momentum equation (as a numerical trick to damp spurious velocities) but not in the energy equation.

use_drag_energy = .true.
drag_coefficient = 0d0
min_q_for_drag = 0d0

v_drag_factor#

v_drag#

q_for_v_drag_full_off#

q_for_v_drag_full_on#

Only when u_flag = .true.. Adds a pseudo drag term of the form -v_drag_factor*(v-v_drag)^2/r, can be used damp velocities in outer layers of a star, useful for smoothing out spurious shocks colliding in ejected layers. Effect is full on for q>q_for_v_drag_full_on and full off for q < q_for_v_drag_full_off.

v_drag_factor = 0d0
v_drag = 0d0
q_for_v_drag_full_off = 0.95d0
q_for_v_drag_full_on = 0.96d0

Rayleigh-Taylor Instability: hydro drag terms#

for hydro comparison tests (e.g., Sedov)

RTI_A#

RTI_B#

RTI_C#

RTI_D#

RTI_C_X_factor#

RTI_C_X0#

RTI_max_alpha#

RTI_min_dm_behind_shock_for_full_on#

RTI_dm_for_center_alpha_nondecreasing#

RTI_energy_floor#

RTI_D_mix_floor#

RTI_min_m_for_D_mix_floor#

RTI_log_max_boost#

RTI_m_full_boost#

RTI_m_no_boost#

Note that these parameters are not exactly the same as used by Paul Duffell. His calibrated D is 2, where mesa has default D = 3 (see MESA IV). Users should try various values since the choice is not clear cut.

RTI_A = 1d-3
RTI_B = 2.5d0
RTI_C = 0.2d0
RTI_D = 3d0
RTI_C_X0_frac = 0.9d0
RTI_C_X_factor = 0d0
RTI_max_alpha = 0.5d0
RTI_min_dm_behind_shock_for_full_on = 0d0
RTI_dm_for_center_eta_nondecreasing = 0.02d0
RTI_energy_floor = 0d0
RTI_D_mix_floor = 0d0
RTI_min_m_for_D_mix_floor = 0d0
RTI_log_max_boost = 3d0
RTI_m_full_boost = 4d0
RTI_m_no_boost = 5d0

retry_for_v_above_clight#

If .true., a retry will be triggered at the end of a step if the maximum velocity exceeds the speed of light. If .false., only a warning is printed.

retry_for_v_above_clight = .true.

solver controls#

the following is from a response on mesa-users to a question about controls for solver tolerances:

The “residual” is the left over difference between the left and right hand sides of the equation we are trying to solve. We do iterations to reduce that, but we are limited by the non-linearity of the problem and the quality of the estimates for the derivatives.

The “correction” is the change in the primary variable that is calculated using good-old Newton’s rule in multiple dimensions — so Jacobian and residuals give a correction that would make the next residual vanish if the problem were linear and the Jacobian was exact, neither of which are true. So the best we can hope for is that the corrections will get smaller next time.

The “norm” is the average; the “max” is the max. Sometimes you mainly care about the norm and will accept a few outliers. But sometimes you don’t want any really bad outliers, so you want to set a low limit for the max residual or correction as well as the norm.

You might want to try for several iterations with strict tolerances, and then relax them if things are still not converged. For example, you might be willing to live with the larger tolerances, but you’d like to give it a good try at the smaller ones before switching. Also, you might be willing to settle for any-old residual if the corrections have become small enough. You can do that too by relaxing the residual tolerances after a few iterations.

Hope that at least helps with the nomenclature.

I agree with Frank that you should consider the effects of smaller timesteps and more grid points as your main technique — tightening up the tolerances for the solver won’t help if you are taking timesteps that are too large or if you have inadequate grid resolution.

tol_correction_norm#

tol_max_correction#

“Correction” for variable x(i,k) is scaled change, dx(i,k)/xscale(i,k). these tolerances are for the magnitude of the scaled corrections.

tol_correction_norm = 3d-5
tol_max_correction = 3d-3

tol_correction_high_T_limit#

For very late stages of massive star evolution, need to relax tolerances. If max T >= this limit, switch scaling factors.

tol_correction_high_T_limit = 1d9

tol_correction_norm_high_T#

tol_max_correction_high_T#

Above tol_correction_high_T_limit use these scaling factors.

tol_correction_norm_high_T = 3d-3
tol_max_correction_high_T = 3d-1

tol_correction_extreme_T_limit#

For very late stages of massive star evolution, need to relax tolerances. If center T >= this limit, switch scaling factors.

tol_correction_extreme_T_limit = 6d9

tol_correction_norm_extreme_T#

tol_max_correction_extreme_T#

For very late stages of massive star evolution, need to relax tolerances. If center T >= this limit, switch scaling factors.

tol_correction_norm_extreme_T = 8d-3
tol_max_correction_extreme_T = 8d-1

tol_bad_max_correction#

if max_correction > tol_max_correction and no more iterations allowed, then still accept the solution if max_correction <= tol_bad_max_correction. but if max_correction > tol_bad_max_correction, then reject the solution.

tol_bad_max_correction = 0d0

bad_max_correction_series_limit#

If have this many steps in a row with max_correction > tol_max_correction, then do a retry with a smaller timestep.

bad_max_correction_series_limit = 2

relax_use_gold_tolerances#

relax_use_gold_tolerances = .false.

relax_solver_iters_timestep_limit#

relax_tol_correction_norm#

relax_tol_max_correction#

relax_tol_residual_norm1#

relax_tol_max_residual1#

relax_iter_for_resid_tol2#

relax_tol_residual_norm2#

relax_tol_max_residual2#

relax_iter_for_resid_tol3#

relax_tol_residual_norm3#

relax_tol_max_residual3#

relax_maxT_for_gold_tolerances#

For use during relax operations. Only used if /= 0.

relax_solver_iters_timestep_limit = 0
relax_tol_correction_norm = 0d0
relax_tol_max_correction = 0d0
relax_tol_residual_norm1 = 0d0
relax_tol_max_residual1 = 0d0
relax_iter_for_resid_tol2 = 3
relax_tol_residual_norm2 = 0d0
relax_tol_max_residual2 = 0d0
relax_iter_for_resid_tol3 = 0
relax_tol_residual_norm3 = 0d0
relax_tol_max_residual3 = 0d0
relax_maxT_for_gold_tolerances = -1d0

include_L_in_correction_limits#

include_v_in_correction_limits#

include_u_in_correction_limits#

include_w_in_correction_limits#

These variables can be excluded from calculation of correction norm and max.

include_L_in_correction_limits = .true.
include_v_in_correction_limits = .true.
include_u_in_correction_limits = .true.
include_w_in_correction_limits = .true.

max_X_for_conv_timescale#

min_X_for_conv_timescale#

max_q_for_conv_timescale#

min_q_for_conv_timescale#

max_q_for_QHSE_timescale#

min_q_for_QHSE_timescale#

max_X_for_conv_timescale = 1d0
min_X_for_conv_timescale = 0d0
max_q_for_conv_timescale = 1d0
min_q_for_conv_timescale = 0d0
max_q_for_QHSE_timescale = 1d0
min_q_for_QHSE_timescale = 0d0

correction_xa_limit#

Ignore correction to abundance when calculating correction norm and max if current mass fraction is less than this limit.

correction_xa_limit = 5d-3

xa_scale#

Scaling for abundance variables is max(xa_scale, current mass fraction).

xa_scale = 1d-5

tol_residual_norm1#

tol_max_residual1#

iter_for_resid_tol2#

“residual” for equation is the difference between left and right sides use tol_residual_norm1 & tol_max_residual1 at iteration number iter_for_resid_tol2, switch to next tolerances.

tol_residual_norm1 = 1d-10
tol_max_residual1 = 1d-9
iter_for_resid_tol2 = 6

tol_residual_norm2#

tol_max_residual2#

iter_for_resid_tol3#

Use tol_residual_norm2 & tol_max_residual2 these apply starting at iteration number iter_for_resid_tol2. at iteration number iter_for_resid_tol3, switch to next tolerances.

tol_residual_norm2 = 1d90
tol_max_residual2 = 1d90
iter_for_resid_tol3 = 15

tol_residual_norm3#

tol_max_residual3#

Use tol_residual_norm3 & tol_max_residual3 these apply starting at iteration number iter_for_resid_tol3.

tol_residual_norm3 = 1d90
tol_max_residual3 = 1d90

If things get worse from one iteration to next, give up. The following are the limits that define “getting worse enough to stop”.

corr_norm_jump_limit#

If correction norm increases by this factor or more, quit.

corr_norm_jump_limit = 1d99

max_corr_jump_limit#

If correction max increases by this factor or more, quit.

max_corr_jump_limit = 1d6

resid_norm_jump_limit#

If residual norm increases by this factor or more, quit.

resid_norm_jump_limit = 1d99

max_resid_jump_limit#

If residual max increases by this factor or more, quit.

max_resid_jump_limit = 1d6

convergence_ignore_equL_residuals#

convergence_ignore_equL_residuals = .false.

convergence_ignore_alpha_RTI_residuals#

convergence_ignore_alpha_RTI_residuals = .false.

trace_solver_damping#

Send solver damping data to screen.

trace_solver_damping = .false.

do_normalize_dqs_as_part_of_set_qs#

normalize_dqs destroys bit-for-bit read as inverse of write for models. ok for create pre ms etc., but not for read model create_pre_ms calls normalize_dqs even if this flag is false.

do_normalize_dqs_as_part_of_set_qs = .false.

use_DGESVX_in_bcyclic use_equilibration_in_DGESVX report_min_rcond_from_DGESXV

FOR DEBUGGING ONLY. NOT FOR GENERAL USE.

use_DGESVX_in_bcyclic = .false.
use_equilibration_in_DGESVX = .false.
report_min_rcond_from_DGESXV = .false.

solver_max_tries_before_reject#

Max number solver iterations before give up.

solver_max_tries_before_reject = 25

max_tries1#

Max tries on 1st model.

max_tries1 = 250

max_tries_for_retry#

Normal number of retries.

max_tries_for_retry = 25

max_tries_after_5_retries#

Increase number of tries after 5 failed ones.

max_tries_after_5_retries = 35

max_tries_after_10_retries#

Increase number of tries after 10 failed ones.

max_tries_after_10_retries = 50

max_tries_after_20_retries#

Increase number of tries after 20 failed ones.

max_tries_after_20_retries = 75

retry_limit#

Only use if > 0. In case the solver fails for some reason, it will retry with a smaller timestep. It does up to this many retries for the current step before terminating.

retry_limit = 100

redo_limit#

Only use if > 0. Do up to this many redo’s for the current step before terminating.

redo_limit = 100

solver_itermin#

Use at least this many iterations in solver for hydro solve.

solver_itermin = 2

solver_itermin_until_reduce_min_corr_coeff#

Use at least this many iterations in solver before try using small min_corr_coeff

solver_itermin_until_reduce_min_corr_coeff = 8

solver_reduced_min_corr_coeff#

For use with solver_itermin_for_reduce_min_corr_coeff.

solver_reduced_min_corr_coeff = 0.1d0

tiny_corr_coeff_limit#

scale_correction_norm#

corr_param_factor#

scale_max_correction#

ignore_too_large_correction#

corr_coeff_limit#

tiny_corr_factor#

ignore_min_corr_coeff_for_scale_max_correction#

ignore_species_in_max_correction#

num_times_solver_reuse_mtx#

see star/private/star_solver for info about these

tiny_corr_coeff_limit = 100
scale_correction_norm = 0.1d0
corr_param_factor = 10
scale_max_correction = 1d99
ignore_too_large_correction = .false.
corr_coeff_limit = 1d-2
tiny_corr_factor = 2
ignore_min_corr_coeff_for_scale_max_correction = .false.
ignore_species_in_max_correction = .false.
num_times_solver_reuse_mtx = 0

min_xa_hard_limit#

min_xa_hard_limit_for_highT#

If solver produces mass fraction < this limit, then reject the trial solution. Can optionally relax this limit at high T.

min_xa_hard_limit = -1d-5
min_xa_hard_limit_for_highT = -3d-5

logT_max_for_xa_hard_limit#

Use min_xa_hard_limit for center logT <= this.

logT_max_for_min_xa_hard_limit = 9.49d0

logT_min_for_xa_hard_limit_for_highT#

Use min_xa_hard_limit_for_highT for center logT >= this. Linear interpolate in logT for intermediate center temperatures.

logT_min_for_min_xa_hard_limit_for_highT = 9.51d0

sum_xa_hard_limit#

sum_xa_hard_limit_for_highT#

If solver produces any cell with abs(sum(xa)-1) > this limit, then reject the trial solution. Can optionally relax this limit at high T.

sum_xa_hard_limit = 5d-4
sum_xa_hard_limit_for_highT = 1d-3

logT_max_for_sum_xa_hard_limit#

Use sum_xa_hard_limit for center logT <= this.

logT_max_for_sum_xa_hard_limit = 9.40d0

logT_min_for_sum_xa_hard_limit_for_highT#

Use sum_xa_hard_limit_for_highT for center logT >= this. Linear interpolate in logT for intermediate center temperatures.

logT_min_for_sum_xa_hard_limit_for_highT = 9.44d0

do_solver_damping_for_neg_xa#

If true, uniformly reduce solver corrections if necessary to avoid neg abundances.

do_solver_damping_for_neg_xa = .true.

scale_max_correction_for_negative_surf_lum#

max_frac_for_negative_surf_lum#

If true, then scales the correction factor in a Newton iteration to prevent the surface from reaching a negative luminosity. If an iteration would require s% L(1) to become negative, then the correction is scaled such that the change in surface luminosity is -max_frac_for_negative_surf_lum*s% L(1)

scale_max_correction_for_negative_surf_lum = .false.
max_frac_for_negative_surf_lum = 0.8

min_chem_eqn_scale#

min_chem_eqn_scale = 1d0

hydro_mtx_max_allowed_{abs}{dlogT | dlogRho | logT | logRho}#

Force retry with smaller timestep if hydro solves change T or Rho by too much or make them too large.

hydro_mtx_max_allowed_abs_dlogT = 99d0
hydro_mtx_max_allowed_abs_dlogRho = 99d0
min_logT_for_hydro_mtx_max_allowed = -1d99
hydro_mtx_max_allowed_logT = 12d0
hydro_mtx_max_allowed_logRho = 12d0
hydro_mtx_min_allowed_logT = 1d0
hydro_mtx_min_allowed_logRho = -1d2

level 1 of gold tolerances for solver solver

use_gold_tolerances#

steps_before_use_gold_tolerances#

gold_solver_iters_timestep_limit#

maxT_for_gold_tolerances#

gold_tol_residual_norm1#

gold_iter_for_resid_tol2#

gold_tol_residual_norm2#

gold_tol_max_residual2#

gold_iter_for_resid_tol3#

gold_tol_residual_norm3#

gold_tol_max_residual3#

use_gold_tolerances = .true.
steps_before_use_gold_tolerances = -1

if >= 0, then after this many steps in run, act as if use_gold_tolerances true this allows a delay before turning on gold tolerances NOTE: if using steps_before_use_gold_tolerances >= 0, then set use_gold_tolerances = false

maxT_for_gold_tolerances = 1d99
gold_tol_residual_norm1 = 1d-11
gold_tol_max_residual1 = 1d-9
gold_iter_for_resid_tol2 = 5
gold_tol_residual_norm2 = 1d-8
gold_tol_max_residual2 = 1d-6
gold_iter_for_resid_tol3 = 10
gold_tol_residual_norm3 = 1d-6
gold_tol_max_residual3 = 1d-4
gold_solver_iters_timestep_limit = 14

level 2 of gold tolerances for solver solver - tighter than level 1

use_gold2_tolerances#

steps_before_use_gold2_tolerances#

gold2_solver_iters_timestep_limit#

gold2_tol_residual_norm1#

gold2_iter_for_resid_tol2#

gold2_tol_residual_norm2#

gold2_tol_max_residual2#

gold2_iter_for_resid_tol3#

gold2_tol_residual_norm3#

gold2_tol_max_residual3#

use_gold2_tolerances = .false.
steps_before_use_gold2_tolerances = -1

if >= 0, then after this many steps in run, act as if use_gold2_tolerances true this allows a delay before turning on level 2 gold tolerances NOTE: if using steps_before_use_gold2_tolerances >= 0, then set use_gold2_tolerances = false

gold2_tol_residual_norm1 = 1d-11
gold2_tol_max_residual1 = 1d-9
gold2_iter_for_resid_tol2 = 5
gold2_tol_residual_norm2 = 1d-10
gold2_tol_max_residual2 = 1d-8
gold2_iter_for_resid_tol3 = 10
gold2_tol_residual_norm3 = 1d-8
gold2_tol_max_residual3 = 1d-5
gold2_solver_iters_timestep_limit = 18

include_rotation_in_total_energy#

previously called include_rotation_in_energy_error_report#

include_rotation_in_total_energy = .false.

artificial viscosity

use_Pvsc_art_visc#

Pvsc_cq#

Pvsc_zsh#

Pvsc is artificial pressure to push back against compression this is the form of artificial viscosity used in RSP if using this, do not set use_artificial_viscosity true.

artificial viscosity controls for the equations see: Appendix C in Stellingwerf 1975 http://adsabs.harvard.edu/abs/1975ApJ…195..441S. In principle, for not too-non-adiabatic convective models artificial viscosity is not needed or should be very small. Hence a large cut-off parameter below (in purely radiative models the default value for cut-off was 0.01)

zsh > 0 delays onset of artificial viscosity can eliminate most/all interior dissipation while still providing for extreme cases. using this parameter the dependence of limiting amplitude on cq is very weak.

use_Pvsc_art_visc = .false.
Pvsc_cq = 4.0d0
Pvsc_zsh = 0.1d0

use_artificial_viscosity#

use_artificial_viscosity has been replaced by use_Pvsc_art_visc.

split burn#

op_split_burn#

op_split_burn = .false.

op_split_burn_min_T#

op_split_burn_eps#

op_split_burn_odescal#

Only do op_split_burn in cells with T >= this limit at start of step.

op_split_burn_min_T = 2d9
op_split_burn_eps = 1d-5
op_split_burn_odescal = 1d-5

op_split_burn_eps_nuc_infall_limit#

turn off op_split_burn nuclear burning if max infall speed exceeds this limit (cm/s).

op_split_burn_eps_nuc_infall_limit = 1d99

timestep controls#

The terminal output during evolution includes a short string for the dt_limit. This is to give you some indication of what is limiting the time steps. Here’s a dictionary mapping those terminal strings to the corresponding control parameters. (There is a similar table in mesa/binary/defaults/binary_controls.defaults.)

 terminal output       related parameter
'burn steps'           burn_steps_limit
'Lnuc'                 delta_lgL_nuc_limit
'Lnuc_cat'             delta_lgL_nuc_cat_limit
'Lnuc_H'               delta_lgL_H_limit
'Lnuc_He'              delta_lgL_He_limit
'lgL_power_phot'       delta_lgL_power_photo_limit
'Lnuc_z'               delta_lgL_z_limit
'bad_X_sum'            (solver found bad mass sum)
'dL/L'                 dL_div_L_limit
'dX'                   dX_limit
'dX/X'                 dX_div_X_limit
'dX_nuc_drop'          dX_nuc_drop_limit
'delta mdot'           delta_mdot_limit
'delta total J'        delta_lg_total_J_limit
'delta_HR'             delta_HR_limit
'delta_mstar'          delta_lg_star_mass_limit
'diff iters'           diffusion_iters_limit
'diff steps'           diffusion_steps_limit
'min_dr_div_cs'        dt_div_min_dr_div_cs_limit
'dt_collapse'          dt_div_dt_cell_collapse_limit
'eps_nuc_cntr'         delta_log_eps_nuc_cntr_limit
'error rate'           limit_for_log_rel_rate_in_energy_conservation
'highT del Ye'         delta_Ye_highT_limit
'hold'                 (recent retry, so no increase in dt)
'lgL'                  delta_lgL_limit
'lgP'                  delta_lgP_limit
'lgP_cntr'             delta_lgP_cntr_limit
'lgR'                  delta_lgR_limit
'lgRho'                delta_lgRho_limit
'lgRho_cntr'           delta_lgRho_cntr_limit
'lgT'                  delta_lgT_limit
'lgT_cntr'             delta_lgT_cntr_limit
'lgT_max'              delta_lgT_max_limit
'lgT_max_hi_T'         delta_lgT_max_at_high_T_limit
'lgTeff'               delta_lgTeff_limit
'dX_div_X_cntr'        delta_dX_div_X_cntr_limit
'lg_XC_cntr'           delta_lg_XC_cntr_limit
'lg_XH_cntr'           delta_lg_XH_cntr_limit
'lg_XHe_cntr'          delta_lg_XHe_cntr_limit
'lg_XNe_cntr'          delta_lg_XNe_cntr_limit
'lg_XO_cntr'           delta_lg_XO_cntr_limit
'lg_XSi_cntr'          delta_lg_XSi_cntr_limit
'XC_cntr'              delta_XC_cntr_limit
'XH_cntr'              delta_XH_cntr_limit
'XHe_cntr'             delta_XHe_cntr_limit
'XNe_cntr'             delta_XNe_cntr_limit
'XO_cntr'              delta_XO_cntr_limit
'XSi_cntr'             delta_XSi_cntr_limit
'log_eps_nuc'          delta_log_eps_nuc_limit
'max_dt'               max_years_for_timestep
'neg_mass_frac'        (solver found neg mass frac)
'adjust_J_q'           adjust_J_q_limit
'solver iters'         solver_iters_timestep_limit
'rel_E_err'            limit_for_rel_error_in_energy_conservation
'varcontrol'           varcontrol_target
'max increase'         max_timestep_factor or max_timestep_factor_at_high_T
'max decrease'         min_timestep_factor
'retry'                (just did a retry)
'b_****'               see binary/defaults/binary_controls.defaults

time_delta_coeff#

time_delta_coeff - smaller forces smaller timesteps giving better time resolution. multiplier for all real number timestep limits and hard limits. does not apply to integer valued limits such as

  • solver_iters_timestep_limit

  • burn_steps_limit

  • diffusion_steps_limit

  • diffusion_iters_limit

does not apply to varcontrol_target. analogous to mesh_delta_coeff for better spatial resolution.

time_delta_coeff = 1d0

max_timestep#

In seconds. max_timestep <= 0 means no upper limit.

max_timestep = 0

max_years_for_timestep#

max_years_for_timestep <= 0 means no upper limit. Note: max_timestep is the control that is used by most of the code. max_years_for_timestep is just provided as a convenience. At the start of each step, the evolve routine checks to see if max_years_for_timestep > 0, and if so, it sets max_timestep = max_years_for_timestep*secyer.

max_years_for_timestep = 0

max_timestep_hi_T_limit#

If max T >= this, then switch to hi_T_max_years_for_timestep. Ignore if <= 0.

max_timestep_hi_T_limit = -1

hi_T_max_years_for_timestep#

Max years for timestep if max_timestep_hi_T_limit is active.

hi_T_max_years_for_timestep = 0

min_timestep_factor#

Lower limit for ratio of new timestep to previous timestep. i.e., allow dt to get smaller by no more than this factor – 0 means no limit.

min_timestep_factor = 0.8d0

force_timestep#

In seconds. force_timestep <= 0 means no forced timestep.

force_timestep = 0

force_timestep_years#

Note: force_timestep is the control that is used by most of the code. force_timestep_years is just provided as a convenience. At the start of each step, the evolve routine checks if force_timestep_years > 0, and if so, it sets force_timestep = force_timestep_years*secyer.

force_timestep_years = 0

force_timestep_min#

In seconds. force_timestep_min <= 0 means no forced lower limit.

force_timestep_min = 0

force_timestep_min_years#

Note: force_timestep_min is the control that is used by most of the code. force_timestep_min_years is just provided as a convenience. At the start of each step, the evolve routine checks if force_timestep_min_years > 0, and if so, it sets force_timestep_min = force_timestep_min_years*secyer.

force_timestep_min_years = 0

force_timestep_min_factor#

If dt is < force_timestep_min, then replace dt by min(dt*force_timestep_min_factor, force_timestep_min)

force_timestep_min_factor = 2d0

max_timestep_factor#

max_timestep_factor_at_high_T#

min_logT_for_max_timestep_factor_at_high_T#

Upper limit for ratio of new timestep to previous timestep. i.e., allow dt to get larger by no more than this factor – 0 means no limit. use max_timestep_factor_at_high_T when max logT > min_logT_for_max_timestep_factor_at_high_T.

max_timestep_factor = 1.2d0
max_timestep_factor_at_high_T = 1.1d0
min_logT_for_max_timestep_factor_at_high_T = 1d99

timestep_factor_for_retries#

Before retry, decrease dt by this.

timestep_factor_for_retries = 0.5d0

retry_hold#

No increases in timestep for retry_hold steps after a retry.

retry_hold = 1

neg_mass_fraction_hold#

No increases in timestep for neg_mass_fraction_hold steps after a retry caused by a negative mass fraction.

neg_mass_fraction_hold = 2

timestep_dt_factor = 0.9#

dt reduction factor exceed timestep limits.

timestep_dt_factor = 0.9d0

use_dt_low_pass_controller#

Enable low pass filter for smoother timestep variations.

use_dt_low_pass_controller = .true.

varcontrol_target#

This is the target value for relative variation in the structure from one model to the next. The default timestep adjustment is to increase or reduce the timestep depending on whether the actual variation was smaller or greater than this value.

varcontrol_target = 1d-3

min_allowed_varcontrol_target#

The run will terminate if varcontrol_target < min_allowed_varcontrol_target It is not usually a good idea to reduce this. Instead, use time_delta_coeff to do time resolution convergence studies.

min_allowed_varcontrol_target = 1d-4

varcontrol_dt_limit_ratio_hard_max#

varcontrol_dt_limit_ratio is the actual varcontrol value divided by the target. if that ratio exceeds this limit, then retry with a smaller timestep. this let’s you prevent large changes from happening in a single step.

varcontrol_dt_limit_ratio_hard_max = 1d99

never_skip_hard_limits#

If true, then don’t skip hard limits even immediately after a retry.

never_skip_hard_limits = .true.

relax_hard_limits_after_retry#

If true, then don’t enforce hard limits immediately after a retry.

relax_hard_limits_after_retry = .false.

limits based on iterations required by various solvers

solver_iters_timestep_limit#

If solver solve uses more solver_iterations than this, reduce the next timestep. NOTE: when using gold tolerances, set gold_solver_iters_timestep_limit.

solver_iters_timestep_limit = 7

burn_steps_limit#

If burn solver uses more steps than this, reduce the next timestep.

burn_steps_limit = 10000

burn_steps_hard_limit#

If burn solver uses more steps than this, retry.

burn_steps_hard_limit = 20000

diffusion_steps_limit#

If diffusion solver uses more steps than this, reduce the next timestep.

diffusion_steps_limit = 500

diffusion_steps_hard_limit#

If diffusion solver uses more steps than this, retry.

diffusion_steps_hard_limit = 700

diffusion_iters_limit#

If use a total number of iters > this, reduce the next timestep.

diffusion_iters_limit = 600

diffusion_iters_hard_limit#

If use a total number of iters > this, retry.

diffusion_iters_hard_limit = 800

limits based on max decrease in mass fraction at any location in star

dX_mix_dist_limit#

Option to ignore decreases in abundance in non-mixed cells near mixing boundaries. Ignore abundance changes if nearest mixing boundary is closer than this in Msun units. This applies to dX, and dX_div_X limits.

dX_mix_dist_limit = 1d-4

dX_limit_species#

Specify which species the dX_limit, dX_div_X_limit, etc array entries apply to. These are limits on magnitude of decrease in any cell abundance during a single timestep. dX here is abs(xa(j,k) - xa_old(j,k)) for any cell k and all species j, eg 'h1', 'he4', etc. Special ‘species’ 'X' (any hydrogen), 'Y' (any helium) and 'Z' (any metals) are allowed here. E.g. 'Z' will trigger the timestep control if any metal isotope abundance individually satisfies the conditions below. Considers all cells except where have convective mixing.

dX_limit_species(1) = 'h1'
dX_limit_species(2) = 'he4'
dX_limit_species(3:) = ''

dX_limit_min_X#

dX limits only apply where xa(j,k) >= this limit.

dX_limit_min_X(:) = 1d99

dX_limit#

If max dX is greater than this, reduce the next timestep by dX_limit/max_dX.

dX_limit(:) = 1d99

dX_hard_limit#

If max dX is greater than this, retry with smaller timestep.

dX_hard_limit(:) = 1d99

dX_decreases_only#

If true, then only consider decreases in abundance. dX_decreases_only applies to dX_div_X also.

dX_decreases_only(:) = .true.

Limit on magnitude of relative decrease in any cell abundance. dX_div_X here is abs(xa(j,k) - xa_old(j,k))/xa(j,k) for any cell k and any species j. Considers all cells except where have convective mixing.

dX_div_X_limit_min_X#

dX_div_X limits only apply where xa(j,k) >= this limit.

dX_div_X_limit_min_X(1) = 1d-3
dX_div_X_limit_min_X(2) = 1d-3
dX_div_X_limit_min_X(3:) = 1d99

dX_div_X_limit#

If max dX_div_X is greater than this, reduce the next timestep by dX_limit/max_dX.

dX_div_X_limit(1) = 0.9d0
dX_div_X_limit(2) = 0.9d0
dX_div_X_limit(3:) = 1d99

dX_div_X_hard_limit#

If max dX_div_X is greater than this, retry with smaller timestep.

dX_div_X_hard_limit(:) = 1d99

dX_div_X_at_high_T_limit#

dX_div_X_at_high_T_hard_limit#

dX_div_X_at_high_T_limit_lgT_min#

dX_div_X_at_high_T_limit(:) = 1d99
dX_div_X_at_high_T_hard_limit(:) = 1d99
dX_div_X_at_high_T_limit_lgT_min(:) = 1d99

Limits on max drop in abundance mass fraction from burning with possible mixing inflow. This considers both nuclear reactions and offsetting effect of mixing inflow.

dX_nuc_drop_min_X_limit#

dX_nuc_drop_limit only for X > dX_nuc_drop_min_X_limit. note that this is the abundance of the species, not the H1 abundance for the cell. i.e., if species abundance is below this limit, then ignore it for the dX_nuc_drop limit.

dX_nuc_drop_min_X_limit = 1d-4

dX_nuc_drop_max_A_limit#

dX_nuc_drop_limit only for species with A <= dX_nuc_drop_max_A_limit.

dX_nuc_drop_max_A_limit = 52

dX_nuc_drop_limit_at_high_T#

Negative means use value for dX_nuc_drop_limit, else use this limit when max logT > 9.45.

dX_nuc_drop_limit_at_high_T = -1

dX_nuc_drop_limit#

If max dX_nuc_drop is greater than dX_nuc_drop_limit, reduce the next timestep by dX_nuc_drop_limit/max_dX_nuc_drop.

dX_nuc_drop_limit = 5d-2

dX_nuc_drop_hard_limit#

If max dX_nuc_drop is greater than dX_nuc_drop_hard_limit, retry with smaller timestep.

dX_nuc_drop_hard_limit = 1d99

dX_nuc_drop_min_yrs_for_dt#

Don’t let dX_nuc_drop change dt to smaller than this.

dX_nuc_drop_min_yrs_for_dt = 1d-9

limits based on relative changes in variables L, P, Rho, T, R, eps_nuc#

limit on magnitude of relative change in L at any grid point

dL_div_L = abs(L(k) - L_old(k))/L(k)

dL_div_L_limit#

If max abs dL_div_L is greater than this, reduce the next timestep.

dL_div_L_limit = -1

dL_div_L_hard_limit#

If max abs dL_div_L is greater than this, retry with smaller timestep.

dL_div_L_hard_limit = -1

dL_div_L_limit_min_L#

In Lsun units. dL_div_L limits only apply where L(k) >= Lsun*dL_limit_min_L

dL_div_L_limit_min_L = 1d99

delta_lgP_limit#

Limit for magnitude of max change in log10 total pressure in any cell.

delta_lgP_limit = 1

delta_lgP_hard_limit#

If max delta_lgP is greater than delta_lgP_hard_limit, retry with smaller timestep.

delta_lgP_hard_limit = -1

delta_lgP_limit_min_lgP#

delta_lgP_limit limits only apply where log10_P(k) >= delta_lgP_limit_min_lgP

delta_lgP_limit_min_lgP = 1d99

delta_lgRho_limit#

Limit for magnitude of max change in log10 density in any cell.

delta_lgRho_limit = 1

delta_lgRho_hard_limit = -1#

If max delta_lgRho is greater than delta_lgRho_hard_limit, retry with smaller timestep.

delta_lgRho_hard_limit = -1

delta_lgRho_limit_min_lgRho#

delta_lgRho_limit limits only apply where log10_Rho(k) >= delta_lgRho_limit_min_lgRho.

delta_lgRho_limit_min_lgRho = 1d99

delta_lgT_limit#

Limit for magnitude of max change in log10 temperature in any cell.

delta_lgT_limit = 0.5d0

delta_lgT_hard_limit#

If max delta_lgT is greater than delta_lgT_hard_limit, retry with smaller timestep.

delta_lgT_hard_limit = -1

delta_lgT_limit_min_lgT#

delta_lgT_limit limits only apply where log10_T(k) >= delta_lgT_limit_min_lgT.

delta_lgT_limit_min_lgT = 1d99

delta_lgE_limit#

Limit for magnitude of max change in log10 internal energy in any cell.

delta_lgE_limit = 0.1d0

delta_lgE_hard_limit#

If max delta_lgE is greater than delta_lgE_hard_limit, retry with smaller timestep.

delta_lgE_hard_limit = -1

delta_lgE_limit_min_lgE#

delta_lgE_limit limits only apply where log10(E(k)) >= delta_lgE_limit_min_lgE.

delta_lgE_limit_min_lgE = 1d99

delta_lgR_limit#

Limit for magnitude of max change in log10 radius at any cell boundary.

delta_lgR_limit = 0.5d0

delta_lgR_hard_limit#

If max delta_lgR is greater than delta_lgR_hard_limit, retry with smaller timestep.

delta_lgR_hard_limit = -1

delta_lgR_limit_min_lgR#

delta_lgR_limit limits only apply where log10_R(k) >= delta_lgR_limit_min_lgR.

delta_lgR_limit_min_lgR = 1d99

delta_Ye_highT_limit#

Limit for magnitude of max change in Ye in high T cells.

delta_Ye_highT_limit = 99

Limit testing for max delta_ye to cells with T >= minT_for_highT_Ye_limit If this high T max delta_Ye is greater than delta_Ye_highT_limit, reduce the next timestep by delta_Ye_highT_limit/max_delta_Ye.

delta_Ye_highT_hard_limit = -1

minT_for_highT_Ye_limit#

Limit testing for max delta_ye to cells with T >= minT_for_highT_Ye_limit. If this high T max delta_Ye is greater than delta_Ye_highT_limit, retry with smaller timestep.

minT_for_highT_Ye_limit = 7d9

delta_log_eps_nuc_limit#

Limit for magnitude of max change in log10 eps_nuc in any cell. Only applies to increases in non-convective zones.

delta_log_eps_nuc_limit = -1

delta_log_eps_nuc_hard_limit#

If max delta_log_eps_nuc is greater than delta_log_eps_nuc_hard_limit, retry with smaller timestep.

delta_log_eps_nuc_hard_limit = -1

limits based on integrated power at each point for each category of nuclear reaction#

lgL_nuc_cat = nuclear reaction energy release for a particular category of reaction (Lsun units). Energy release here excludes neutrinos.

delta_lgL_nuc_cat_limit#

Limit for magnitude of change in lgL_nuc for category.

delta_lgL_nuc_cat_limit = -1

delta_lgL_nuc_cat_hard_limit#

If max delta is greater than delta_lgL_nuc_cat_hard_limit, retry with smaller timestep.

delta_lgL_nuc_cat_hard_limit = -1

lgL_nuc_cat_burn_min#

Ignore changes in lgL_nuc for category if value is less than this.

lgL_nuc_cat_burn_min = -1

lgL_nuc_mix_dist_limit#

Ignore if nearest boundary is closer than this. Ignore changes in lgL in cells near mixing boundaries.

lgL_nuc_mix_dist_limit = 1d-6

L_H_burn = integrated power at surface from PP and CNO (in Lsun units)

values for lgL_H are log10(max(1, L_H_burn))

delta_lgL_H_limit#

limit for magnitude of change in lgL_H

delta_lgL_H_limit = -1

delta_lgL_H_hard_limit#

if max delta is greater than delta_lgL_H_hard_limit, retry with smaller timestep

delta_lgL_H_hard_limit = -1

lgL_H_burn_min#

ignore changes in lgL_H if value is less than this

lgL_H_burn_min = 1.5d0

lgL_H_drop_factor#

when L_H is dropping, multiply limits by this factor

lgL_H_drop_factor = 1

lgL_H_burn_relative_limit#

ignore changes in lgL_H if max(lgL_He,lgL_z) - lgL_H > this

lgL_H_burn_relative_limit = 3

L_He_burn = integrated power at surface from triple alpha (in Lsun units)

values for lgL_He are log10(max(1, L_He_burn))

delta_lgL_He_limit#

Limit for magnitude of change in lgL_He.

delta_lgL_He_limit = 0.025d0

delta_lgL_He_hard_limit#

If max delta is greater than delta_lgL_He_hard_limit, retry with smaller timestep.

delta_lgL_He_hard_limit = -1

lgL_He_burn_min#

Ignore changes in lgL_He if value is less than this.

lgL_He_burn_min = 2.5d0

lgL_He_drop_factor#

When L_He is dropping, multiply limits by this factor.

lgL_He_drop_factor = 1

lgL_He_burn_relative_limit#

Ignore changes in lgL_He if max(lgL_H,lgL_z) - lgL_He > this.

lgL_He_burn_relative_limit = 3

L_z_burn = integrated power at surface from nuclear burning other than H, He, or C (in Lsun units) excluding photodistintegrations

values for lgL_z are log10(max(1, L_z_burn))

delta_lgL_z_limit#

Limit for magnitude of change in lgL_z.

delta_lgL_z_limit = -1

delta_lgL_z_hard_limit#

If max delta is greater than delta_lgL_z_hard_limit, retry with smaller timestep.

delta_lgL_z_hard_limit = -1

lgL_z_burn_min#

Ignore changes in lgL_z if value is less than this.

lgL_z_burn_min = 2.5d0

lgL_z_drop_factor#

When L_z is dropping, multiply limits by this factor.

lgL_z_drop_factor = 1

lgL_z_burn_relative_limit#

Ignore changes in lgL_z if max(lgL_H,lgL_He) - lgL_z > this.

lgL_z_burn_relative_limit = 3

limits based on total integrated power at surface for all nuclear reactions#

excluding photodistintegrations

L_nuc = nuclear reaction total energy release for all nuclear reactions (Lsun units)

delta_lgL_nuc_limit#

delta_lgL_nuc_hard_limit#

delta_lgL_nuc_at_high_T_limit#

delta_lgL_nuc_at_high_T_hard_limit#

delta_lgL_nuc_at_high_T_limit_lgT_min#

max_lgT_for_lgL_nuc_limit#

When L_nuc is dropping, multiply limits by lgL_nuc_drop_factor.

ignore changes in lgL_nuc if max logT > max_lgT_for_lgL_nuc_limit.

limit for magnitude of change in lgL_nuc retry if max delta is greater than delta_lgL_nuc_hard_limit, ignore changes in lgL_nuc if value is less than lgL_nuc_burn_min

if max logT is > delta_lgL_nuc_at_high_T_limit_lgT_min then

use delta_lgL_nuc_at_high_T_limit and delta_lgL_nuc_at_high_T_hard_limit

else

use delta_lgL_nuc_limit and delta_lgL_nuc_hard_limit

at extreme temperatures, can have numerical jitters in this as a result of almost cancelling large positive and negative contributions of forward and reverse rates such as photodisintegration and rebuilding of he4 so if doing advanced burning, may want to turn off lgL_nuc limits above some max logT and turn on lgL_power_photo limits at a somewhat smaller max logT.

delta_lgL_nuc_limit = -1
delta_lgL_nuc_hard_limit = -1
delta_lgL_nuc_at_high_T_limit = -1
delta_lgL_nuc_at_high_T_hard_limit = -1
delta_lgL_nuc_at_high_T_limit_lgT_min = 1d99
max_lgT_for_lgL_nuc_limit = 1d99
lgL_nuc_burn_min = 0.5d0
lgL_nuc_drop_factor = 10

limits based on total integrated power at surface for photodistintegrations#

L_photo = nuclear reaction total energy release for all photodistintegrations (Lsun units) note that photodistintegrations consume energy so the total released is negative. photodisintegrations become large during late burning at high temperatures. values for lgL_photo are based on L_by_category(iphoto)

delta_lgL_power_photo_limit#

delta_lgL_power_photo_hard_limit#

min_lgT_for_lgL_power_photo_limit#

lgL_power_photo_burn_min#

lgL_power_photo_drop_factor#

Limit for magnitude of change in lgL_photo. If max delta is greater than delta_lgL_power_photo_hard_limit, retry with smaller timestep. ignore delta_lgL_power_photo_limit when max logT < min_lgT_for_lgL_power_photo_limit. Ignore changes in lgL_photo if value is less than lgL_power_photo_burn_min. When L_photo is dropping, multiply limits by lgL_power_photo_drop_factor.

delta_lgL_power_photo_limit = -1
delta_lgL_power_photo_hard_limit = -1
min_lgT_for_lgL_power_photo_limit = 9d0
lgL_power_photo_burn_min = 10d0
lgL_power_photo_drop_factor = 10

limits based on changes at photosphere

delta_lgTeff_limit#

delta_lgTeff_hard_limit#

Limit for magnitude of max change in log10 temperature at photosphere.

delta_lgTeff_limit = 0.01d0
delta_lgTeff_hard_limit = -1

delta_lgL_limit_L_min#

delta_lgL_limit#

delta_lgL_hard_limit#

Limit for magnitude of change in log10(L_surf/Lsun). Only apply this limit when L_surf >= delta_lgL_limit_L_min (in Lsun units).

delta_lgL_limit_L_min = -100
delta_lgL_limit = 0.1d0
delta_lgL_hard_limit = -1

dt_div_min_dr_div_cs_limit#

dt_div_min_dr_div_cs_hard_limit#

limit for dt compared to explicit solver timescale (negative means no limit)

min_dr_div_cs = min over all cells of dr/csound (seconds)
dt_div_min_dr_div_cs_limit = -1
dt_div_min_dr_div_cs_hard_limit = -1

dt_div_dt_cell_collapse_limit#

dt_div_dt_cell_collapse_hard_limit#

limit for dt compared to cell_collapse timescale (negative means no limit)

dt_cell_collapse = min over shells k that have v(k+1) > v(k) of
(r(k)-r(k+1))/(v(k+1)-v(k)), the time for the cell to collapse
to zero thickness at current velocities.  only for v_flag true.
dt_div_dt_cell_collapse_limit = -1
dt_div_dt_cell_collapse_hard_limit = -1

min_k_for_dt_div_min_dr_div_cs_limit#

min_q_for_dt_div_min_dr_div_cs_limit#

max_q_for_dt_div_min_dr_div_cs_limit#

check_remnant_only_for_dt_div_min_dr_div_cs_limit#

min_k_for_dt_div_min_dr_div_cs_limit = 20
min_q_for_dt_div_min_dr_div_cs_limit = 0.005d0
max_q_for_dt_div_min_dr_div_cs_limit = 0.995d0
check_remnant_only_for_dt_div_min_dr_div_cs_limit = .false.

min_abs_du_div_cs_for_dt_div_min_dr_div_cs_limit#

min_abs_u_div_cs_for_dt_div_min_dr_div_cs_limit#

only use dt_div_min_dr_div_cs_limit at cells where abs_du_div_cs > min_abs_du_div_cs_for_dt_div_min_dr_div_cs_limit and and abs_u_div_cs > min_abs_u_div_cs_for_dt_div_min_dr_div_cs_limit

allow focus on regions near shock face.

min_abs_u_div_cs_for_dt_div_min_dr_div_cs_limit = 0.8d0
min_abs_du_div_cs_for_dt_div_min_dr_div_cs_limit = 0.01d0

limits based on changes in location on HR diagram

delta_HR_ds_L#

delta_HR_ds_Teff#

dlgL = log10(L/L_prev)
dlgTeff = log10(Teff/Teff_prev)
delta_HR_ds_L = 1
delta_HR_ds_Teff = 1

delta_HR_limit#

delta_HR_hard_limit#

limit for dHR (negative means no limit)

dHR = sqrt((delta_HR_ds_L*dlgL)**2 + (delta_HR_ds_Teff*dlgTeff)**2)
delta_HR_limit = -1
delta_HR_hard_limit = -1

delta_lgT_max_limit#

delta_lgT_max_hard_limit#

delta_lgT_max_limit_lgT_min#

delta_lgT_max_at_high_T_limit#

delta_lgT_max_at_high_T_hard_limit#

delta_lgT_max_at_high_T_limit_lgT_min#

delta_lgT_max_limit_only_after_near_zams#

limit for magnitude of change in max over all cells of log10 T this is for off center flashes in degenerate material (e.g., He or Ne) Only apply this limit when lgT_max >= delta_lgT_max_limit_lgT_min.

similarly, use at_high_T limits only when lgT_max >= delta_lgT_max_at_high_T_limit_lgT_min. this can be useful since higher order temperature sensitivity of rates at high T may require smaller limits for changes.

delta_lgT_max_limit = -1
delta_lgT_max_hard_limit = -1
delta_lgT_max_limit_lgT_min = 9d0
delta_lgT_max_at_high_T_limit = -1
delta_lgT_max_at_high_T_hard_limit = -1
delta_lgT_max_at_high_T_limit_lgT_min = -1
delta_lgT_max_limit_only_after_near_zams = .false.

limits based on changes at center

delta_lgT_cntr_limit#

delta_lgT_cntr_hard_limit#

delta_lgT_cntr_limit_only_after_near_zams#

limit for magnitude of change in log10 temperature at center

delta_lgT_cntr_limit = 0.01d0
delta_lgT_cntr_hard_limit = -1
delta_lgT_cntr_limit_only_after_near_zams = .true.

delta_lgP_cntr_limit#

delta_lgP_cntr_hard_limit#

limit for magnitude of change in log10 pressure at center

delta_lgP_cntr_limit = -1
delta_lgP_cntr_hard_limit = -1

delta_lgRho_cntr_limit#

delta_lgRho_cntr_hard_limit#

limit for magnitude of change in log10 density at center

delta_lgRho_cntr_limit = 0.05d0
delta_lgRho_cntr_hard_limit = -1

dX_div_X_cntr is max(abs(xa(j,nz)-xa_old(j,nz))/xa(j,nz)) for any species j Small timesteps as the center carbon is exhausted.

delta_dX_div_X_cntr_min#

Ignore changes in dX_div_X_cntr if value is less than this.

delta_dX_div_X_cntr_min = -5

delta_dX_div_X_cntr_max#

Ignore changes in dX_div_X_cntr if value is more than this.

delta_dX_div_X_cntr_max = 0

delta_dX_div_X_cntr_limit#

If max delta is greater than delta_dX_div_X_cntr_limit, reduce the next timestep by delta_dX_div_X_cntr_limit/max_delta.

delta_dX_div_X_cntr_limit = 0.1d0

delta_dX_div_X_cntr_hard_limit#

If max delta is greater than delta_dX_div_X_cntr_hard_limit, retry with smaller timestep.

delta_dX_div_X_cntr_hard_limit = -1

lg_XH_cntr is log10(h1 mass fraction at center). Small timesteps as the center hydrogen is exhausted.

delta_lg_XH_cntr_min#

Ignore changes in lg_XH_cntr if value is less than this.

delta_lg_XH_cntr_min = -6

delta_lg_XH_cntr_max#

Ignore changes in lg_XH_cntr if value is more than this.

delta_lg_XH_cntr_max = 0

delta_lg_XH_cntr_limit#

If max delta is greater than this, reduce the next timestep by delta_lg_XH_cntr_limit/max_delta.

delta_lg_XH_cntr_limit = 0.05d0

delta_lg_XH_cntr_hard_limit#

If max delta is greater than delta_lg_XH_cntr_hard_limit, retry with smaller timestep.

delta_lg_XH_cntr_hard_limit = -1

lg_XHe_cntr is log10(he4 mass fraction at center) small timesteps as the center helium is exhausted.

delta_lg_XHe_cntr_min#

Ignore changes in lg_XHe_cntr if value is less than this.

delta_lg_XHe_cntr_min = -6

delta_lg_XHe_cntr_max#

Ignore changes in lg_XHe_cntr if value is more than this.

delta_lg_XHe_cntr_max = 0

delta_lg_XHe_cntr_limit#

If max delta is greater than delta_lg_XHe_cntr_limit, reduce the next timestep by delta_lg_XHe_cntr_limit/max_delta.

delta_lg_XHe_cntr_limit = 0.1d0

delta_lg_XHe_cntr_hard_limit#

If max delta is greater than delta_lg_XHe_cntr_hard_limit, retry with smaller timestep.

delta_lg_XHe_cntr_hard_limit = -1

lg_XC_cntr is log10(c12 mass fraction at center). Small timesteps as the center carbon is exhausted.

delta_lg_XC_cntr_min#

Ignore changes in lg_XC_cntr if value is less than this.

delta_lg_XC_cntr_min = -5

delta_lg_XC_cntr_max#

Ignore changes in lg_XC_cntr if value is more than this.

delta_lg_XC_cntr_max = 0

delta_lg_XC_cntr_limit#

If max delta is greater than delta_lg_XC_cntr_limit, reduce the next timestep by delta_lg_XC_cntr_limit/max_delta.

delta_lg_XC_cntr_limit = 0.1d0

delta_lg_XC_cntr_hard_limit#

If max delta is greater than delta_lg_XC_cntr_hard_limit, retry with smaller timestep.

delta_lg_XC_cntr_hard_limit = -1

lg_XO_cntr is log10(o16 mass fraction at center) Small timesteps as the center oxygen is exhausted.

delta_lg_XNe_cntr_min#

Ignore changes in lg_XNe_cntr if value is less than this.

delta_lg_XNe_cntr_min = -5

delta_lg_XNe_cntr_max#

Ignore changes in lg_XNe_cntr if value is more than this.

delta_lg_XNe_cntr_max = 0

delta_lg_XNe_cntr_limit#

If max delta is greater than delta_lg_XNe_cntr_limit, reduce the next timestep by delta_lg_XNe_cntr_limit/max_delta.

delta_lg_XNe_cntr_limit = 1d99

delta_lg_XNe_cntr_hard_limit#

If max delta is greater than delta_lg_XNe_cntr_hard_limit, retry with smaller timestep.

delta_lg_XNe_cntr_hard_limit = -1

delta_lg_XO_cntr_min#

Ignore changes in lg_XO_cntr if value is less than this.

delta_lg_XO_cntr_min = -5

delta_lg_XO_cntr_max#

Ignore changes in lg_XO_cntr if value is more than this.

delta_lg_XO_cntr_max = 0

delta_lg_XO_cntr_limit#

If max delta is greater than delta_lg_XO_cntr_limit, reduce the next timestep by delta_lg_XO_cntr_limit/max_delta.

delta_lg_XO_cntr_limit = 1d99

delta_lg_XO_cntr_hard_limit#

If max delta is greater than delta_lg_XO_cntr_hard_limit, retry with smaller timestep.

delta_lg_XO_cntr_hard_limit = -1

delta_lg_XSi_cntr_min#

Ignore changes in lg_XSi_cntr if value is less than this.

delta_lg_XSi_cntr_min = -5

delta_lg_XSi_cntr_max#

Ignore changes in lg_XSi_cntr if value is more than this.

delta_lg_XSi_cntr_max = 0

delta_lg_XSi_cntr_limit#

If max delta is greater than delta_lg_XSi_cntr_limit, reduce the next timestep by delta_lg_XSi_cntr_limit/max_delta.

delta_lg_XSi_cntr_limit = 1d99

delta_lg_XSi_cntr_hard_limit#

If max delta is greater than delta_lg_XSi_cntr_hard_limit, retry with smaller timestep.

delta_lg_XSi_cntr_hard_limit = -1

delta_XH_cntr_limit#

If max delta is greater than this, reduce the next timestep by delta_XH_cntr_limit/max_delta.

delta_XH_cntr_limit = 0.01d0

delta_XH_cntr_hard_limit#

If max delta is greater than delta_XH_cntr_hard_limit, retry with smaller timestep.

delta_XH_cntr_hard_limit = -1

delta_XHe_cntr_limit#

If max delta is greater than delta_XHe_cntr_limit, reduce the next timestep by delta_XHe_cntr_limit/max_delta.

delta_XHe_cntr_limit = 0.01d0

delta_XHe_cntr_hard_limit#

If max delta is greater than delta_XHe_cntr_hard_limit, retry with smaller timestep.

delta_XHe_cntr_hard_limit = -1

delta_XC_cntr_limit#

If max delta is greater than delta_XC_cntr_limit, reduce the next timestep by delta_XC_cntr_limit/max_delta.

delta_XC_cntr_limit = 0.01d0

delta_XC_cntr_hard_limit#

If max delta is greater than delta_XC_cntr_hard_limit, retry with smaller timestep.

delta_XC_cntr_hard_limit = -1

delta_XNe_cntr_limit#

If max delta is greater than delta_XNe_cntr_limit, reduce the next timestep by delta_XNe_cntr_limit/max_delta.

delta_XNe_cntr_limit = 0.01d0

delta_XNe_cntr_hard_limit#

If max delta is greater than delta_XNe_cntr_hard_limit, retry with smaller timestep.

delta_XNe_cntr_hard_limit = -1

delta_XO_cntr_limit#

If max delta is greater than delta_XO_cntr_limit, reduce the next timestep by delta_XO_cntr_limit/max_delta.

delta_XO_cntr_limit = 0.01d0

delta_XO_cntr_hard_limit#

If max delta is greater than delta_XO_cntr_hard_limit, retry with smaller timestep.

delta_XO_cntr_hard_limit = -1

delta_XSi_cntr_limit#

If max delta is greater than delta_XSi_cntr_limit, reduce the next timestep by delta_XSi_cntr_limit/max_delta.

delta_XSi_cntr_limit = 0.01d0

delta_XSi_cntr_hard_limit#

If max delta is greater than delta_XSi_cntr_hard_limit, retry with smaller timestep.

delta_XSi_cntr_hard_limit = -1

delta_dX_div_X_drop_only#

delta_lg_XH_drop_only#

delta_lg_XHe_drop_only#

delta_lg_XC_drop_only#

delta_lg_XNe_drop_only#

delta_lg_XO_drop_only#

delta_lg_XSi_drop_only#

delta_XH_drop_only#

delta_XHe_drop_only#

delta_XC_drop_only#

delta_XNe_drop_only#

delta_XO_drop_only#

delta_XSi_drop_only#

If true, then only limit drops in abundance.

delta_dX_div_X_drop_only = .false.
delta_lg_XH_drop_only = .false.
delta_lg_XHe_drop_only = .false.
delta_lg_XC_drop_only = .false.
delta_lg_XNe_drop_only = .false.
delta_lg_XO_drop_only = .false.
delta_lg_XSi_drop_only = .false.
delta_XH_drop_only = .false.
delta_XHe_drop_only = .false.
delta_XC_drop_only = .false.
delta_XNe_drop_only = .false.
delta_XO_drop_only = .false.
delta_XSi_drop_only = .false.

limits based on changes in mass of the star#

delta_lg_star_mass_limit#

delta_lg_star_mass_hard_limit#

Limit for magnitude of change in log10(M/Msun).

delta_lg_star_mass_limit = 5d-3
delta_lg_star_mass_hard_limit = -1

limit for change in mdot in Msun/yr + delta_mdot_atol tolerance for absolute changes + delta_mdot_rtol tolerance for relative changes

delta_mdot_atol = 1d-3
delta_mdot_rtol = 0.5d0

delta_mdot_limit#

delta_mdot_hard_limit#

delta_mot = abs(mdot - mdot_old)/ (delta_mdot_atol*Msun/secyer + &
   delta_mdot_rtol*max(abs(mdot),abs(mdot_old)))

ignore if < 0

delta_mdot_limit = -1
delta_mdot_hard_limit = -1

adjust_J_q_limit#

adjust_J_q_hard_limit#

limit for mass coordinate down to which angular momentum is adjusted when using do_adjust_J_lost

adjust_J_q_limit = 0.99
adjust_J_q_hard_limit = 0.98

limit_for_rel_error_in_energy_conservation#

hard_limit_for_rel_error_in_energy_conservation#

rel_error_in_energy_conservation = abs(error_in_energy_conservation/total_energy)
limit_for_rel_error_in_energy_conservation = 1d-4
hard_limit_for_rel_error_in_energy_conservation = 1d-3

report_min_dr_div_cs#

If true, produce terminal output about minimum of cell dr/cs

report_min_dr_div_cs = .false.

report_dt_hard_limit_retries#

If true, produce terminal output about choice of timestep.

report_dt_hard_limit_retries = .false.

report_solver_dt_info#

If true, produce terminal output about choice of timestep based on varcontrol_target.

report_solver_dt_info = .false.

debugging controls#

report_solver_progress#

Set true to see info about solver iterations.

report_solver_progress = .false.

report_ierr#

If true, produce terminal output when have some internal error.

report_ierr = .false.

report_bad_negative_xa#

If true, produce terminal output when have bad negative mass fraction error.

report_bad_negative_xa = .false.

stop_for_bad_nums#

If true, then stop for bad numbers (NaNs or infinity). this replaces old control stop_for_NaNs

stop_for_bad_nums = .false.

show_mesh_changes#

When show_mesh_changes is true, the terminal output includes the mesh_call_number.

show_mesh_changes = .false.

trace_evolve#

Send evolve output to screen.

trace_evolve = .false.

variety of output from the solver

solver_numerical_jacobian = .false.
solver_jacobian_nzlo = 1
solver_jacobian_nzhi = -1
solver_check_everything = .false.
solver_inspect_soln_flag = .false.
solver_test_partials_dx_0 = -1d0
solver_test_partials_k = -1
solver_test_partials_k_low = -1
solver_test_partials_k_high = -1
solver_test_eos_partials = .false.
solver_test_kap_partials = .false.
solver_test_net_partials = .false.
solver_test_atm_partials = .false.
solver_test_partials_var_name = ''
solver_test_partials_sink_name = ''
solver_test_partials_equ_name = ''
solver_test_partials_show_dx_var_name = ''
solver_show_correction_info = .false.
solver_test_partials_call_number = -1
solver_test_partials_iter_number = -1
solver_test_partials_write_eos_call_info = .false.
solver_epsder_struct = 1d-5
solver_epsder_chem = 1d-5

energy_conservation_dump_model_number = -1

solver_save_photo_call_number#

Saves a photo when solver_call_number = solver_save_photo_call_number - 1

solver_save_photo_call_number = -1

xa_clip_limit#

Abundances smaller than this limit are set to 0.

xa_clip_limit = 1d-99

fill_arrays_with_NaNs#

initialize arrays with NaNs to trap reads of uninitialized entries.

fill_arrays_with_NaNs = .false.

zero_when_allocate#

initialize arrays with zeros.

zero_when_allocate = .false.

miscellaneous controls#

warn_rates_for_high_temp#

If true then when any zone tries to evaluate a rate above max_safe_logT_for_rates it generates a warning message. The code will cap the rate at the value for max_safe_logT_for_rates whether the warning is on or not. 10d0 is a sensible default for the max temperature, as that is where the partition tables and polynomial fits to the rates are valid until. warning messages include the text “rates have been truncated” and “WARNING: evaluating rates”.

warn_rates_for_high_temp = .true.
max_safe_logT_for_rates = 10d0

warn_when_large_rel_run_E_err#

message includes the text “WARNING: rel_run_E_err”

rel_run_E_err = abs(cumulative_energy_error/total_energy) you can turn off this warning message by setting this to a large number.

warn_when_large_rel_run_E_err = 0.1d0

absolute_cumulative_energy_err#

cumulative energy error is the sum of the absolute value of per-step errors. Set this to .false. to allow positive and negative errors to cancel when integrating over multiple steps.

absolute_cumulative_energy_err = .true.

warn_when_stop_checking_residuals obsolete#

warning_limit_for_max_residual#

message includes the text “WARNING: max_residual > warning_limit_for_max_residual”

warning_limit_for_max_residual = 1d0

warn_when_large_virial_thm_rel_err#

message includes the text “WARNING: virial_thm_rel_err” only applies to models with no velocities, rotation, or mass change.

warn_when_large_virial_thm_rel_err = 1d-2

warn_when_get_a_bad_eos_result#

message includes the text “WARNING eos:”

warn_when_get_a_bad_eos_result = .true.

relax_dY#

Change Y by this amount per step when relaxing Y.

relax_dY = 0.005d0

relax_dlnZ#

Change lnZ by this amount per step when relaxing Z. Default is ln10/10.

relax_dlnZ = 2.3025850929940459d-1

eps_mdot_leak_frac_factor#

eps_mdot_leak_frac_factor = 1d0

zams_filename#

Default is for Z=0.02, Y=0.28.

zams_filename = 'zams_z2m2_y28.data'
set_rho_to_dm_div_dV = .true.

steps_before_start_stress_test#

stress_test_relax#

steps_before_start_stress_test = 0
stress_test_relax = .false.

use_other_{hook}#

Logicals to deploy the use_other routines.

use_other_kap = .false.
use_other_diffusion = .false.
use_other_mlt_results = .false.
use_other_adjust_mdot = .false.
use_other_j_for_adjust_J_lost = .false.
use_other_alpha_mlt = .false.
use_other_am_mixing = .false.
use_other_brunt = .false.
use_other_brunt_smoothing = .false.
use_other_build_initial_model = .false.
use_other_cgrav = .false.
use_other_mesh_delta_coeff_factor = .false.
use_other_energy_implicit = .false.
use_other_energy = .false.
use_other_pressure = .false.
use_other_momentum_implicit = .false.
use_other_momentum = .false.
use_other_eps_grav = .false.
use_other_mesh_functions = .false.
use_other_D_mix = .false.
use_other_close_gaps = .false.
use_other_neu = .false.
use_other_net_get = .false.
use_other_solver_monitor = .false.
use_other_opacity_factor = .false.
use_other_diffusion_factor = .false.
use_other_diffusion_coefficients = .false.
use_other_pgstar_plots = .false.
use_other_gradr_factor = .false.
use_other_eval_fp_ft = .false.
use_other_torque = .false.
use_other_screening = .false.
use_other_rate_get = .false.
use_other_net_derivs = .false.
use_other_split_burn = .false.
use_other_torque_implicit = .false.
use_other_wind = .false.
use_other_after_struct_burn_mix = .false.
use_other_before_struct_burn_mix = .false.
use_other_surface_PT = .false.
use_other_remove_surface = .false.
use_other_set_pgstar_controls = .false.
use_other_accreting_state = .false.
use_other_eval_i_rot = .false.
use_other_export_pulse_data = .false.
use_other_get_pulse_data = .false.
use_other_edit_pulse_data = .false.
use_other_astero_freq_corr = .false.
use_other_timestep_limit = .false.

mixing diffusion coeffs#

sig_term_limit#

Limit on coefficients in convective mixing equations. Consider a diffusion eqn of form:

x(k) - x0(k) = c1*(x(k-1) - x(k)) - c2*(x(k) - x(k+1))

Simplify for c1=c2=c, x(k-1)=x(k+1)=x0(k)=x0, x(k)=x0+dx Then eqn becomes

(1+2*c)*(x0+dx) - 2*c*x0 = x0

If 2*c >> 1, then eqn becomes ill-conditioned, so we enforce c <= sig_term_limit In physical terms c is dt*sig/dm, where sig = (4 pi r^2 rho)^2*D and D = diffusion coeff (cm^2/s), so c can get large when dt/dm is large.

sig_term_limit = 1d13

am_sig_term_limit#

Limit on coefficients in angular momentum transport equations. Necessary for numerical stability. Plays same role as sig_term_limit for material mixing.

am_sig_term_limit = 1d13

sig_min_factor_for_high_Tcenter#

High center T limit to avoid negative mass fractions. If Tcenter >= Tcenter_min_for_sig_min_factor_full_on, then okay to reduce sig by as much as this factor as needed to prevent causing negative abundances. Inactive when >= 1d0.

sig_min_factor_for_high_Tcenter = 0.01d0

Tcenter_min_for_sig_min_factor_full_on#

If Tcenter >= this, factor = sig_min_factor_for_neg_abundances, this should be > Tcenter_max_for_sig_min_factor_full_off.

Tcenter_min_for_sig_min_factor_full_on = 3.2d9

Tcenter_max_for_sig_min_factor_full_off#

If Tcenter <= this, factor = 1, so has no effect this should be < Tcenter_min_for_sig_min_factor_full_on. For T > full_off and < full_on, factor changes linearly with Tcenter.

Tcenter_max_for_sig_min_factor_full_off = 2.8d9

max_delta_m_to_bdy_for_sig_min_factor#

sig_min factor goes to 1 as distance (in Msun units) from boundary of mixing region reaches this value

max_delta_m_to_bdy_for_sig_min_factor = 0.5d0

delta_m_upper_for_sig_min_factor#

okay to change sig min factor to 1 for mix region larger than this

delta_m_upper_for_sig_min_factor = 0.3d0

delta_m_lower_for_sig_min_factor#

don’t change sig min factor for mix region smaller than this

delta_m_lower_for_sig_min_factor = 0.1d0

extra params as a convenience for developing new features note: the parameter num_x_ctrls is defined in star_def.inc

x_ctrl(1:num_x_ctrls) = 0d0
x_integer_ctrl(1:num_x_ctrls) = 0
x_logical_ctrl(1:num_x_ctrls) = .false.
x_character_ctrl(1:num_x_ctrls) = ''

One can split controls inlist into pieces using the following parameters. BTW: it works recursively, so the extras can read extras too.

read_extra_controls_inlist(1..5)#

extra_controls_inlist_name(1..5)#

If read_extra_controls_inlist(1) is true, then read &controls from this namelist file.

read_extra_controls_inlist(:) = .false.
extra_controls_inlist_name(:) = 'undefined'

save_controls_namelist#

dumps all values for &controls namelist to file

save_controls_namelist = .false.

controls_namelist_name#

if empty, uses a default name

controls_namelist_name = ''