Quick Start Tutorial
(Each numbered section can be run independently)
Total runtime: ~4 min
Time to read: ~30 min
Note: Units are cgs unless otherwise specified. It is not necessary to use Astropy units.
Reminder: Relaxation codes are inherently finicky! Have patience and try ramping parameters in different orders or taking smaller steps in parameter space if the solver gets stuck.
Be aware that extremely high and low gravity planets push the physical and numerical limits of this model, so caution is suggested with such solutions.
[28]:
from wind_ae.wrapper.relax_wrapper import wind_simulation as wind_sim
from wind_ae.wrapper.wrapper_utils.plots import energy_plot, six_panel_plot, quick_plot, spectrum_plot
from wind_ae.wrapper.wrapper_utils import constants as const
from wind_ae.wrapper.wrapper_utils.system import system
from wind_ae.wrapper.wrapper_utils.spectrum import spectrum
import wind_ae.McAstro.atoms.atomic_species as McAtom
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
%config InlineBackend.figure_format='retina'
1. Ramping to a new planet
Load Initial Guess
Load initial guess from which the relaxation code will ramp to your desired parameters. calc_postfacto takes the time to compute all of the post-facto variables. Post facto calculations are not necessary to run the ramping.
Note: sim.load_planet('saves/filename') will load the guess parameters into the input files in Inputs/ and may interrupt ongoing runs of the relaxation code. To do computations with or plot a solution without ovewritting the Inputs/ files, use sim.load_uservars('saves/filename'). This is useful when checking the outputs of currently running calculations (e.g., sim.load_uservars('wind_ae/saves/windsoln.csv')) or to perform analysis while Wind-AE is solving for another
planet (sim.load_uservars('saves/[path-to-a-prior-solution]')).
[29]:
sim = wind_sim()
sim.load_planet('../saves/HD209458b/HD209_13.6-2000eV_H-He.csv')
sim.system.print_system()
Atmosphere Compositionrom 4 to 2. Remaking C code...
Species: HI, HeI
Mass frac: 8.00e-01, 2.00e-01
Loaded Planet:
System parameters (cgs) # Normalized units
Mp: 1.330000e+30 g # 0.70 MJ
Rp: 1.000000e+10 cm # 1.40 RJ
Mstar: 1.988416e+33 g # 1.00 Msun
semimajor: 7.480000e+11 cm # 0.05 au
Ftot: 1.095713e+03 erg/cm^2/s # 2.35e+02 FuvEarth
Lstar: 6.807420e+33 erg/s # 1.78e+00 Lsun
To run relaxation code (NOT necessary for ramping)
To run relaxation code for whatever solution is loaded into the wind_ae/inputs/ files, run sim.run_wind() which returns 0 if successfully converges on solution.
Any time a simulation is run through the python wrapper, the solution generated at each step is saved in wind_ae/saves/windsoln.csv and the solution is copied into wind_ae/inputs/guess.inp to be the guess for for the next iteration.
Ramping to “new planet”
New planets can be created by supplying a planetary mass (\(M_\mathrm{p}\)), planetary radius (\(R_\mathrm{p}\)), stellar mass (\(M_\star\)), semimajor-axis (\(a\)), stellar bolometric luminosity (\(L_\star\)) and the total integrated ionizing flux (\(F_\mathrm{tot}\)) to the system object.
To run relaxation code for whatever solution is loaded into the inputs/ files, run sim.run_wind() which returns 0 if successful.
A ramping output of 5 also means that it ramped successfully, but that one or more of the polishing steps could not be completed. Effects on mass loss rates and wind structure are minimal, but confirm that the result is acceptable by plotting the energy_plot() and six_panel_plot() or quick_plot().
Helpful Hint: Take advantage of the McAstro cgs constants object (const). These do not contain units and no functions accept astropy units.
[30]:
example = system(0.8*const.Mjupiter, #planet mass
1.5*const.Rjupiter, #planet optical transit radius
1.988416e+33, #stellar mass
748000000000.0, #semimajor axis
1095.712977126, #stellar flux at planet's semimajor axis
6.80742e+33, #stellar bolometric luminosity -
# (used to physically approximate lower boundary conditions,
# e.g., temp and density at 1 microbar)
name='example planet')
example.print_system(norm='Jupiter')
example planet:
System parameters (cgs) # Normalized units
Mp: 1.518819e+30 g # 0.80 MJ
Rp: 1.072380e+10 cm # 1.50 RJ
Mstar: 1.988416e+33 g # 1.00 Msun
semimajor: 7.480000e+11 cm # 0.05 au
Ftot: 1.095713e+03 erg/cm^2/s # 2.35e+02 FuvEarth
Lstar: 6.807420e+33 erg/s # 1.78e+00 Lsun
[31]:
%%time
sim.ramp_to(system=example,intermediate_converge_bcs=True,integrate_out=True,final_polish=True,make_plot=False)
print("\n\n Mass loss rate using sim.ramp_to(intermediate_converge_bcs=True):",
f" {sim.windsoln.Mdot:.4e} g/s\n")
Ramping Ftot from 1.096e+03 to 1.096e+03.
Ftot already done.
Ramping Mp from 1.330e+30 to 1.519e+30 g AND Rp from 1.000e+10 to 1.072e+10 cm.
Final: Mp:1.518819e+30 & Rp:1.072380e+10, M_delta:0.07652
Polishing up boundary conditions...
Successfully ramped base boundary conditions.
Ramping Mstar from 1.988e+33 to 1.988e+33 AND semimajor from 7.480e+11 to 7.480e+11 AND Lstar from 6.807e+33 to 6.807e+33.
Mstar already done.
semimajor already done.
Lstar already done.
Polishing up boundary conditions...
Successfully ramped base boundary conditions.
Successfully converged Rmax to 6.103840 (Ncol also converged).
Mass loss rate using sim.ramp_to(intermediate_converge_bcs=True): 1.9050e+10 g/s
CPU times: user 8.03 s, sys: 2.95 s, total: 11 s
Wall time: 37.1 s
To save time at a small accuracy cost, set intermediate_converge_bcs=False, though if simulation gets stuck, polishing boundary conditions between each ramping function by setting intermediate_converge_bcs=True may be advisable.
Since many of the use cases of this model require multiple iterations with changing just a few parameters (e.g., grids, evolutions, etc.) it is possible to slightly speed up computations by ramp only the relevant variables using:
ramp_Ftot(goal_Ftot,converge_bcs=)- to ramp just stellar fluxramp_grav(system=example,converge_bcs=)- to ramp just planet mass and/or radiusramp_star(system=example,converge_bcs=)- to ramp semimajor axis, Mstar, or Lstar
Speed comparisons* follow:
\(*\) Need to reload HD209 to compare apples to apples for speed. This is because when sim runs, it loads most recent solution to be the next guess. If we ramped without reloading the HD209 solution as a guess, we would be ramping from the recently solved-for example planet with 0.8 :math:`M_J`, 1.5 :math:`R_J` to the same planet with 0.8 :math:`M_J`, 1.5 :math:`R_J`.
[32]:
#Just turning off intermediate polishing
sim.load_planet('../saves/HD209458b/HD209_13.6-2000eV_H-He.csv',
calc_postfacto=False,print_atmo=False)
#expediting because we don't need to do the post facto calculations/
%time sim.ramp_to(system=example,intermediate_converge_bcs=False,integrate_out=True,final_polish=True)
print("\n\nMass loss rate using sim.ramp_to(intermediate_converge_bcs=False):",
f" {sim.windsoln.Mdot:.4e} g/s\n")
Ramping Ftot from 1.096e+03 to 1.096e+03.
Ftot already done.
Ramping Mp from 1.330e+30 to 1.519e+30 g AND Rp from 1.000e+10 to 1.072e+10 cm.
Final: Mp:1.518819e+30 & Rp:1.072380e+10, M_delta:0.07652
Ramping Mstar from 1.988e+33 to 1.988e+33 AND semimajor from 7.480e+11 to 7.480e+11 AND Lstar from 6.807e+33 to 6.807e+33.
Mstar already done.
semimajor already done.
Lstar already done.
Polishing up boundary conditions...
Successfully ramped base boundary conditions.
Successfully converged Rmax to 6.127091 (Ncol also converged).
CPU times: user 4.54 s, sys: 1.52 s, total: 6.05 s
Wall time: 15.4 s
Mass loss rate using sim.ramp_to(intermediate_converge_bcs=False): 1.9680e+10 g/s
[33]:
#Turning off all polishing (Fastest)
sim.load_planet('../saves/HD209458b/HD209_13.6-2000eV_H-He.csv',
calc_postfacto=False,print_atmo=False)
#expediting because we don't need to do the post facto calculations/
%time sim.ramp_to(system=example,intermediate_converge_bcs=False,integrate_out=False,final_polish=False)
print("\n\nMass loss rate using sim.ramp_to(intermediate_converge_bcs=False):",
f" {sim.windsoln.Mdot:.4e} g/s\n")
Ramping Ftot from 1.096e+03 to 1.096e+03.
Ftot already done.
Ramping Mp from 1.330e+30 to 1.519e+30 g AND Rp from 1.000e+10 to 1.072e+10 cm.
Final: Mp:1.518819e+30 & Rp:1.072380e+10, M_delta:0.07652
Ramping Mstar from 1.988e+33 to 1.988e+33 AND semimajor from 7.480e+11 to 7.480e+11 AND Lstar from 6.807e+33 to 6.807e+33.
Mstar already done.
semimajor already done.
Lstar already done.
CPU times: user 25.6 ms, sys: 167 ms, total: 192 ms
Wall time: 4.17 s
Mass loss rate using sim.ramp_to(intermediate_converge_bcs=False): 1.9323e+10 g/s
Since we are just ramping \(M_P\) and \(R_P\), these values are handled by sim.ramp_grav() so we can actually forgo the sim.ramp_to(system) entirely. We can always polish the boundary conditions to self-convergence later by running sim.polish_bcs()
[34]:
sim.load_planet('../saves/HD209458b/HD209_13.6-2000eV_H-He.csv',print_atmo=False)
#expediting because we don't need to do the post facto calculations
%time sim.ramp_grav(system=example,converge_bcs=False,make_plot=False)
print("\n\nMass loss rate using sim.ramp_grav(converge=False):",
f" {sim.windsoln.Mdot:.4e} g/s\n")
Ramping Mp from 1.330e+30 to 1.519e+30 g AND Rp from 1.000e+10 to 1.072e+10 cm.
Final: Mp:1.518819e+30 & Rp:1.072380e+10, M_delta:0.07652
CPU times: user 33.2 ms, sys: 171 ms, total: 204 ms
Wall time: 4.35 s
Mass loss rate using sim.ramp_grav(converge=False): 1.9323e+10 g/s
A note about Wind-AE’s ramping scheme
In the case of this example, Wind-AE likely could have ramped in just one step, since the goal was so close in parameter space. Future versions of Wind-AE may include a smarter ramping algorithm, but, for the time being, the ramping algorithm is designed to take adaptive step sizes that allow to more robustly navigate the large jumps in parameter space of which the relaxation method is infamously intolerant.
To run manually outside of the ramping scheme, load your guess
sim = wind_sim()
sim.load_planet('path/to/planet.csv')
Then edit the parameters directly in the wind_ae/inputs/ directory (PAY ATTENTION TO UNITS) or write the input values using
sim.inputs.write_bcs()
sim.inputs.write_planet()
sim.inputs.write_flags()
sim.inputs.write_physics_params()
...
etc.
Finally, run using the command line (cd wind_ae, then ./bin/relax_ae) or in Python via
sim.run_wind()
Solution outputs (windsoln object)
To interact with wind solution object, use sim.windsoln.
Dependent variables returned by the relaxation code:
'r','v','T','rho','Ys_[species]','Ncol_[species]',E, etc.Dependent variables calculated post facto by running
sim.windsoln.add_user_vars()or loading solution from save file assim.load_uservars('saves/filename'):'mu','heat_ion'(etc.), . See below.Parameters of the solution:
sim.windsoln.Rp,sim.windsoln.E_wl,sim.windsoln.Ys_rmin,sim.windsoln.species_list, etc.Tuples of the solution parameters:
sim.windsoln.planet_tuple,sim.windsoln.bcs_tuple, etc.
[35]:
sim.windsoln.add_user_vars()
sim.windsoln.soln.columns
[35]:
Index(['r', 'rho', 'v', 'T', 'Ys_HI', 'Ys_HeI', 'Ncol_HI', 'Ncol_HeI', 'q',
'z', 'n_e', 'n_H', 'n_HI', 'n_HII', 'n_He', 'n_HeI', 'n_HeII', 'n_tot',
'mu', 'P', 'ram', 'cs', 'Mach', 'Hsc', 'heat_ion', 'heat_no2nd',
'ion_rate_HI', 'ion_rate_HeI', 'euv_ion_rate_HI', 'euv_ion_rate_HeI',
'heating_rate_HI', 'heating_rate_HeI', 'DlnP', 'mfp_hb', 'mfp_Co',
'mfp_mx', 'Kn_hb_HI', 'Kn_Co_HI', 'Kn_mx_HI', 'recomb_HI', 'advec_HI',
'recomb_HeI', 'advec_HeI', 'cool_lyman', 'cool_PdV', 'cool_cond',
'cool_rec', 'boloheat', 'bolocool', 'cool_grav', 'e_therm',
'heat_advect', 'cum_heat', 'sp_kinetic', 'sp_enthalpy', 'sp_grav',
'bern_rogue', 'bern_isen', 'bern', 'ad_prof', 'static_prof',
'static_rho', 'static_P', 'static_T', 'del_ad', 'Gamma_2', 'Gamma_3'],
dtype='object')
Get Mass Loss Rate
To account for geometric and circulation effects, the substellar \(\dot{M}\) is multiplied by a factor of 0.33 in the post-facto calculations in the python wrapper.
Given the uncertainties inherent in 1D modeling of mass loss, \(\dot{M}\) differences of a few to a few tens of percent difference are negligible.
[36]:
Mdot = sim.windsoln.Mdot
print(f"Mass loss rate of {Mdot:.3e} g/s")
Mass loss rate of 1.932e+10 g/s
Plots
[37]:
energy_plot(sim.windsoln)
[38]:
#Compare to an existing solution
#load JUST the post-facto calculations from an existing solution (does not load solution to be run or ramped)
sim_HD209 = wind_sim()
sim_HD209.load_uservars('../saves/HD209458b/HD209_13.6-2000eV_H-He.csv')
#perform the post-facto calculations for the planet currently in the windsoln object
# sim.windsoln.add_user_vars()
#alternatives:
#sim.load_uservars('saves/windsoln.csv')
#sim.load_planet('saves/windsoln.csv',expedite=False)
#To plot comparisons
ax = six_panel_plot(sim_HD209.windsoln,
label='HD 209548b',c='k')
R = sim.windsoln.Rp/const.Rjupiter
M = sim.windsoln.Mp/const.Mjupiter
six_panel_plot(sim.windsoln,
first_plotted=False,ax=ax,
label=f'Example Planet: {R:.2f}R$_J$, {M:.2f}M$_J$',
c='c')
plt.show()
***** HD 209548b: Mdot = 2.01e+10 g/s *****
***** Example Planet: 1.50R$_J$, 0.80M$_J$: Mdot = 1.93e+10 g/s *****
Save Planet
To polish self consistently polish the boundary conditions, run sim.polish_bcs(converge_Rmax=) or set sim.save_planet(polish=True).
[39]:
sim.polish_bcs()
Polishing up boundary conditions...
Successfully ramped base boundary conditions.
Successfully converged Rmax to 6.127091 (Ncol also converged).
[39]:
0
[40]:
sim.save_planet('../saves/Tutorial_examples/example_13.6-2000eV_solar_H-He.csv',
overwrite=True)
Saving ../saves/Tutorial_examples/example_13.6-2000eV_solar_H-He.csv
2. Ramping Stellar Spectrum
Changing spectrum range
Fnorm flux in units of ergs/s/cm2 at the planet in the range defined by norm_spec_range. Set the full range of your spectrum in goal_spec_range. If Fnorm=0 (default), it will normalize to the current flux value in the norm_spec_range in your spectrum. If goal_spec_range=[] (default), it will maintain the same range bounds.
[41]:
sim = wind_sim()
sim.load_planet('../saves/Tutorial_examples/example_13.6-2000eV_solar_H-He.csv',
print_atmo=False)
current_wl_range = sim.windsoln.spec_resolved
eV_range = (const.hc/current_wl_range) / const.eV
print(f"Current spectrum range is {current_wl_range}cm ([{eV_range[1]:.1f},{eV_range[0]:.1f}]eV)")
Current spectrum range is [6.00e-08 9.12e-06]cm ([13.6,2066.4]eV)
If we want to maintain or change the current flux over the specified wavelength (or energy) range, we can do so.
Note that in the case below, we are keeping the flux over the EUV range (13.6-100 eV) constant by increasing Ftot as the spectrum range increases from 13.6-2000eV to 10.3-2000 eV.
[42]:
sim.ramp_spectrum(Fnorm=0,norm_spec_range=[13.6,100],
goal_spec_range=[10.3,2000],units='eV')
Spectrum will be normalized such that sum(Flux[12.40, 91.16]nm) = 979 ergs/s/cm2.
Goal: [ 0.61992127 120.37306285] nm
Ramped spectrum wavelength range, now normalizing spectrum.
Ramping Ftot from 1.096e+03 to 1.458e+03.
Final: Ftot:1.458103e+03, delta:0.001327
[42]:
0
[43]:
spectrum_plot(sim.windsoln,
xaxis='energy',
semimajor_au=sim.windsoln.semimajor/const.au,
highlight_euv=True)
#----- Alternative method --------
# sim.load_planet('../wind_ae/saves/windsoln.csv')
# #re-loading to update spectrum range for plotting purposes
# sim.spectrum.plot(xaxis='energy',
# semimajor_au=sim.windsoln.semimajor/const.au,
# highlight_euv=True)
plt.show()
[44]:
sim.save_planet('../saves/Tutorial_examples/example_10.3-2000eV_solar_H-He.csv',
overwrite=True)
Saving ../saves/Tutorial_examples/example_10.3-2000eV_solar_H-He.csv
Ramping to custom spectrum
Wind-AE includes an automatic smoothing and binning algorithm for spectra that ensures that the flux at the ionization energies for all of the species present in the simulation is accurate.
To ramp to a custom stellar spectrum first format your XUV spectrum using format_user_spectrum(). Once formatted, it will be stored in McAstro/stars/spectrum/additional_spectra/ and will be accessed automatically by the code.
The MUSCLES Survey is a good source for X-ray and UV spectra.
[45]:
## Flux, NOT flux density
# sim.format_user_spectrum(wl = ,
# flux_1au = ,
# spectrum_name='example_star',
# comment='Spectrum for some star',
# overwrite=True)
Let’s ramp to a different spectrum. We’ll also increase the flux just to show we can.
Because we did not set make_plot=False, the intermediate ramping step solutions will be plotted. This can be turned on for any ramping function via plot=True.
[46]:
sim = wind_sim()
sim.load_planet('../saves/Tutorial_examples/example_13.6-2000eV_solar_H-He.csv',print_atmo=False)
sim.ramp_to_user_spectrum(spectrum_filename='hd189733',
species_list = ['HI','HeI'],
updated_F=2000,
goal_spec_range=[13.6,2000])
*** REMINDER: Goal spectrum (orange) oversmoothed? Lower the hires_savgol_window value. ***
...Ramping Ncol Iter 1: Current average diff 1.09e+01
Success ramping to hd189733 spectrum shape!
Now ramping flux and spectral range.
Goal: [ 0.61992127 91.16489319] nm
Ramped spectrum wavelength range, now normalizing spectrum.
Ramping Ftot from 1.096e+03 to 2.018e+03.
Trying: Ftot:2.018338e+03, delta:0.05005
/Users/m/Research/relaxed-wind_xray-and-multispecies/wind_ae/wrapper/wrapper_utils/plots.py:118: UserWarning: Attempt to set non-positive ylim on a log-scaled axis will be ignored.
ax[i][j].set_ylim([min(ylims[0], 0.9*ylo/norm[i][j]),
Final: Ftot:2.018338e+03, delta:0.05005
Success! Ramped to user-input stellar spectrum hd189733
[46]:
0
3. Metals
Adding Metals
Wind-AE’s speed \(\propto n_{species}\) so it is recommended to add metals last. If you are attempting to ramp far in parameter space, it may also be advisable to ramp with only H/He and add metals as a last step.
Future versions will include the ability to include multiple ionization states of the same element, but for now, select one ionization state per species.
This computation takes, on average, 5 minutes to add C I, N I, O I, Ne I, Mg II, Ca I.
If the ramper is having trouble ramping in metals, try raising the total column density (sim.inputs.write_bcs(*sim.windsoln.bcs_tuple[0:5],sim.windsoln.Ncol_sp*5,sim.windsoln.bcs_tuple[-1])+sim.run_wind()) or raising the stellar flux (sim.ramp_var('Ftot',goal_F,converge=False)).
[47]:
sim = wind_sim()
sim.load_planet('../saves/Tutorial_examples/example_13.6-2000eV_solar_H-He.csv')
#printing the MASS fractions for Z x solar metallicity.
sim.metallicity(['HI','HeI','CI','NI','OI','NeI','MgII','CaI','SiII'],Z=1)
Atmosphere Composition
Species: HI, HeI
Mass frac: 8.00e-01, 2.00e-01
[47]:
array([7.87038072e-01, 1.99987365e-01, 2.29481799e-03, 7.90297090e-04,
6.69189203e-03, 1.74756285e-03, 6.47899249e-04, 6.35653783e-05,
7.38528061e-04])
[48]:
%%time
#default, Z=1, but can ramp in with higher Z or using custom MASS fractions
sim.add_metals(['CI','OI'],Z=1)
Successfully ramped mass fracs. Integrating out & converging Ncol_sp... g mass fractions.
Successfully converged Rmax to 6.167191. It is recommended to run converge_Ncol_sp(integrate_out=True) as Ncol will have updated.
CPU times: user 6.08 s, sys: 1.15 s, total: 7.23 s
Wall time: 26.5 s
[48]:
0
[49]:
sim.save_planet('../saves/Tutorial_examples/example_13.6-2000eV_solar_H-He-C-O.csv',
polish=True, #can skip to save time
overwrite=True)
Polishing up boundary conditions...
Successfully ramped base boundary conditions.
Successfully converged Rmax to 6.167191 (Ncol also converged).
Saving ../saves/Tutorial_examples/example_13.6-2000eV_solar_H-He-C-O.csv
[50]:
quick_plot(sim.windsoln)
*****254.32Mearth_16.81Rearth_1095F_scaled-solar Mdot = 1.82e+10 g/s ******
/Users/m/Research/relaxed-wind_xray-and-multispecies/wind_ae/wrapper/wrapper_utils/plots.py:263: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
minrad = float(radius[soln.soln['T'] == np.min(soln.soln['T'])])
[50]:
array([[<Axes: ylabel='Density (g/cm$^3$)'>,
<Axes: ylabel='Temperature (1000 K)'>],
[<Axes: xlabel='Radius (R$_p$)', ylabel='Velocity (10 km/s)'>,
<Axes: xlabel='Radius (R$_p$)', ylabel='Ionization Fraction'>]],
dtype=object)
Ramping Metallicity
Can ramp in units of \(Z_{\odot}\) or to a custom mass fraction array.
[51]:
sim = wind_sim()
sim.load_planet('../saves/Tutorial_examples/example_13.6-2000eV_solar_H-He-C-O.csv')
Atmosphere Composition
Species: HI, HeI, CI, OI
Mass frac: 7.91e-01, 2.00e-01, 2.29e-03, 6.69e-03
[52]:
sim.ramp_metallicity(Z=10) #can turn off polishing to save time
Starting metallicity: 1 xSolar
Success! Ramped to goal metallicity, Z = 10 x Solar
Integrating out and attempting to converge Ncol_sp...
Successfully converged Rmax to 6.078954 (Ncol also converged).
[52]:
0
[53]:
sim.save_planet('../saves/Tutorial_examples/example_13.6-2000eV_solar_H-He-C-O_10Z.csv',
polish=False,
overwrite=True)
Saving ../saves/Tutorial_examples/example_13.6-2000eV_solar_H-He-C-O_10Z.csv
4. Saving in an easy output format
Be nice to your collaborators and don’t send them the complicated Wind-AE solution files. Instead, send a csv with just the columns you care about.
[54]:
sim.easy_output_file(outputs=['v','T','rho','heat_ion'],
output_file='../saves/simplified_outputs/example.dat',
comments='Tutorial example planet w/ 10xZ_solar H,He,C,O. Scaled solar spectrum.',
overwrite=True)