yt
3.0¶
# Run this cell to load the slideshow css.
%run talktools
def fizzbuzz(number):
if number % 3 == 0 and number % 5 == 0:
return 'FizzBuzz'
elif number % 3 == 0:
return 'Fizz'
elif number % 5 == 0:
return 'Buzz'
else:
return number
fb = []
for number in range(1, 100):
fb.append(fizzbuzz(number))
print fb
[1, 2, 'Fizz', 4, 'Buzz', 'Fizz', 7, 8, 'Fizz', 'Buzz', 11, 'Fizz', 13, 14, 'FizzBuzz', 16, 17, 'Fizz', 19, 'Buzz', 'Fizz', 22, 23, 'Fizz', 'Buzz', 26, 'Fizz', 28, 29, 'FizzBuzz', 31, 32, 'Fizz', 34, 'Buzz', 'Fizz', 37, 38, 'Fizz', 'Buzz', 41, 'Fizz', 43, 44, 'FizzBuzz', 46, 47, 'Fizz', 49, 'Buzz', 'Fizz', 52, 53, 'Fizz', 'Buzz', 56, 'Fizz', 58, 59, 'FizzBuzz', 61, 62, 'Fizz', 64, 'Buzz', 'Fizz', 67, 68, 'Fizz', 'Buzz', 71, 'Fizz', 73, 74, 'FizzBuzz', 76, 77, 'Fizz', 79, 'Buzz', 'Fizz', 82, 83, 'Fizz', 'Buzz', 86, 'Fizz', 88, 89, 'FizzBuzz', 91, 92, 'Fizz', 94, 'Buzz', 'Fizz', 97, 98, 'Fizz']
yt
?¶yt
is a python library for analysis and visualization environment of volumetric data
website('yt-project.org', height=600)
yt
relies on the incredibly rich scientific python ecosystem¶"If I have seen further it is by standing on the shoulders of giants."
Patch AMR | Oct AMR | SPH and N-body | Hexahedral mesh | Misc |
---|---|---|---|---|
enzo | Ramses | Gadget | MOAB | Unigrid ndarray |
FLASH | ART | PKDGRAV | - | Patch AMR |
Boxlib | ARTIO | Gasoline | - | hexahedral mesh |
Chombo | - | - | - | Particle data |
Athena | - | - | - | FITS files |
Piernik | - | - | - |
The molecular outflow in the Sculptor galaxy, imaged in CO 1→0 emission by ALMA. Visualization by Erik Rosolowsky.
A 48003 simulation of reionization with enzo. Courtesy Skillman & Norman.
A giant star undergoing a close passage with a supermassive black hole. Macleod et al. (2013)
HTML('<iframe width="800" height="640" frameborder="0" allowFullScreen webkitallowfullscreen mozallowfullscreen src="https://sketchfab.com/models/486421c20a63434d9267020c11a4bc31/embed"></iframe><p style="font-size: 13px; font-weight: normal; margin: 5px; color: #4A4A4A;"><a href="https://sketchfab.com/models/486421c20a63434d9267020c11a4bc31" style="font-weight: bold; color: #1CAAD9;">Stellar Winds isodensity contour, rho=1e-24 g/cc</a>by <a href="https://sketchfab.com/jnaiman" style="font-weight: bold; color: #1CAAD9;">jnaiman</a>on <a href="https://sketchfab.com" style="font-weight: bold; color: #1CAAD9;">Sketchfab</a></p>')
website('http://yt-project.org/gallery.html')
yt
scripts look a little strange"¶yt
3.0 includes a rethought API and user interface. If you're used to yt
2.x, things might look a bit different.
Some things have changed (for the better) so bear with me. I'll try to point out differences for any yt users in the room.
yt
script¶# load yt, as simple as:
import yt
# yt.load() automagically figures out the output type
ds = yt.load('HiresGalaxy/DD0100/DD0100')
# Visualize a slice through the midplane of the galaxy
slc = yt.SlicePlot(ds, 'z', 'density', width=20*yt.units.kiloparsec)
# Customize the plot a bit
slc.set_figure_size(4)
import yt
ds = yt.load('big_cosmo/DD0020/DD0020')
print "Domain width is %s" % ds.domain_width.in_cgs()
print "Current redshift is %s" % ds.current_redshift
slc = yt.SlicePlot(ds, 0, 'density')
slc.set_figure_size(5)
Domain width is [ 2.41197986e+25 2.41197986e+25 2.41197986e+25] cm Current redshift is 5.36022168546
# A data object that selects all of the data in the simulation, regardless of AMR level.
ad = ds.all_data()
# All of the gas density values, equivalent to ad['density']
dens = ad[('gas', 'density')]
print dens
print len(dens)
print type(dens)
[ 3.47961414e-28 2.57449887e-28 1.85669820e-28 ..., 8.01621880e-28 7.76357441e-28 5.59644410e-28] g/cm**3 161766900 <class 'yt.units.yt_array.YTArray'>
%matplotlib inline
from matplotlib import pyplot as plt
# Each value corresponds to a single cell - this dataset includes a large amount of data on 4 AMR levels
grid_level = ad['grid_level'].v
plt.hist(grid_level, bins=[-0.5, 0.5, 1.5, 2.5, 3.5, 4.5], log=True, facecolor='green')
plt.xlabel('AMR Level')
plt.ylabel('Number of cells')
<matplotlib.text.Text at 0x2b79fca21290>
max_val, max_loc = ds.find_max('matter_density')
print max_val
print max_loc
9.67700879564e-24 g/cm**3 [ 0.4833374 0.50494385 0.48797607] code_length
from yt.units import Mpc
sphere = ds.sphere(center=max_loc, radius=1*Mpc)
prj = yt.ProjectionPlot(ds, 0, ('gas', 'density'), center=max_loc, width=2.2*Mpc, data_source=sphere)
prj.set_zlim(('gas', 'density'), 1e-5, 1e-2)
prj.set_figure_size(5)
# Note: using a "deposit" field along with a weight field
# ('deposit', 'all_cic') is the cloud-in-cell deposition of all of the particles in the simulatin
prj = yt.ProjectionPlot(ds, 0, ('deposit', 'all_cic'), center=max_loc, width=2.2*Mpc, data_source=sphere, weight_field=('deposit', 'all_cic'))
prj.set_figure_size(5)
/pfs/goldbaum/yt-hg/yt/data_objects/construction_data_containers.py:302: RuntimeWarning: invalid value encountered in divide np.divide(nvals, nwvals[:,None], nvals) /pfs/goldbaum/anaconda/lib/python2.7/site-packages/matplotlib/colors.py:986: RuntimeWarning: invalid value encountered in less_equal mask |= resdat <= 0
yt
recognizes two types of fields:¶"on-disk" fields
derived fields
from pprint import pprint
ds = yt.load("GalaxyClusterMerger/fiducial_1to3_b0.273d_hdf5_plt_cnt_0175")
pprint(ds.field_list)
[('flash', 'pden'), ('flash', 'gpot'), ('flash', 'dens'), ('flash', 'temp'), ('flash', 'clr1'), ('flash', 'clr2'), ('flash', 'velx'), ('flash', 'vely'), ('flash', 'velz')]
ds.derived_field_list
[('flash', 'clr1'), ('flash', 'clr2'), ('flash', 'dens'), ('flash', 'gpot'), ('flash', 'pden'), ('flash', 'temp'), ('flash', 'velx'), ('flash', 'vely'), ('flash', 'velz'), ('gas', 'angular_momentum_magnitude'), ('gas', 'angular_momentum_x'), ('gas', 'angular_momentum_y'), ('gas', 'angular_momentum_z'), ('gas', 'averaged_density'), ('gas', 'cell_mass'), ('gas', 'cutting_plane_velocity_x'), ('gas', 'cutting_plane_velocity_y'), ('gas', 'cutting_plane_velocity_z'), ('gas', 'cylindrical_radial_velocity'), ('gas', 'cylindrical_radial_velocity_absolute'), ('gas', 'cylindrical_tangential_velocity'), ('gas', 'cylindrical_tangential_velocity_absolute'), ('gas', 'density'), ('gas', 'density_gradient_magnitude'), ('gas', 'density_gradient_x'), ('gas', 'density_gradient_y'), ('gas', 'density_gradient_z'), ('gas', 'dynamical_time'), ('gas', 'entropy'), ('gas', 'gravitational_potential'), ('gas', 'kT'), ('gas', 'kinetic_energy'), ('gas', 'radial_velocity'), ('gas', 'radial_velocity_absolute'), ('gas', 'shear'), ('gas', 'specific_angular_momentum_magnitude'), ('gas', 'specific_angular_momentum_x'), ('gas', 'specific_angular_momentum_y'), ('gas', 'specific_angular_momentum_z'), ('gas', 'sz_kinetic'), ('gas', 'szy'), ('gas', 'tangential_over_velocity_magnitude'), ('gas', 'tangential_velocity'), ('gas', 'temperature'), ('gas', 'velocity_divergence'), ('gas', 'velocity_divergence_absolute'), ('gas', 'velocity_magnitude'), ('gas', 'velocity_x'), ('gas', 'velocity_y'), ('gas', 'velocity_z'), ('gas', 'vorticity_magnitude'), ('gas', 'vorticity_squared'), ('gas', 'vorticity_stretching_magnitude'), ('gas', 'vorticity_stretching_x'), ('gas', 'vorticity_stretching_y'), ('gas', 'vorticity_stretching_z'), ('gas', 'vorticity_x'), ('gas', 'vorticity_y'), ('gas', 'vorticity_z'), ('gas', 'xray_emissivity'), ('index', 'cell_volume'), ('index', 'contours'), ('index', 'cylindrical_r'), ('index', 'cylindrical_theta'), ('index', 'cylindrical_z'), ('index', 'disk_angle'), ('index', 'dx'), ('index', 'dy'), ('index', 'dz'), ('index', 'grid_indices'), ('index', 'grid_level'), ('index', 'height'), ('index', 'ones'), ('index', 'ones_over_dx'), ('index', 'radius'), ('index', 'spherical_phi'), ('index', 'spherical_r'), ('index', 'spherical_theta'), ('index', 'x'), ('index', 'y'), ('index', 'z'), ('index', 'zeros')]
fi = ds.field_info[('gas', 'cell_mass')]
print fi.get_source()
def _cell_mass(field, data): return data[ftype, "density"] * data["index", "cell_volume"]
print ds.field_info[('gas', 'specific_angular_momentum_x')].get_source()
def _specific_angular_momentum_x(field, data): xv, yv, zv = obtain_velocities(data, ftype) center = data.get_field_parameter('center') v_vec = obtain_rvec(data) v_vec = np.rollaxis(v_vec, 0, len(v_vec.shape)) v_vec = data.pf.arr(v_vec, input_units = data["index", "x"].units) rv = v_vec - center return yv * rv[...,2] - zv * rv[...,1]
def three_times_density(field, data):
return 3*data['density']
ds = yt.load('HiresGalaxy/DD0100/DD0100')
ds.add_field(("gas", "three_times_density"), units="g/cm**3", function=three_times_density)
ad = ds.all_data()
print ad['density']
print ad['three_times_density']
[ 1.28845778e-30 1.28909307e-30 1.28917292e-30 ..., 1.38437060e-25 1.17462462e-25 1.06280662e-25] g/cm**3 [ 3.86537335e-30 3.86727922e-30 3.86751875e-30 ..., 4.15311179e-25 3.52387385e-25 3.18841985e-25] g/cm**3
from yt.units import cm, meter, kilogram, gram, second, joule
print kilogram
print type(kilogram)
1.0 kg <class 'yt.units.yt_array.YTQuantity'>
print kilogram*meter**2/second**2
1.0 kg*m**2/s**2
print kilogram*meter**2/second**2 == joule
print joule
True 1.0 J
print kilogram + gram
1.001 kg
from yt.utilities.physical_constants import G, kboltz
print "Newton's constant: ", G
print "Newton's constant in MKS: ", G.in_mks(), "\n"
print "Boltzmann constant: ", kboltz
print "Boltzmann constant in MKS: ", kboltz.in_mks()
Newton's constant: 6.67384e-08 cm**3/(g*s**2) Newton's constant in MKS: 6.67384e-11 m**3/(kg*s**2) Boltzmann constant: 1.3806488e-16 erg/K Boltzmann constant in MKS: 1.3806488e-23 kg*m**2/(K*s**2)
import numpy as np
rho = 1e-23*gram/cm**3
t_ff = np.sqrt(3*np.pi/(32*G*rho))
print t_ff
print t_ff.in_units('Myr')
6.64312889844e+14 s 21.0508051894 Myr
ds = yt.load('big_cosmo/DD0020/DD0020')
prj = yt.ProjectionPlot(ds, 0, 'matter_density', width=(32, 'Mpccm/h'), weight_field='matter_density')
prj.set_figure_size(6)
prj.set_unit('matter_density', 'Msun/pc**3')
%matplotlib inline
import yt
import numpy as np
from yt.data_objects.particle_filters import add_particle_filter
from matplotlib import pyplot as plt
from pprint import pprint
ds = yt.load('HiresGalaxy/DD0100/DD0100')
pprint(sorted(ds.field_list))
[('all', 'creation_time'), ('all', 'dynamical_time'), ('all', 'metallicity_fraction'), ('all', 'particle_index'), ('all', 'particle_mass'), ('all', 'particle_position_x'), ('all', 'particle_position_y'), ('all', 'particle_position_z'), ('all', 'particle_type'), ('all', 'particle_velocity_x'), ('all', 'particle_velocity_y'), ('all', 'particle_velocity_z'), ('enzo', 'Cooling_Time'), ('enzo', 'Dark_Matter_Density'), ('enzo', 'Density'), ('enzo', 'GasEnergy'), ('enzo', 'Metal_Density'), ('enzo', 'Temperature'), ('enzo', 'TotalEnergy'), ('enzo', 'x-velocity'), ('enzo', 'y-velocity'), ('enzo', 'z-velocity'), ('io', 'creation_time'), ('io', 'dynamical_time'), ('io', 'metallicity_fraction'), ('io', 'particle_index'), ('io', 'particle_mass'), ('io', 'particle_position_x'), ('io', 'particle_position_y'), ('io', 'particle_position_z'), ('io', 'particle_type'), ('io', 'particle_velocity_x'), ('io', 'particle_velocity_y'), ('io', 'particle_velocity_z')]
def formed_star(pfilter, data):
filter = data["all", "creation_time"] > 0
return filter
add_particle_filter("formed_star", function=formed_star, filtered_type='all',
requires=["creation_time"])
# Load the data
ds = yt.load("HiresGalaxy/DD0100/DD0100")
ds.add_particle_filter('formed_star')
# Create a selector to access the data
ad = ds.all_data()
# Access the creation time and particle mass datasets.
masses = ad['formed_star', 'particle_mass'].in_units('Msun')
formation_time = ad['formed_star', 'creation_time'].in_units('yr')
# Create a plot of the star formation rate as a function of time
time_range = [0, 1e8] # years
n_bins = 1000
hist, bins = np.histogram(formation_time, bins=n_bins, range=time_range,)
time = (bins[:-1] + bins[1:])/2
inds = np.digitize(formation_time, bins=bins)
sfr = np.array([masses[inds == j].sum()/(bins[j+1]-bins[j]) for j in range(len(time))])
sfr[sfr == 0] = np.nan
plt.plot(time/1e6, sfr)
plt.xlabel('Time [Myr]')
plt.ylabel('SFR [M$_\odot$ yr$^{-1}$]')
<matplotlib.text.Text at 0x2b7b3aeafd10>
# Can create deposition fields based on these new particle types
yt.SlicePlot(ds, 2, ('deposit', 'formed_star_cic'), width=(20, 'kpc'))
from IPython.display import YouTubeVideo
YouTubeVideo('YWi4hO76Wss', width=1024, height=768)
YouTubeVideo('cyzzxkX4gjw', width=1024, height=768)
%%writefile parallel_test.py
import yt
yt.enable_parallelism()
ds = yt.load('big_cosmo/DD0020/DD0020')
ad = ds.all_data()
tm = ad.quantities.total_mass()
if yt.is_root():
print tm
Overwriting parallel_test.py
cat parallel_test.py
import yt yt.enable_parallelism() ds = yt.load('big_cosmo/DD0020/DD0020') ad = ds.all_data() tm = ad.quantities.total_mass() if yt.is_root(): print tm
%%time
!python parallel_test.py
[1.52960144299e+48 g, 0.832720588232 code_mass] CPU times: user 271 ms, sys: 66 ms, total: 337 ms Wall time: 1min 44s
%%time
!mpirun -n 4 python parallel_test.py
[1.52960144299e+48 g, 0.832720588232 code_mass] CPU times: user 118 ms, sys: 31 ms, total: 149 ms Wall time: 26.9 s
website('http://yt-project.org/docs/dev-3.0/analyzing/parallel_computation.html')
import yt
import glob
from IPython.html.widgets import interact
from IPython.display import display
from yt.units import kpc
fns = sorted(glob.glob('/trove/goldbaum/nofeedback/DD????/DD????'))
time_series = yt.DatasetSeries(fns)
ds = time_series[-1]
prj = yt.ProjectionPlot(ds, 2, 'density')
@interact(width=[1, 100])
def show_plot(width=100):
prj.set_width(width*kpc)
display(prj)
@interact(index=[0, len(fns)-1])
def show_plot(index=57):
ds = time_series[index]
prj = yt.ProjectionPlot(ds, 2, 'density', width=(20, 'kpc'))
prj.set_figure_size(5)
prj.set_zlim('density', 1e-4, 1e-1)
display(prj)