In [1]:
from yt.mods import *
import numpy as np
import pdb
import matplotlib.pyplot as plt
import yt.units as units
from yt.units import kpc
from yt.frontends.sph.data_structures import ParticleDataset
ParticleDataset.filter_bbox = True
ParticleDataset._skip_cache = True
In [2]:
fname = '/Volumes/pegasus/gadgetruns/m13_mr_Dec16_2013/snapdir_257/snapshot_257.0.hdf5'
unit_base = {'UnitLength_in_cm' : 3.08568e+21,
'UnitMass_in_g' : 1.989e+43,
'UnitVelocity_in_cm_per_s' : 100000}
bbox_lim = 1.e5
bbox = [[-bbox_lim,bbox_lim],
[-bbox_lim,bbox_lim],
[-bbox_lim,bbox_lim]]
ds1 = load(fname,unit_base=unit_base,bounding_box=bbox)
ds1.index
ad1= ds1.all_data()
In [3]:
print 'left edge: ',ds1.domain_left_edge
print 'right edge: ',ds1.domain_right_edge
print 'center: ',ds1.domain_center
print 'total SFR is [in entire volume]: ',np.sum(ad1[("PartType0","StarFormationRate")])
left edge: [-100000. -100000. -100000.] code_length right edge: [ 100000. 100000. 100000.] code_length center: [ 0. 0. 0.] code_length total SFR is [in entire volume]: 237.399819416 dimensionless
In [4]:
px = ProjectionPlot(ds1, 'y', ('deposit', 'all_density'))
px.show()
In [5]:
density = ad1[("PartType0","density")]
wdens = np.where(density == np.max(density))
coordinates = ad1[("PartType0","Coordinates")]
maxdens_coordinates = coordinates[wdens]
center = maxdens_coordinates[0]
In [6]:
bbox_lim = ds1.quan(5.e3,'code_length')
bbox = [[center[0]-bbox_lim,center[0]+bbox_lim],
[center[1]-bbox_lim,center[1]+bbox_lim],
[center[2]-bbox_lim,center[2]+bbox_lim]]
print bbox
[[24397.5429688 code_length, 34397.5429688 code_length], [29971.9140625 code_length, 39971.9140625 code_length], [23197.0175781 code_length, 33197.0175781 code_length]]
In [7]:
ds2 = load(fname,unit_base=unit_base,bounding_box=bbox)
ds2.index
ds2.periodicity=(False,False,False)
ad2= ds2.all_data()
np.sum(ad2["PartType0","StarFormationRate"])
Out[7]:
219.527018242 dimensionless
In [8]:
print ds2.domain_center
print ds2.domain_width
print 'right edge: ',ds2.domain_right_edge
print 'left edge: ',ds2.domain_left_edge
[ 29397.54296875 34971.9140625 28197.01757812] code_length [ 10000. 10000. 10000.] code_length right edge: [ 34397.54296875 39971.9140625 33197.01757812] code_length left edge: [ 24397.54296875 29971.9140625 23197.01757812] code_length
In [9]:
px = ProjectionPlot(ds2, 'x', ('deposit', 'PartType0_smoothed_density'))
px.show()
In [10]:
fname = '/Volumes/pegasus/gadgetruns/m13_mr_Dec16_2013/snapdir_239/snapshot_239.0.hdf5'
unit_base = {'UnitLength_in_cm' : 3.08568e+21,
'UnitMass_in_g' : 1.989e+43,
'UnitVelocity_in_cm_per_s' : 100000}
bbox_lim = 1.e5
bbox = [[-bbox_lim,bbox_lim],
[-bbox_lim,bbox_lim],
[-bbox_lim,bbox_lim]]
ds1 = load(fname,unit_base=unit_base,bounding_box=bbox)
ds1.index
ad1= ds1.all_data()
density = ad1[("PartType0","density")]
wdens = np.where(density == np.max(density))
coordinates = ad1[("PartType0","Coordinates")]
maxdens_coordinates = coordinates[wdens]
center = maxdens_coordinates[0]
bbox_lim = ds1.quan(5.e3,'code_length')
bbox2 = [[center[0]-bbox_lim,center[0]+bbox_lim],
[center[1]-bbox_lim,center[1]+bbox_lim],
[center[2]-bbox_lim,center[2]+bbox_lim]]
print bbox2
ds2 = load(fname,unit_base=unit_base,bounding_box=bbox2)
ds2.index
ds2.periodicity=(False,False,False)
ad2= ds2.all_data()
print "zoomed SFR: ",np.sum(ad2["PartType0","StarFormationRate"])
print ds2.domain_center
print ds2.domain_width
print 'right edge: ',ds2.domain_right_edge
print 'left edge: ',ds2.domain_left_edge
px = ProjectionPlot(ds2, 'x', ('deposit', 'PartType0_smoothed_density'))
px.show()
[[25125.9824219 code_length, 35125.9824219 code_length], [29278.6796875 code_length, 39278.6796875 code_length], [23595.5566406 code_length, 33595.5566406 code_length]] zoomed SFR:
--------------------------------------------------------------------------- YTDomainOverflow Traceback (most recent call last) <ipython-input-10-2daffd009d27> in <module>() 42 print 'left edge: ',ds2.domain_left_edge 43 ---> 44 px = ProjectionPlot(ds2, 'x', ('deposit', 'PartType0_smoothed_density')) 45 px.show() /Users/desika/yt-x86_64/src/yt-hg/yt/visualization/plot_window.pyc in __init__(self, pf, axis, fields, center, width, axes_unit, weight_field, max_level, origin, fontsize, field_parameters, data_source, proj_style, window_size, aspect) 1213 proj = pf.proj(fields, axis, weight_field=weight_field, 1214 center=center, data_source=data_source, -> 1215 field_parameters = field_parameters, style = proj_style) 1216 PWViewerMPL.__init__(self, proj, bounds, fields=fields, origin=origin, 1217 fontsize=fontsize, window_size=window_size, aspect=aspect) /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/construction_data_containers.pyc in __init__(self, field, axis, weight_field, center, pf, data_source, style, field_parameters) 227 self.data_source = data_source 228 self.weight_field = weight_field --> 229 self.get_data(field) 230 231 @property /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/construction_data_containers.pyc in get_data(self, fields) 271 with self.data_source._field_parameter_state(self.field_parameters): 272 for chunk in parallel_objects(self.data_source.chunks( --> 273 chunk_fields, "io", local_only = True)): 274 mylog.debug("Adding chunk (%s) to tree (%0.3e GB RAM)", chunk.ires.size, 275 get_memory_usage()/1024.) /Users/desika/yt-x86_64/src/yt-hg/yt/utilities/parallel_tools/parallel_analysis_interface.pyc in parallel_objects(objects, njobs, storage, barrier, dynamic) 448 # this will prevent intermediate objects from being created. 449 oiter = itertools.islice(enumerate(objects), my_new_id, None, njobs) --> 450 for obj_id, obj in oiter: 451 result_id = obj_id * njobs + my_new_id 452 if storage is not None: /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in chunks(self, fields, chunking_style, **kwargs) 527 for chunk in self.index._chunk(self, chunking_style, **kwargs): 528 with self._chunked_read(chunk): --> 529 self.get_data(fields) 530 # NOTE: we yield before releasing the context 531 yield self /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in get_data(self, fields) 631 632 fields_to_generate += gen_fluids + gen_particles --> 633 self._generate_fields(fields_to_generate) 634 for field in self.field_data.keys(): 635 if field not in ofields: /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in _generate_fields(self, fields_to_generate) 651 fi = self.pf._get_field_info(*field) 652 try: --> 653 fd = self._generate_field(field) 654 if type(fd) == np.ndarray: 655 fd = self.pf.arr(fd, fi.units) /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in _generate_field(self, field) 255 tr = self._generate_particle_field(field) 256 else: --> 257 tr = self._generate_fluid_field(field) 258 if tr is None: 259 raise YTCouldNotGenerateField(field, self.pf) /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in _generate_fluid_field(self, field) 273 finfo.check_available(gen_obj) 274 except NeedsGridType as ngt_exception: --> 275 rv = self._generate_spatial_fluid(field, ngt_exception.ghost_zones) 276 else: 277 rv = finfo(gen_obj) /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in _generate_spatial_fluid(self, field, ngz) 289 o = self._current_chunk.objs[0] 290 with o._activate_cache(): --> 291 ind += o.select(self.selector, self[field], rv, ind) 292 else: 293 chunks = self.index._chunk(self, "spatial", ngz = ngz) /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in __getitem__(self, key) 218 return self.field_data[f] 219 else: --> 220 self.get_data(f) 221 # fi.units is the unit expression string. We depend on the registry 222 # hanging off the dataset to define this unit object. /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in get_data(self, fields) 631 632 fields_to_generate += gen_fluids + gen_particles --> 633 self._generate_fields(fields_to_generate) 634 for field in self.field_data.keys(): 635 if field not in ofields: /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in _generate_fields(self, fields_to_generate) 651 fi = self.pf._get_field_info(*field) 652 try: --> 653 fd = self._generate_field(field) 654 if type(fd) == np.ndarray: 655 fd = self.pf.arr(fd, fi.units) /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in _generate_field(self, field) 255 tr = self._generate_particle_field(field) 256 else: --> 257 tr = self._generate_fluid_field(field) 258 if tr is None: 259 raise YTCouldNotGenerateField(field, self.pf) /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in _generate_fluid_field(self, field) 275 rv = self._generate_spatial_fluid(field, ngt_exception.ghost_zones) 276 else: --> 277 rv = finfo(gen_obj) 278 return rv 279 /Users/desika/yt-x86_64/src/yt-hg/yt/fields/derived_field.pyc in __call__(self, data) 176 "for %s" % (self.name,)) 177 with self.unit_registry(data): --> 178 dd = self._function(self, data) 179 for field_name in data.keys(): 180 if field_name not in original_fields: /Users/desika/yt-x86_64/src/yt-hg/yt/fields/particle_fields.pyc in _vol_weight(field, data) 519 rv = data.smooth(pos, [mass, hsml, dens, quan], 520 method="volume_weighted", --> 521 create_octree = True)[0] 522 rv[np.isnan(rv)] = 0.0 523 # Now some quick unit conversions. /Users/desika/yt-x86_64/src/yt-hg/yt/data_objects/octree_subset.pyc in smooth(self, positions, fields, index_fields, method, create_octree, nneighbors) 172 positions[:,0], positions[:,1], positions[:,2], 173 self.pf.domain_left_edge, --> 174 self.pf.domain_right_edge,False) 175 morton.sort() 176 particle_octree = ParticleOctreeContainer([1, 1, 1], /Users/desika/yt-x86_64/src/yt-hg/yt/utilities/lib/geometry_utils.so in yt.utilities.lib.geometry_utils.compute_morton (yt/utilities/lib/geometry_utils.c:5974)() YTDomainOverflow: Particle bounds (25860.53125 code_length, 29278.6796875 code_length, 23615.6601562 code_length) and (35113.1171875 code_length, 39271.1523438 code_length, 33594.4140625 code_length) exceed domain bounds [ 25125.98242188 29278.6796875 23595.55664062] code_length and [ 35125.98242188 39278.6796875 33595.55664062] code_length
219.07782428 dimensionless [ 30125.98242188 34278.6796875 28595.55664062] code_length [ 10000. 10000. 10000.] code_length right edge: [ 35125.98242188 39278.6796875 33595.55664062] code_length left edge: [ 25125.98242188 29278.6796875 23595.55664062] code_length
In [ ]: