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
Load the snapshot on a coarse box and take a projection plot
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()

Find the density peak to serve as a center for the zoomed box
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]]
Load the zoomed snap, and ensure that the bulk of the SFR is still encapsulated
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()

SAME STUFF BUT DIFFERENT SNAPSHOT; compressed for readability
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 [ ]: