NuclearBurning.jl
This notebook provides a simple example of a star with simplified microphysics undergoing nuclear burning. Import all necessary Jems modules. We will also do some benchmarks, so we import BenchmarkTools as well.
using BenchmarkTools
using Jems.Chem
using Jems.Constants
using Jems.EOS
using Jems.Opacity
using Jems.NuclearNetworks
using Jems.Turbulence
using Jems.StellarModels
using Jems.Evolution
using Jems.ReactionRatesModel creation
We start by creating the stellar model. In this example we consider a model with 6 independent variables, two of which correspond to composition. The independent variables here are $\ln(P)$, $\ln(T)$, $\ln(r)$, the luminosity $L$ and the mass fractions of Hydrogen and Helium.
The Evolution module has pre-defined equations corresponding to these variables, which we provide here. For now, only a simple (fully ionized) ideal gas law EOS is available. Similarly, only a simple simple electron scattering opacity equal to $\kappa=0.2(1+X)\;[\mathrm{cm^2\;g^{-1}}]$ is available.
varnames = [:lnρ, :lnT, :lnr, :lum]
varscaling = [:log, :log, :log, :maxval]
structure_equations = [Evolution.equationHSE, Evolution.equationT,
Evolution.equationContinuity, Evolution.equationLuminosity]
remesh_split_functions = [StellarModels.split_lnr_lnρ, StellarModels.split_lum,
StellarModels.split_lnT, StellarModels.split_xa]
net = NuclearNetwork([:H1, :He4, :C12, :N14, :O16], [(:kipp_rates, :kipp_pp), (:kipp_rates, :kipp_cno)])
nz = 1000
nextra = 100
eos = EOS.IdealEOS(true)
opacity = Opacity.SimpleElectronScatteringOpacity()
turbulence = Turbulence.BasicMLT(1.0)
sm = StellarModel(varnames, varscaling, structure_equations, Evolution.equation_composition,
nz, nextra, remesh_split_functions, net, eos, opacity, turbulence);Initialize StellarModel and evaluate equations and jacobian
We do not have a working initial condition yet. We require pressure, temperature profiles. One simple available initial condition is that of an n=1 polytrope. This sets the pressure and density and computes the temperature from the EOS. The luminosity is initialized by assuming pure radiative transport for the temperature gradient produced by the polytrope.
The normal evolution loop will store the information at the end of the step into an attribute of type StellarStepInfo, stored at sm.esi (end step info). After initializing our polytrope we can mimic that behavior by calling set_end_step_info!(sm). We then 'cycle' this info into the information of a hypothetical previous step with cycle_step_info, so now sm.psi contains our initial condition. Finally we call set_start_step_info to use sm.psi (previous step info) to populate the information needed before the Newton solver in sm.ssi (start step info). At last we are in position to evaluate the equations and compute the Jacobian.
n = 3
StellarModels.n_polytrope_initial_condition!(n, sm, nz, 0.7154, 0.0142, 0.0, Chem.abundance_lists[:ASG_09], MSUN,
100 * RSUN; initial_dt=10 * SECYEAR)
StellarModels.evaluate_stellar_model_properties!(sm, sm.props)
StellarModels.cycle_props!(sm);
StellarModels.copy_scalar_properties!(sm.start_step_props, sm.prv_step_props)
StellarModels.copy_mesh_properties!(sm, sm.start_step_props, sm.prv_step_props) # or do StellarModels.remesher!(sm);
StellarModels.evaluate_stellar_model_properties!(sm, sm.start_step_props)
StellarModels.copy_scalar_properties!(sm.props, sm.start_step_props)
StellarModels.copy_mesh_properties!(sm, sm.props, sm.start_step_props)Benchmarking
The previous code leaves everything ready to solve the linearized system. For now we make use of a the serial Thomas algorithm for tridiagonal block matrices. We first show how long it takes to evaluate the Jacobian matrix. This requires two steps, the first is to evaluate properties across the model (for example, the EOS) and then evaluate all differential equations and fill the Jacobian.
@benchmark begin
StellarModels.evaluate_stellar_model_properties!($sm, $sm.props)
end
@benchmark begin
Evolution.eval_jacobian_eqs!($sm)
endBenchmarkTools.Trial: 2339 samples with 1 evaluation per sample.
Range (min … max): 1.919 ms … 4.688 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.058 ms ┊ GC (median): 0.00%
Time (mean ± σ): 2.134 ms ± 227.621 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▄▇█▅▂▃
▅████████▇▇▇▇▅▆▅▅▄▄▄▄▃▄▃▃▃▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂ ▃
1.92 ms Histogram: frequency by time 3.04 ms <
Memory estimate: 592 bytes, allocs estimate: 7.To get an idea of how much a complete iteration of the solver takes, we need to Perform the jacobian evaluation as a setup for the benchmark. This is because the solver destroys the Jacobian to perform in-place operations.
@benchmark begin
Evolution.thomas_algorithm!($sm)
end setup=(Evolution.eval_jacobian_eqs!($sm))BenchmarkTools.Trial: 1203 samples with 1 evaluation per sample.
Range (min … max): 1.678 ms … 2.596 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.731 ms ┊ GC (median): 0.00%
Time (mean ± σ): 1.737 ms ± 40.495 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▇█▃▃▁
▂▂▂▂▃▃▃▃▂▂▂▂▂▂▄▅▆▇██████▇▅▅▄▄▄▃▃▃▃▃▂▂▂▃▂▃▂▂▂▂▂▂▂▂▂▁▁▁▂▂▁▂▂ ▃
1.68 ms Histogram: frequency by time 1.82 ms <
Memory estimate: 125.00 KiB, allocs estimate: 2000.Evolving our model
We can now evolve our star! We will initiate a $1M_\odot$ star with a radius of $100R_\odot$ using an n=1 polytrope (it would be much better to use n=3 or n=3/2 polytropes, for now I only use this because there is a simple analytical solution). The star is expected to contract until it ignites hydrogen. We set a few options for the simulation with a toml file, which we generate dynamically. These simulation should complete in about a thousand steps once it reaches the max_center_T limit.
Output is stored in HDF5 files, and easy to use functions are provided with the Evolution module to turn these HDF5 files into DataFrame objects. HDF5 output is compressed by default.
open("example_options.toml", "w") do file
write(file,
"""
[remesh]
do_remesh = true
[solver]
newton_max_iter_first_step = 1000
initial_model_scale_max_correction = 0.2
newton_max_iter = 50
scale_max_correction = 0.1
[timestep]
dt_max_increase = 1.5
delta_R_limit = 0.01
delta_Tc_limit = 0.01
delta_Xc_limit = 0.005
[termination]
max_model_number = 2000
max_center_T = 1e8
[plotting]
do_plotting = true
wait_at_termination = false
plotting_interval = 1
window_specs = ["HR", "Kippenhahn", "profile", "TRhoProfile"]
window_layout = [[1, 1], # arrangement of plots
[1, 2],
[2, 1],
[2, 2]
]
profile_xaxis = 'mass'
profile_yaxes = ['log10_T']
profile_alt_yaxes = ['X','Y']
history_xaxis = 'age'
history_yaxes = ['R_surf']
history_alt_yaxes = ['T_center']
max_log_eps = 5.0
[io]
profile_interval = 50
terminal_header_interval = 100
terminal_info_interval = 100
""")
end
StellarModels.set_options!(sm.opt, "./example_options.toml")
rm(sm.opt.io.hdf5_history_filename; force=true)
rm(sm.opt.io.hdf5_profile_filename; force=true)
n = 3
StellarModels.n_polytrope_initial_condition!(n, sm, nz, 0.7154, 0.0142, 0.0, Chem.abundance_lists[:ASG_09],
1 * MSUN, 100 * RSUN; initial_dt=10 * SECYEAR)
@time Evolution.do_evolution_loop!(sm);Found first model
model logdt logL logTeff logPs logρs H_cntr iters
mass age logR logTc logPc logρc He_cntr zones
-------------------------------------------------------------------------------
1 1.000 1.205 3.069 0.7508 -10.45 0.7154 23
1.000 10.00 1.988 5.080 9.096 -4.116 0.2704 1004
model logdt logL logTeff logPs logρs H_cntr iters
mass age logR logTc logPc logρc He_cntr zones
-------------------------------------------------------------------------------
100 2.617 1.103 3.243 1.549 -9.826 0.7154 4
1.000 2.361e+04 1.589 5.462 10.54 -3.049 0.2704 1004
model logdt logL logTeff logPs logρs H_cntr iters
mass age logR logTc logPc logρc He_cntr zones
-------------------------------------------------------------------------------
200 3.049 1.102 3.456 2.404 -9.184 0.7154 4
1.000 9.494e+04 1.161 5.892 12.26 -1.760 0.2704 1004
model logdt logL logTeff logPs logρs H_cntr iters
mass age logR logTc logPc logρc He_cntr zones
-------------------------------------------------------------------------------
300 3.477 1.102 3.670 3.260 -8.542 0.7154 4
1.000 2.866e+05 0.7331 6.321 13.98 -0.4712 0.2704 1009
model logdt logL logTeff logPs logρs H_cntr iters
mass age logR logTc logPc logρc He_cntr zones
-------------------------------------------------------------------------------
400 3.906 1.102 3.884 4.116 -7.901 0.7154 4
1.000 8.008e+05 0.3052 6.750 15.70 0.8152 0.2704 1019
model logdt logL logTeff logPs logρs H_cntr iters
mass age logR logTc logPc logρc He_cntr zones
-------------------------------------------------------------------------------
500 4.374 1.089 4.095 4.971 -7.256 0.7154 4
1.000 2.204e+06 -0.1223 7.181 17.38 2.064 0.2704 1031
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (535, 50, 3.883769839288261e13, -0.0021559253489612856, 358, 6, 2162.206351543362, 358, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (537, 50, 4.438967367461355e13, -0.002464121888457657, 357, 6, 2471.3178172787666, 357, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (539, 50, 5.083198628338282e13, -0.0028217420779620872, 356, 6, 2830.003529559306, 356, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (542, 50, 4.301682902705433e13, 0.004775827403219949, 354, 5, 4800.695930114982, 352, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (544, 50, 4.311344142345696e13, 0.00478655353391499, 352, 5, 4811.677025679843, 350, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (545, 50, 6.509314411934924e13, -0.003613395367327257, 351, 6, 3624.108670579656, 351, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (547, 50, 7.503935883018642e13, -0.004165521196319101, 350, 6, 4187.57519650951, 348, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (548, 50, 5.6679028170934305e13, -0.003146318104437457, 349, 6, 3155.696463513475, 349, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (551, 50, 4.325624748508586e13, 0.004802408191682634, 347, 5, 4828.147660656283, 345, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (552, 50, 3.2654015695736395e13, 0.0036253240071877795, 346, 5, 3636.213135619884, 346, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (555, 50, 8.683039309248997e13, -0.004820055082426776, 344, 6, 4846.242638997615, 342, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (556, 50, 3.272916260304412e13, 0.0036336669898600206, 343, 5, 3644.670947823182, 343, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (559, 50, 4.350859041679295e13, 0.004830423884970897, 341, 5, 4857.006001094588, 339, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (560, 50, 6.559972285548941e13, -0.003641516136159673, 340, 6, 3652.6323089405687, 340, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (564, 50, 4.379473737063507e13, 0.004862192578629084, 337, 5, 4889.4272331449165, 335, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (565, 50, 6.5936149783107586e13, -0.003660191582216953, 336, 6, 3671.4877796567284, 336, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (571, 50, 8.83127793413039e13, -0.004902344049667823, 330, 6, 4917.704171149259, 330, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (572, 50, 3.3154855306772742e13, 0.0036809283739694656, 330, 5, 3692.479506698632, 330, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (585, 50, 8.881606428247361e13, -0.004930281976149509, 318, 6, 4946.287727960346, 318, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (589, 50, 8.898246974032345e13, -0.004939519334685919, 315, 6, 4955.696779215043, 315, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (592, 50, 8.892978783448883e13, -0.004936594901443982, 313, 6, 4952.859733638613, 313, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (595, 50, 8.903764994862914e13, -0.004942582451574148, 311, 6, 4958.963201202845, 311, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (598, 50, 8.914021192487538e13, -0.004948275784948129, 309, 6, 4964.7709519815135, 309, 5)
model logdt logL logTeff logPs logρs H_cntr iters
mass age logR logTc logPc logρc He_cntr zones
-------------------------------------------------------------------------------
600 6.174 1.025 4.108 5.088 -7.152 0.4857 4
1.000 1.461e+08 -0.1808 7.310 17.57 2.206 0.5001 1033
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (600, 50, 4.4474319750939984e13, -0.004937641379203493, 308, 6, 4954.154103124135, 308, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (602, 50, 8.900425384178662e13, 0.004940728595236678, 307, 5, 4957.304426862038, 307, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (605, 50, 8.93434238738333e13, 0.004959556314178731, 305, 5, 4976.289410861738, 305, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (607, 50, 8.915512612034461e13, 0.0049491036891083145, 304, 5, 4965.854059362167, 304, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (609, 50, 8.92044587400021e13, 0.004951842199618755, 303, 5, 4968.6541612084675, 303, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (611, 50, 8.925450477686486e13, 0.004954620312739569, 302, 5, 4971.493912643433, 302, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (613, 50, 8.930301691077208e13, 0.004957313277141444, 301, 5, 4974.248137733503, 301, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (615, 50, 8.935001904477452e13, 0.004959922419710245, 300, 5, 4976.918163366888, 300, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (617, 50, 8.939555259874575e13, 0.004962450039710839, 299, 5, 4979.506291322965, 299, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (619, 50, 8.943965765824828e13, 0.004964898362338948, 298, 5, 4982.01474831688, 298, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (621, 50, 8.948237278836447e13, 0.00496726952838669, 297, 5, 4984.445675513768, 297, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (623, 50, 4.476186752934859e13, -0.004969565595630398, 296, 6, 4986.801129989137, 296, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (625, 50, 4.478189003470768e13, -0.004971788540276432, 295, 6, 4989.083086135971, 295, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (626, 50, 6.737271585341871e13, 0.003739937018599688, 295, 5, 3752.9636078934664, 295, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (628, 50, 4.832114954036136e12, -0.004291780223685158, 294, 1, 4298.663481037591, 294, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (629, 50, 3.6394456770131606e12, -0.0032324771099940014, 294, 1, 3234.267156531768, 294, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (634, 50, 5.756994555816382e12, -0.005113238326807218, 291, 1, 5017.454473722435, 291, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (636, 50, 5.775903677150523e12, -0.005130033000311694, 290, 1, 5000.674281640743, 290, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (637, 50, 4.350066099149867e12, -0.0038636348335340173, 290, 1, 3761.3479434608757, 290, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (642, 50, 5.934516903338314e12, -0.005270909844890705, 287, 1, 5025.216104717741, 287, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (644, 50, 5.956267008384105e12, -0.005290227818818587, 286, 1, 5009.0271073209715, 286, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (649, 50, 6.121455845551771e12, -0.005436944979166266, 283, 1, 5031.752801664377, 283, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (654, 50, 6.27250841346462e12, -0.005571106610226746, 280, 1, 5036.215212733044, 280, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (659, 50, 6.431027032018131e12, -0.00571189942634665, 277, 1, 5040.125739694274, 277, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (663, 50, 6.551786949357047e12, -0.005819155778891846, 275, 1, 5043.071373025258, 275, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (665, 50, 6.584576525139198e12, -0.005848278772492675, 274, 1, 5028.948761037793, 274, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (681, 50, 7.26153908449943e12, -0.006449542308661398, 264, 1, 5053.357105966197, 264, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (683, 50, 7.308133182239589e12, -0.006490926180759605, 263, 1, 5040.884617861525, 263, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (687, 50, 7.481872087072307e12, -0.006645237326694359, 261, 1, 5054.789022747512, 261, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (691, 50, 7.643487131553261e12, -0.006788780321501218, 259, 1, 5055.756791958109, 259, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (693, 50, 7.697853106097898e12, -0.006837067006857732, 258, 1, 5043.9926466139095, 258, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (694, 50, 5.79599609015837e12, -0.005147878647944966, 258, 1, 3793.2586333971585, 258, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (698, 50, 7.883795412449002e12, -0.007002216950668052, 256, 1, 5053.657808998066, 256, 6)
model logdt logL logTeff logPs logρs H_cntr iters
mass age logR logTc logPc logρc He_cntr zones
-------------------------------------------------------------------------------
700 5.672 1.135 4.074 4.844 -7.363 0.1092 4
1.000 2.311e+08 -0.05878 7.369 17.54 2.301 0.8766 1033
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (702, 50, 8.063176505838143e12, -0.007161539366718499, 254, 1, 5055.583799338941, 253, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (709, 50, 8.4150807182658955e12, -0.007474093093998398, 250, 1, 5050.67314515272, 249, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (711, 50, 8.477321527582129e12, -0.007529374037657747, 249, 1, 5038.1547412477175, 248, 6)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (716, 50, 8.708237820839814e12, -0.007734468906086369, 246, 1, 5029.980224322171, 245, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (718, 50, 8.736383779588859e12, -0.00775946753924904, 245, 1, 5009.285089159549, 244, 5)
(sm.props.model_number, i, rel_corr, signed_max_corr, corr_nz, corr_equ, max_res, res_nz, res_equ) = (722, 50, 8.691829447502206e12, 0.007719895342985522, 241, 1, 0.004916445299819093, 241, 5)
model logdt logL logTeff logPs logρs H_cntr iters
mass age logR logTc logPc logρc He_cntr zones
-------------------------------------------------------------------------------
800 4.907 1.319 3.975 4.264 -7.844 6.543e-175 4
1.000 3.021e+08 0.2314 7.712 20.25 4.749 0.9858 1033
model logdt logL logTeff logPs logρs H_cntr iters
mass age logR logTc logPc logρc He_cntr zones
-------------------------------------------------------------------------------
900 4.285 1.360 3.764 3.378 -8.519 -1.018e-321 4
1.000 3.059e+08 0.6742 7.950 21.03 5.285 0.9858 1047
model logdt logL logTeff logPs logρs H_cntr iters
mass age logR logTc logPc logρc He_cntr zones
-------------------------------------------------------------------------------
962 4.040 1.376 3.632 2.832 -8.932 -2.520e-322 5
1.000 3.068e+08 0.9473 8.000 21.20 5.402 0.9858 1061
Reached maximum central temperature
632.021013 seconds (455.19 M allocations: 36.282 GiB, 1.00% gc time, 2.87% compilation time: 3% of which was recompilation)Plotting with Makie
Now that our simulation is complete we can analyze the results. We make use of the Makie package for this. I'm not a fan of the Makie defaults, so I adjust them. I normally also adjust the fonts to be consistent with \LaTeX, but I avoid that here so we don't need to distribute those fonts together with Jems.
using CairoMakie, LaTeXStrings, MathTeXEngine
basic_theme = Theme(fonts=(regular=texfont(:text), bold=texfont(:bold),
italic=texfont(:italic), bold_italic=texfont(:bolditalic)),
fontsize=30, size=(1000, 750), linewidth=7,
Axis=(xlabelsize=40, ylabelsize=40, titlesize=40, xgridvisible=false, ygridvisible=false,
spinewidth=2.5, xminorticksvisible=true, yminorticksvisible=true, xtickalign=1, ytickalign=1,
xminortickalign=1, yminortickalign=1, xticksize=14, xtickwidth=2.5, yticksize=14,
ytickwidth=2.5, xminorticksize=7, xminortickwidth=2.5, yminorticksize=7, yminortickwidth=2.5,
xticklabelsize=35, yticklabelsize=35, xticksmirrored=true, yticksmirrored=true),
Legend=(patchsize=(70, 10), framevisible=false, patchlabelgap=20, rowgap=10))
set_theme!(basic_theme)Compare against polytropes
Below we see how the profile of the star compares to different polytropes. We make use of the facility tools to obtain DataFrame objects out of the hdf5 output. In particular, get_profile_names_from_hdf5 will provide the names of all profiles contained within the hdf5 file, while get_profile_dataframe_from_hdf5 is used to obtain one DataFrame corresponding to one stellar profile. The animation is constructed using the Observable type that makie provides. Note that the zero points of the polytropes are arbitrary.
profile_names = StellarModels.get_profile_names_from_hdf5("profiles.hdf5")
f = Figure();
ax = Axis(f[1, 1]; xlabel=L"\log_{10}(\rho/\mathrm{[g\;cm^{-3}]})", ylabel=L"\log_{10}(P/\mathrm{[dyn]})")
pname = Observable(profile_names[1])
profile = @lift(StellarModels.get_profile_dataframe_from_hdf5("profiles.hdf5", $pname))
#To see why this is done this way, see https://docs.makie.org/stable/explanations/nodes/index.html#problems_with_synchronous_updates
#the main issue is that remeshing changes the size of the arrays
log10_ρ = @lift($profile[!, "log10_ρ"])
log10_P = Observable(rand(length(log10_ρ.val)))
profile_line = lines!(ax, log10_ρ, log10_P; label="real profile")
xvals = LinRange(-13, 4, 100)
lines!(ax, xvals, (1 + 1 / 1) .* xvals .+ 20; label="n=1")
lines!(ax, xvals, (1 + 1 / (1.5)) .* xvals .+ 15; label="n=1.5")
lines!(ax, xvals, (1 + 1 / 3) .* xvals .+ 15; label="n=3")
axislegend(ax; position=:rb)
model_number_str = @lift("model number=$(parse(Int,$pname))")
profile_text = text!(ax, -10, 20; text=model_number_str)
record(f, "rho_P_evolution.gif", profile_names[1:end]; framerate=2) do profile_name
profile = StellarModels.get_profile_dataframe_from_hdf5("profiles.hdf5", profile_name)
log10_P.val = profile[!, "log10_P"]
pname[] = profile_name
end"rho_P_evolution.gif"
Check nuclear burning
We see that the structure evolves towards an n=3 polytrope. Deviations near the core are due to the non-homogeneous composition as hydrogen is burnt. We can similarly visualize how the hydrogen mass fraction changes in the simulation. In here, only one frame shows the hydrogen that was burnt. To better visualize that you can adjust profile_interval in the IO options (and probably adjust the framerate).
profile_names = StellarModels.get_profile_names_from_hdf5("profiles.hdf5")
f = Figure();
ax = Axis(f[1, 1]; xlabel=L"\mathrm{Mass}\;[M_\odot]", ylabel=L"Abundance")
pname = Observable(profile_names[1])
profile = @lift(StellarModels.get_profile_dataframe_from_hdf5("profiles.hdf5", $pname))
mass = @lift($profile[!, "mass"])
X = Observable(rand(length(mass.val)))
Y = Observable(rand(length(mass.val)))
model_number_str = @lift("model number=$(parse(Int,$pname))")
profile_line = lines!(ax, mass, X; label="X")
profile_line = lines!(ax, mass, Y; label="Y")
profile_text = text!(ax, 0.7, 0.0; text=model_number_str)
axislegend(ax; position=:rb)
record(f, "X_evolution.gif", profile_names[1:end]; framerate=2) do profile_name
profile = StellarModels.get_profile_dataframe_from_hdf5("profiles.hdf5", profile_name)
X.val = profile[!, "X"]
Y.val = profile[!, "Y"]
pname[] = profile_name
end"X_evolution.gif"
Plot a funny HR diagram
Finally, we can also access the history data of the simulation. We use this to plot a simple HR diagram. As our microphysics are very simplistic, and the initial condition is not very physical, this looks a bit funny!
f = Figure();
ax = Axis(f[1, 1]; xlabel=L"\log_{10}(T_\mathrm{eff}/[K])", ylabel=L"\log_{10}(L/L_\odot)", xreversed=true)
history = StellarModels.get_history_dataframe_from_hdf5("history.hdf5")
lines!(ax, log10.(history[!, "T_surf"]), log10.(history[!, "L_surf"]))
f
Perform some cleanup
Internally we want to prevent storing any of the hdf5 files into our git repos, so I remove them. You can also take advantage of julia as a scripting language to post-process your simulation output in a similar way.
rm("history.hdf5")
rm("profiles.hdf5")
rm("example_options.toml")This page was generated using Literate.jl.