[yt-svn] commit/yt-3.0: 17 new changesets
Bitbucket
commits-noreply at bitbucket.org
Thu Jan 10 06:56:53 PST 2013
17 new commits in yt-3.0:
https://bitbucket.org/yt_analysis/yt-3.0/commits/58aeb5270b15/
changeset: 58aeb5270b15
branch: yt-3.0
user: MatthewTurk
date: 2012-11-21 15:22:20
summary: Removing a stray save_data call that saved field lists.
affected #: 1 file
diff -r 57db505f0622e3f3de950413d8556293174aa461 -r 58aeb5270b1589dfcee71c41517b3218e0a1e52d yt/frontends/enzo/data_structures.py
--- a/yt/frontends/enzo/data_structures.py
+++ b/yt/frontends/enzo/data_structures.py
@@ -440,7 +440,6 @@
else:
field_list = None
field_list = self.comm.mpi_bcast(field_list)
- self.save_data(list(field_list),"/","DataFields",passthrough=True)
self.field_list = list(field_list)
def _generate_random_grids(self):
https://bitbucket.org/yt_analysis/yt-3.0/commits/c57573bb2025/
changeset: c57573bb2025
branch: yt-3.0
user: MatthewTurk
date: 2012-12-07 16:00:03
summary: Merge
affected #: 1 file
diff -r b739b6a59d3db78d816cc6b8769b3708aa525912 -r c57573bb2025a974faa1e6a346eae38c9135df53 yt/frontends/enzo/data_structures.py
--- a/yt/frontends/enzo/data_structures.py
+++ b/yt/frontends/enzo/data_structures.py
@@ -440,7 +440,6 @@
else:
field_list = None
field_list = self.comm.mpi_bcast(field_list)
- self.save_data(list(field_list),"/","DataFields",passthrough=True)
self.field_list = list(field_list)
def _generate_random_grids(self):
https://bitbucket.org/yt_analysis/yt-3.0/commits/a1d1f0d4be06/
changeset: a1d1f0d4be06
branch: yt-3.0
user: MatthewTurk
date: 2012-12-07 20:36:40
summary: Adding fp_utils include to oct_container.pxd. Ensuring that in the particle
octree container all octs have only a single domain.
affected #: 3 files
diff -r c57573bb2025a974faa1e6a346eae38c9135df53 -r a1d1f0d4be060b1372732c33064e0515f0623db5 yt/geometry/oct_container.pxd
--- a/yt/geometry/oct_container.pxd
+++ b/yt/geometry/oct_container.pxd
@@ -24,6 +24,7 @@
"""
cimport numpy as np
+from fp_utils cimport *
cdef struct ParticleArrays
diff -r c57573bb2025a974faa1e6a346eae38c9135df53 -r a1d1f0d4be060b1372732c33064e0515f0623db5 yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -729,6 +729,7 @@
cdef np.float64_t dds[3], cp[3], pp[3]
cdef int ind[3]
self.max_domain = max(self.max_domain, domain_id)
+ cdef int mid, mad
if self.root_mesh[0][0][0] == NULL: self.allocate_root()
for p in range(no):
level = 0
@@ -740,7 +741,7 @@
cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
if cur == NULL:
raise RuntimeError
- if cur.sd.np == 32:
+ if self._check_refine(cur, cp) == 1:
self.refine_oct(cur, cp)
while cur.sd.np < 0:
for i in range(3):
@@ -753,7 +754,7 @@
cp[i] += dds[i]/2.0
cur = cur.children[ind[0]][ind[1]][ind[2]]
level += 1
- if cur.sd.np == 32:
+ if self._check_refine(cur, cp) == 1:
self.refine_oct(cur, cp)
# Now we copy in our particle
pi = cur.sd.np
@@ -763,6 +764,16 @@
cur.sd.domain_id[pi] = domain_id
cur.sd.np += 1
+ cdef int _check_refine(self, Oct *cur, np.float64_t cp[3]):
+ cdef int mid = 16384
+ cdef int mad = -16384
+ for i in range(imax(cur.sd.np, 0)):
+ mid = imin(cur.sd.domain_id[i], mid)
+ mad = imax(cur.sd.domain_id[i], mad)
+ if cur.sd.np == 32 or mid < mad:
+ return 1
+ return 0
+
cdef void refine_oct(self, Oct *o, np.float64_t pos[3]):
cdef int i, j, k, m, ind[3]
cdef Oct *noct
diff -r c57573bb2025a974faa1e6a346eae38c9135df53 -r a1d1f0d4be060b1372732c33064e0515f0623db5 yt/geometry/setup.py
--- a/yt/geometry/setup.py
+++ b/yt/geometry/setup.py
@@ -9,6 +9,7 @@
config = Configuration('geometry',parent_package,top_path)
config.add_extension("oct_container",
["yt/geometry/oct_container.pyx"],
+ include_dirs=["yt/utilities/lib/"],
libraries=["m"],
depends=["yt/utilities/lib/fp_utils.pxd",
"yt/geometry/oct_container.pxd",
https://bitbucket.org/yt_analysis/yt-3.0/commits/3874ec822d9e/
changeset: 3874ec822d9e
branch: yt-3.0
user: MatthewTurk
date: 2012-12-07 22:26:04
summary: Wiping out most of the Gadget frontend, which leads down a rabbit hole of refactoring
affected #: 3 files
diff -r a1d1f0d4be060b1372732c33064e0515f0623db5 -r 3874ec822d9eda570abbcf277e85f2b7cbcd7175 yt/frontends/gadget/api.py
--- a/yt/frontends/gadget/api.py
+++ b/yt/frontends/gadget/api.py
@@ -28,14 +28,3 @@
"""
-from .data_structures import \
- GadgetGrid, \
- GadgetHierarchy, \
- GadgetStaticOutput
-
-from .fields import \
- GadgetFieldInfo, \
- add_gadget_field
-
-from .io import \
- IOHandlerGadget
diff -r a1d1f0d4be060b1372732c33064e0515f0623db5 -r 3874ec822d9eda570abbcf277e85f2b7cbcd7175 yt/frontends/gadget/data_structures.py
--- a/yt/frontends/gadget/data_structures.py
+++ b/yt/frontends/gadget/data_structures.py
@@ -1,13 +1,11 @@
"""
-Data structures for Gadget.
+Data structures for a generic SPH/Gadget frontend.
Author: Matthew Turk <matthewturk at gmail.com>
-Affiliation: UCSD
-Author: Chris Moody <cemoody at ucsc.edu>
-Affiliation: UCSC
+Affiliation: Columbia University
Homepage: http://yt-project.org/
License:
- Copyright (C) 2007-2011 Matthew Turk. All Rights Reserved.
+ Copyright (C) 2007-2012 Matthew Turk. All Rights Reserved.
This file is part of yt.
@@ -30,172 +28,31 @@
from itertools import izip
from yt.funcs import *
-from yt.data_objects.grid_patch import \
- AMRGridPatch
-from yt.geometry.grid_geometry_handler import \
- GridGeometryHandler
+from yt.geometry.oct_geometry_handler import \
+ OctreeGeometryHandler
+from yt.geometry.geometry_handler import \
+ GeometryHandler, YTDataChunk
from yt.data_objects.static_output import \
StaticOutput
-from yt.utilities.definitions import \
- sec_conversion
-from .fields import GadgetFieldInfo, KnownGadgetFields
from yt.data_objects.field_info_container import \
FieldInfoContainer, NullFunc
-class GadgetGrid(AMRGridPatch):
- _id_offset = 0
- def __init__(self, hierarchy, id, dimensions, start,
- level, parent_id, particle_count):
- AMRGridPatch.__init__(self, id, filename = hierarchy.filename,
- hierarchy = hierarchy)
- self.Parent = [] # Only one parent per grid
- self.Children = []
- self.Level = level
- self.ActiveDimensions = dimensions.copy()
- self.NumberOfParticles = particle_count
- self.start_index = start.copy().astype("int64")
- self.stop_index = self.start_index + dimensions.copy()
- self.id = id
- self._parent_id = parent_id
-
- try:
- padd = '/data/grid_%010i/particles' % id
- self.particle_types = self.hierarchy._handle[padd].keys()
- except:
- self.particle_types = ()
- self.filename = hierarchy.filename
-
- def __repr__(self):
- return 'GadgetGrid_%05i'%self.id
-
-class GadgetHierarchy(GridGeometryHandler):
- grid = GadgetGrid
-
- def __init__(self, pf, data_style='gadget_hdf5'):
- self.filename = pf.filename
- self.directory = os.path.dirname(pf.filename)
- self.data_style = data_style
- self._handle = h5py.File(pf.filename)
- GridGeometryHandler.__init__(self, pf, data_style)
- self._handle.close()
- self._handle = None
-
-
- def _initialize_data_storage(self):
+class GadgetDomainFile(object):
+ def select(self, selector):
pass
- def _detect_fields(self):
- #this adds all the fields in
- #/particle_types/{Gas/Stars/etc.}/{position_x, etc.}
- self.field_list = []
- for ptype in self._handle['particle_types'].keys():
- for field in self._handle['particle_types'+'/'+ptype]:
- if field not in self.field_list:
- self.field_list += field,
-
- def _setup_classes(self):
- dd = self._get_data_reader_dict()
- GridGeometryHandler._setup_classes(self, dd)
- self.object_types.sort()
+ def count(self, selector):
+ pass
- def _count_grids(self):
- self.num_grids = len(self._handle['/grid_dimensions'])
-
- def _parse_hierarchy(self):
- f = self._handle # shortcut
- npa = np.array
- DLE = self.parameter_file.domain_left_edge
- DRE = self.parameter_file.domain_right_edge
- DW = (DRE - DLE)
-
- self.grid_levels.flat[:] = f['/grid_level'][:].astype("int32")
- LI = f['/grid_left_index'][:]
- print LI
- self.grid_dimensions[:] = f['/grid_dimensions'][:]
- self.grid_left_edge[:] = (LI * DW + DLE)
- dxs = 1.0/(2**(self.grid_levels+1)) * DW
- self.grid_right_edge[:] = self.grid_left_edge \
- + dxs *(1 + self.grid_dimensions)
- self.grid_particle_count.flat[:] = f['/grid_particle_count'][:].astype("int32")
- grid_parent_id = f['/grid_parent_id'][:]
- self.max_level = np.max(self.grid_levels)
-
- args = izip(xrange(self.num_grids), self.grid_levels.flat,
- grid_parent_id, LI,
- self.grid_dimensions, self.grid_particle_count.flat)
- self.grids = np.empty(len(args), dtype='object')
- for gi, (j,lvl,p, le, d, n) in enumerate(args):
- self.grids[gi] = self.grid(self,j,d,le,lvl,p,n)
-
- def _populate_grid_objects(self):
- for g in self.grids:
- if g._parent_id >= 0:
- parent = self.grids[g._parent_id]
- g.Parent = parent
- parent.Children.append(g)
- for g in self.grids:
- g._prepare_grid()
- g._setup_dx()
-
- def _setup_derived_fields(self):
- self.derived_field_list = []
-
-class GadgetStaticOutput(StaticOutput):
- _hierarchy_class = GadgetHierarchy
- _fieldinfo_fallback = GadgetFieldInfo
- _fieldinfo_known = KnownGadgetFields
-
- def __init__(self, filename,storage_filename=None) :
- self.storage_filename = storage_filename
- self.filename = filename
-
- StaticOutput.__init__(self, filename, 'gadget_infrastructure')
- self._set_units()
-
- def _set_units(self):
- self.units = {}
- self.time_units = {}
- self.time_units['1'] = 1
- self.units['1'] = 1.0
- self.units['cm'] = 1.0
- self.units['unitary'] = 1.0 / \
- (self.domain_right_edge - self.domain_left_edge).max()
- for unit in sec_conversion.keys():
- self.time_units[unit] = 1.0 / sec_conversion[unit]
-
- def _parse_parameter_file(self):
- fileh = h5py.File(self.filename)
- sim_param = fileh['/simulation_parameters'].attrs
- self.refine_by = sim_param['refine_by']
- self.dimensionality = sim_param['dimensionality']
- self.num_ghost_zones = sim_param['num_ghost_zones']
- self.field_ordering = sim_param['field_ordering']
- self.domain_dimensions = sim_param['domain_dimensions']
- self.current_time = sim_param['current_time']
- self.domain_left_edge = sim_param['domain_left_edge']
- self.domain_right_edge = sim_param['domain_right_edge']
- self.unique_identifier = sim_param['unique_identifier']
- self.cosmological_simulation = sim_param['cosmological_simulation']
- self.current_redshift = sim_param['current_redshift']
- self.omega_lambda = sim_param['omega_lambda']
- self.hubble_constant = sim_param['hubble_constant']
- fileh.close()
-
-
- @classmethod
- def _is_valid(self, *args, **kwargs):
- format = 'Gadget Infrastructure'
- add1 = 'griddded_data_format'
- add2 = 'data_software'
- try:
- fileh = h5py.File(args[0],'r')
- if add1 in fileh['/'].items():
- if add2 in fileh['/'+add1].attrs.keys():
- if fileh['/'+add1].attrs[add2] == format:
- fileh.close()
- return True
- fileh.close()
- except:
- pass
- return False
+class GadgetDomainSubset(object):
+ def __init__(self, domain, mask, cell_count):
+ self.mask = mask
+ self.domain = domain
+ self.oct_handler = domain.pf.h.oct_handler
+ self.cell_count = cell_count
+ level_counts = self.oct_handler.count_levels(
+ self.domain.pf.max_level, self.domain.domain_id, mask)
+ level_counts[1:] = level_counts[:-1]
+ level_counts[0] = 0
+ self.level_counts = np.add.accumulate(level_counts)
diff -r a1d1f0d4be060b1372732c33064e0515f0623db5 -r 3874ec822d9eda570abbcf277e85f2b7cbcd7175 yt/frontends/ramses/data_structures.py
--- a/yt/frontends/ramses/data_structures.py
+++ b/yt/frontends/ramses/data_structures.py
@@ -29,8 +29,6 @@
import cStringIO
from yt.funcs import *
-from yt.data_objects.grid_patch import \
- AMRGridPatch
from yt.geometry.oct_geometry_handler import \
OctreeGeometryHandler
from yt.geometry.geometry_handler import \
https://bitbucket.org/yt_analysis/yt-3.0/commits/16a2374a9a6b/
changeset: 16a2374a9a6b
branch: yt-3.0
user: MatthewTurk
date: 2012-12-07 23:19:24
summary: Have ParticleOctree collect domains on a per-oct basis, not per-particle.
Mandate only a single domain per Oct.
affected #: 2 files
diff -r 3874ec822d9eda570abbcf277e85f2b7cbcd7175 -r 16a2374a9a6bb093a891ad937d284b1c1b5d1690 yt/geometry/oct_container.pxd
--- a/yt/geometry/oct_container.pxd
+++ b/yt/geometry/oct_container.pxd
@@ -65,5 +65,4 @@
Oct *oct
ParticleArrays *next
np.float64_t **pos
- np.int64_t *domain_id
np.int64_t np
diff -r 3874ec822d9eda570abbcf277e85f2b7cbcd7175 -r 16a2374a9a6bb093a891ad937d284b1c1b5d1690 yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -311,25 +311,6 @@
count[cur.my_octs[i - cur.offset].domain - 1] += 1
return count
- def check(self, int curdom):
- cdef int dind, pi
- cdef Oct oct
- cdef OctAllocationContainer *cont = self.domains[curdom - 1]
- cdef int nbad = 0
- for pi in range(cont.n_assigned):
- oct = cont.my_octs[pi]
- for i in range(2):
- for j in range(2):
- for k in range(2):
- if oct.children[i][j][k] != NULL and \
- oct.children[i][j][k].level != oct.level + 1:
- if curdom == 61:
- print pi, oct.children[i][j][k].level,
- print oct.level
- nbad += 1
- print "DOMAIN % 3i HAS % 9i BAD OCTS (%s / %s / %s)" % (curdom, nbad,
- cont.n - cont.n_assigned, cont.n_assigned, cont.n)
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -580,7 +561,6 @@
for i in range(3):
free(o.sd.pos[i])
free(o.sd.pos)
- free(o.sd.domain_id)
free(o)
@cython.boundscheck(False)
@@ -689,6 +669,7 @@
malloc(sizeof(ParticleArrays))
cdef int i, j, k
my_oct.ind = my_oct.domain = -1
+ my_oct.domain = -1
my_oct.local_ind = self.nocts - 1
my_oct.pos[0] = my_oct.pos[1] = my_oct.pos[2] = -1
my_oct.level = -1
@@ -705,13 +686,11 @@
self.last_sd = sd
sd.oct = my_oct
sd.next = NULL
- sd.domain_id = <np.int64_t *> malloc(sizeof(np.int64_t) * 32)
sd.pos = <np.float64_t **> malloc(sizeof(np.float64_t*) * 3)
for i in range(3):
sd.pos[i] = <np.float64_t *> malloc(sizeof(np.float64_t) * 32)
for i in range(32):
sd.pos[0][i] = sd.pos[1][i] = sd.pos[2][i] = 0.0
- sd.domain_id[i] = -1
sd.np = 0
return my_oct
@@ -741,7 +720,7 @@
cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
if cur == NULL:
raise RuntimeError
- if self._check_refine(cur, cp) == 1:
+ if self._check_refine(cur, cp, domain_id) == 1:
self.refine_oct(cur, cp)
while cur.sd.np < 0:
for i in range(3):
@@ -754,23 +733,22 @@
cp[i] += dds[i]/2.0
cur = cur.children[ind[0]][ind[1]][ind[2]]
level += 1
- if self._check_refine(cur, cp) == 1:
+ if self._check_refine(cur, cp, domain_id) == 1:
self.refine_oct(cur, cp)
# Now we copy in our particle
pi = cur.sd.np
cur.level = level
for i in range(3):
cur.sd.pos[i][pi] = pp[i]
- cur.sd.domain_id[pi] = domain_id
+ cur.domain = domain_id
cur.sd.np += 1
- cdef int _check_refine(self, Oct *cur, np.float64_t cp[3]):
- cdef int mid = 16384
- cdef int mad = -16384
- for i in range(imax(cur.sd.np, 0)):
- mid = imin(cur.sd.domain_id[i], mid)
- mad = imax(cur.sd.domain_id[i], mad)
- if cur.sd.np == 32 or mid < mad:
+ cdef int _check_refine(self, Oct *cur, np.float64_t cp[3], int domain_id):
+ if cur.children[0][0][0] != NULL:
+ return 0
+ elif cur.sd.np == 32:
+ return 1
+ elif cur.domain >= 0 and cur.domain != domain_id:
return 1
return 0
@@ -787,7 +765,7 @@
noct.pos[2] = (o.pos[2] << 1) + k
noct.parent = o
o.children[i][j][k] = noct
- for m in range(32):
+ for m in range(o.sd.np):
for i in range(3):
if o.sd.pos[i][m] < pos[i]:
ind[i] = 0
@@ -797,12 +775,11 @@
k = noct.sd.np
for i in range(3):
noct.sd.pos[i][k] = o.sd.pos[i][m]
- noct.sd.domain_id[k] = o.sd.domain_id[k]
+ noct.domain = o.domain
noct.sd.np += 1
o.sd.np = -1
for i in range(3):
free(o.sd.pos[i])
- free(o.sd.domain_id)
free(o.sd.pos)
def recursively_count(self):
@@ -844,8 +821,7 @@
m = 1
break
if m == 0: continue
- for i in range(o.sd.np):
- dmask[o.sd.domain_id[i]] = 1
+ dmask[o.domain] = 1
return dmask.astype("bool")
@cython.boundscheck(False)
https://bitbucket.org/yt_analysis/yt-3.0/commits/60c1c1533c74/
changeset: 60c1c1533c74
branch: yt-3.0
user: MatthewTurk
date: 2012-12-07 23:41:03
summary: Use qsort to sort octs once they have been inserted into the particle octree
container's oct_list.
affected #: 1 file
diff -r 16a2374a9a6bb093a891ad937d284b1c1b5d1690 -r 60c1c1533c748cbbe190a280d95b42ba16f8026d yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -23,7 +23,7 @@
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
-from libc.stdlib cimport malloc, free
+from libc.stdlib cimport malloc, free, qsort
cimport numpy as np
import numpy as np
from oct_container cimport Oct, OctAllocationContainer, OctreeContainer
@@ -523,10 +523,21 @@
local_filled += 1
return local_filled
+cdef public int compare_octs(void *vo1, void *vo2) nogil:
+ cdef Oct *o1 = (<Oct**> vo1)[0]
+ cdef Oct *o2 = (<Oct**> vo2)[0]
+ if o1.domain < o2.domain: return -1
+ elif o1.domain == o2.domain:
+ if o1.level < o2.level: return -1
+ if o1.level > o2.level: return 1
+ else: return 0
+ elif o1.domain > o2.domain: return 1
+
cdef class ParticleOctreeContainer(OctreeContainer):
cdef ParticleArrays *first_sd
cdef ParticleArrays *last_sd
cdef Oct** oct_list
+ cdef np.int64_t *dom_offsets
def allocate_root(self):
cdef int i, j, k
@@ -548,6 +559,8 @@
for j in range(self.nn[1]):
for k in range(self.nn[2]):
self.visit_free(self.root_mesh[i][j][k])
+ free(self.oct_list)
+ free(self.dom_offsets)
cdef void visit_free(self, Oct *o):
cdef int i, j, k
@@ -649,10 +662,13 @@
def finalize(self):
self.oct_list = <Oct**> malloc(sizeof(Oct*)*self.nocts)
- cdef i = 0
+ cdef np.int64_t i = 0
cdef ParticleArrays *c = self.first_sd
while c != NULL:
self.oct_list[i] = c.oct
+ if c.oct.domain == -1:
+ assert(c.oct.children[0][0][0] == NULL)
+ assert(c.np == 0)
if c.np >= 0:
for j in range(3):
free(c.pos[j])
@@ -661,6 +677,15 @@
# We should also include a shortening of the domain IDs here
c = c.next
i += 1
+ assert(i == self.nocts)
+ qsort(self.oct_list, self.nocts, sizeof(Oct*), &compare_octs)
+ cdef int cur_dom = -1
+ self.dom_offsets = <np.int64_t *>malloc(sizeof(np.int64_t) *
+ self.max_domain + 1)
+ for i in range(self.nocts):
+ if self.oct_list[i].domain > cur_dom:
+ cur_dom = self.oct_list[i].domain
+ self.dom_offsets[cur_dom] = i
cdef Oct* allocate_oct(self):
self.nocts += 1
@@ -815,7 +840,7 @@
for oi in range(self.nocts):
m = 0
o = self.oct_list[oi]
- if o.sd.np <= 0: continue
+ if o.sd.np <= 0 or o.domain == -1: continue
for i in range(8):
if mask[oi, i] == 1:
m = 1
https://bitbucket.org/yt_analysis/yt-3.0/commits/fb30c2b133fe/
changeset: fb30c2b133fe
branch: yt-3.0
user: MatthewTurk
date: 2012-12-08 00:29:08
summary: Adding "count_levels" and making it work for particle octrees. Ready to set up
subsetting and IO.
affected #: 1 file
diff -r 60c1c1533c748cbbe190a280d95b42ba16f8026d -r fb30c2b133fe253bf631201d3ecb7c814e5a2ffd yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -666,26 +666,23 @@
cdef ParticleArrays *c = self.first_sd
while c != NULL:
self.oct_list[i] = c.oct
- if c.oct.domain == -1:
- assert(c.oct.children[0][0][0] == NULL)
- assert(c.np == 0)
if c.np >= 0:
for j in range(3):
free(c.pos[j])
free(c.pos)
c.pos = NULL
- # We should also include a shortening of the domain IDs here
c = c.next
i += 1
- assert(i == self.nocts)
qsort(self.oct_list, self.nocts, sizeof(Oct*), &compare_octs)
cdef int cur_dom = -1
self.dom_offsets = <np.int64_t *>malloc(sizeof(np.int64_t) *
- self.max_domain + 1)
+ self.max_domain + 2)
for i in range(self.nocts):
+ self.oct_list[i].local_ind = i
if self.oct_list[i].domain > cur_dom:
cur_dom = self.oct_list[i].domain
self.dom_offsets[cur_dom] = i
+ self.dom_offsets[cur_dom + 1] = self.nocts
cdef Oct* allocate_oct(self):
self.nocts += 1
@@ -727,6 +724,28 @@
c = c.next
return total
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def count_levels(self, int max_level, int domain_id,
+ np.ndarray[np.uint8_t, ndim=2, cast=True] mask):
+ cdef np.ndarray[np.int64_t, ndim=1] level_count
+ cdef Oct *o
+ cdef int oi, i
+ level_count = np.zeros(max_level+1, 'int64')
+ ndo
+ cdef np.int64_t ndo, doff
+ ndo = self.dom_offsets[domain_id + 1] \
+ - self.dom_offsets[domain_id]
+ doff = self.dom_offsets[domain_id]
+ print "Domain %s Total %s" % (domain_id, ndo)
+ for oi in range(ndo):
+ o = self.oct_list[oi + doff]
+ for i in range(8):
+ if mask[o.local_ind, i] == 0: continue
+ level_count[o.level] += 1
+ return level_count
+
def add(self, np.ndarray[np.float64_t, ndim=2] pos, np.int64_t domain_id):
cdef int no = pos.shape[0]
cdef int p, i, level
@@ -803,6 +822,7 @@
noct.domain = o.domain
noct.sd.np += 1
o.sd.np = -1
+ o.domain = -1
for i in range(3):
free(o.sd.pos[i])
free(o.sd.pos)
https://bitbucket.org/yt_analysis/yt-3.0/commits/5f02c6f858e6/
changeset: 5f02c6f858e6
branch: yt-3.0
user: MatthewTurk
date: 2012-12-18 18:47:05
summary: Merging from tip
affected #: 6 files
diff -r 58aeb5270b1589dfcee71c41517b3218e0a1e52d -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -347,8 +347,8 @@
else:
w = np.ones(chunk.size, dtype="float64")
icoords = chunk.icoords
- i1 = icoords[:,0]
- i2 = icoords[:,1]
+ i1 = icoords[:,x_dict[self.axis]]
+ i2 = icoords[:,y_dict[self.axis]]
ilevel = chunk.ires
tree.add_chunk_to_tree(i1, i2, ilevel, v, w)
diff -r 58aeb5270b1589dfcee71c41517b3218e0a1e52d -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -258,7 +258,7 @@
self._last_freq = field
self._last_finfo = self.field_info[fname]
return self._last_finfo
- raise YTFieldNotFound((fname, ftype), self)
+ raise YTFieldNotFound((ftype, fname), self)
def _reconstruct_pf(*args, **kwargs):
pfs = ParameterFileStore()
diff -r 58aeb5270b1589dfcee71c41517b3218e0a1e52d -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c yt/frontends/enzo/data_structures.py
--- a/yt/frontends/enzo/data_structures.py
+++ b/yt/frontends/enzo/data_structures.py
@@ -394,7 +394,7 @@
if last != g.filename:
if handle is not None: handle.close()
handle = h5py.File(g.filename)
- node = handle["/Grid%08i/ActiveParticles/" % g.id]
+ node = handle["/Grid%08i/Particles/" % g.id]
for ptype in (str(p) for p in node):
for field in (str(f) for f in node[ptype]):
_fields[ptype].append(field)
diff -r 58aeb5270b1589dfcee71c41517b3218e0a1e52d -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c yt/frontends/enzo/io.py
--- a/yt/frontends/enzo/io.py
+++ b/yt/frontends/enzo/io.py
@@ -59,7 +59,7 @@
size = 0
for chunk in chunks:
data = self._read_chunk_data(chunk, pfields, 'active',
- "/ActiveParticles/%s" % ptypes[0])
+ "/Particles/%s" % ptypes[0])
for g in chunk.objs:
if g.NumberOfActiveParticles == 0: continue
x, y, z = (data[g.id].pop(fn) for ft, fn in pfields)
@@ -73,7 +73,7 @@
ind = 0
for chunk in chunks:
data = self._read_chunk_data(chunk, read_fields, 'active',
- "/ActiveParticles/%s" % ptypes[0])
+ "/Particles/%s" % ptypes[0])
for g in chunk.objs:
if g.NumberOfActiveParticles == 0: continue
x, y, z = (data[g.id][fn] for ft, fn in pfields)
diff -r 58aeb5270b1589dfcee71c41517b3218e0a1e52d -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c yt/utilities/parallel_tools/parallel_analysis_interface.py
--- a/yt/utilities/parallel_tools/parallel_analysis_interface.py
+++ b/yt/utilities/parallel_tools/parallel_analysis_interface.py
@@ -435,7 +435,8 @@
to_share = {}
# If our objects object is slice-aware, like time series data objects are,
# this will prevent intermediate objects from being created.
- for obj_id, obj in enumerate(objects):
+ oiter = itertools.islice(enumerate(objects), my_new_id, None, njobs)
+ for obj_id, obj in oiter:
result_id = obj_id * njobs + my_new_id
if storage is not None:
rstore = ResultsStorage()
diff -r 58aeb5270b1589dfcee71c41517b3218e0a1e52d -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -1082,7 +1082,7 @@
_frb_generator = FixedResolutionBuffer
def __init__(self, pf, axis, fields, center='c', width=None, axes_unit=None,
- weight_field=None, max_level=None, origin='center-window', fontsize=15):
+ weight_field=None, origin='center-window', fontsize=15):
r"""Creates a projection plot from a parameter file
Given a pf object, an axis to project along, and a field name
@@ -1141,8 +1141,6 @@
the center of the plot window.
weight_field : string
The name of the weighting field. Set to None for no weight.
- max_level: int
- The maximum level to project to.
fontsize : integer
The size of the fonts for the axis, colorbar, and tick labels.
@@ -1163,7 +1161,7 @@
(bounds, center, units) = GetWindowParameters(axis, center, width, pf)
if axes_unit is None and units != ('1', '1'):
axes_unit = units
- proj = pf.h.proj(fields, axis, weight_field=weight_field,max_level=max_level,center=center)
+ proj = pf.h.proj(fields, axis, weight_field=weight_field, center=center)
PWViewerMPL.__init__(self,proj,bounds,origin=origin)
self.set_axes_unit(axes_unit)
https://bitbucket.org/yt_analysis/yt-3.0/commits/32e76c66e4e4/
changeset: 32e76c66e4e4
branch: yt-3.0
user: MatthewTurk
date: 2012-12-18 18:55:58
summary: Merging some particle octree work
affected #: 6 files
diff -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c -r 32e76c66e4e435ba4945c2d62e0390b957294f24 yt/frontends/gadget/api.py
--- a/yt/frontends/gadget/api.py
+++ b/yt/frontends/gadget/api.py
@@ -28,14 +28,3 @@
"""
-from .data_structures import \
- GadgetGrid, \
- GadgetHierarchy, \
- GadgetStaticOutput
-
-from .fields import \
- GadgetFieldInfo, \
- add_gadget_field
-
-from .io import \
- IOHandlerGadget
diff -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c -r 32e76c66e4e435ba4945c2d62e0390b957294f24 yt/frontends/gadget/data_structures.py
--- a/yt/frontends/gadget/data_structures.py
+++ b/yt/frontends/gadget/data_structures.py
@@ -1,13 +1,11 @@
"""
-Data structures for Gadget.
+Data structures for a generic SPH/Gadget frontend.
Author: Matthew Turk <matthewturk at gmail.com>
-Affiliation: UCSD
-Author: Chris Moody <cemoody at ucsc.edu>
-Affiliation: UCSC
+Affiliation: Columbia University
Homepage: http://yt-project.org/
License:
- Copyright (C) 2007-2011 Matthew Turk. All Rights Reserved.
+ Copyright (C) 2007-2012 Matthew Turk. All Rights Reserved.
This file is part of yt.
@@ -30,172 +28,31 @@
from itertools import izip
from yt.funcs import *
-from yt.data_objects.grid_patch import \
- AMRGridPatch
-from yt.geometry.grid_geometry_handler import \
- GridGeometryHandler
+from yt.geometry.oct_geometry_handler import \
+ OctreeGeometryHandler
+from yt.geometry.geometry_handler import \
+ GeometryHandler, YTDataChunk
from yt.data_objects.static_output import \
StaticOutput
-from yt.utilities.definitions import \
- sec_conversion
-from .fields import GadgetFieldInfo, KnownGadgetFields
from yt.data_objects.field_info_container import \
FieldInfoContainer, NullFunc
-class GadgetGrid(AMRGridPatch):
- _id_offset = 0
- def __init__(self, hierarchy, id, dimensions, start,
- level, parent_id, particle_count):
- AMRGridPatch.__init__(self, id, filename = hierarchy.filename,
- hierarchy = hierarchy)
- self.Parent = [] # Only one parent per grid
- self.Children = []
- self.Level = level
- self.ActiveDimensions = dimensions.copy()
- self.NumberOfParticles = particle_count
- self.start_index = start.copy().astype("int64")
- self.stop_index = self.start_index + dimensions.copy()
- self.id = id
- self._parent_id = parent_id
-
- try:
- padd = '/data/grid_%010i/particles' % id
- self.particle_types = self.hierarchy._handle[padd].keys()
- except:
- self.particle_types = ()
- self.filename = hierarchy.filename
-
- def __repr__(self):
- return 'GadgetGrid_%05i'%self.id
-
-class GadgetHierarchy(GridGeometryHandler):
- grid = GadgetGrid
-
- def __init__(self, pf, data_style='gadget_hdf5'):
- self.filename = pf.filename
- self.directory = os.path.dirname(pf.filename)
- self.data_style = data_style
- self._handle = h5py.File(pf.filename)
- GridGeometryHandler.__init__(self, pf, data_style)
- self._handle.close()
- self._handle = None
-
-
- def _initialize_data_storage(self):
+class GadgetDomainFile(object):
+ def select(self, selector):
pass
- def _detect_fields(self):
- #this adds all the fields in
- #/particle_types/{Gas/Stars/etc.}/{position_x, etc.}
- self.field_list = []
- for ptype in self._handle['particle_types'].keys():
- for field in self._handle['particle_types'+'/'+ptype]:
- if field not in self.field_list:
- self.field_list += field,
-
- def _setup_classes(self):
- dd = self._get_data_reader_dict()
- GridGeometryHandler._setup_classes(self, dd)
- self.object_types.sort()
+ def count(self, selector):
+ pass
- def _count_grids(self):
- self.num_grids = len(self._handle['/grid_dimensions'])
-
- def _parse_hierarchy(self):
- f = self._handle # shortcut
- npa = np.array
- DLE = self.parameter_file.domain_left_edge
- DRE = self.parameter_file.domain_right_edge
- DW = (DRE - DLE)
-
- self.grid_levels.flat[:] = f['/grid_level'][:].astype("int32")
- LI = f['/grid_left_index'][:]
- print LI
- self.grid_dimensions[:] = f['/grid_dimensions'][:]
- self.grid_left_edge[:] = (LI * DW + DLE)
- dxs = 1.0/(2**(self.grid_levels+1)) * DW
- self.grid_right_edge[:] = self.grid_left_edge \
- + dxs *(1 + self.grid_dimensions)
- self.grid_particle_count.flat[:] = f['/grid_particle_count'][:].astype("int32")
- grid_parent_id = f['/grid_parent_id'][:]
- self.max_level = np.max(self.grid_levels)
-
- args = izip(xrange(self.num_grids), self.grid_levels.flat,
- grid_parent_id, LI,
- self.grid_dimensions, self.grid_particle_count.flat)
- self.grids = np.empty(len(args), dtype='object')
- for gi, (j,lvl,p, le, d, n) in enumerate(args):
- self.grids[gi] = self.grid(self,j,d,le,lvl,p,n)
-
- def _populate_grid_objects(self):
- for g in self.grids:
- if g._parent_id >= 0:
- parent = self.grids[g._parent_id]
- g.Parent = parent
- parent.Children.append(g)
- for g in self.grids:
- g._prepare_grid()
- g._setup_dx()
-
- def _setup_derived_fields(self):
- self.derived_field_list = []
-
-class GadgetStaticOutput(StaticOutput):
- _hierarchy_class = GadgetHierarchy
- _fieldinfo_fallback = GadgetFieldInfo
- _fieldinfo_known = KnownGadgetFields
-
- def __init__(self, filename,storage_filename=None) :
- self.storage_filename = storage_filename
- self.filename = filename
-
- StaticOutput.__init__(self, filename, 'gadget_infrastructure')
- self._set_units()
-
- def _set_units(self):
- self.units = {}
- self.time_units = {}
- self.time_units['1'] = 1
- self.units['1'] = 1.0
- self.units['cm'] = 1.0
- self.units['unitary'] = 1.0 / \
- (self.domain_right_edge - self.domain_left_edge).max()
- for unit in sec_conversion.keys():
- self.time_units[unit] = 1.0 / sec_conversion[unit]
-
- def _parse_parameter_file(self):
- fileh = h5py.File(self.filename)
- sim_param = fileh['/simulation_parameters'].attrs
- self.refine_by = sim_param['refine_by']
- self.dimensionality = sim_param['dimensionality']
- self.num_ghost_zones = sim_param['num_ghost_zones']
- self.field_ordering = sim_param['field_ordering']
- self.domain_dimensions = sim_param['domain_dimensions']
- self.current_time = sim_param['current_time']
- self.domain_left_edge = sim_param['domain_left_edge']
- self.domain_right_edge = sim_param['domain_right_edge']
- self.unique_identifier = sim_param['unique_identifier']
- self.cosmological_simulation = sim_param['cosmological_simulation']
- self.current_redshift = sim_param['current_redshift']
- self.omega_lambda = sim_param['omega_lambda']
- self.hubble_constant = sim_param['hubble_constant']
- fileh.close()
-
-
- @classmethod
- def _is_valid(self, *args, **kwargs):
- format = 'Gadget Infrastructure'
- add1 = 'griddded_data_format'
- add2 = 'data_software'
- try:
- fileh = h5py.File(args[0],'r')
- if add1 in fileh['/'].items():
- if add2 in fileh['/'+add1].attrs.keys():
- if fileh['/'+add1].attrs[add2] == format:
- fileh.close()
- return True
- fileh.close()
- except:
- pass
- return False
+class GadgetDomainSubset(object):
+ def __init__(self, domain, mask, cell_count):
+ self.mask = mask
+ self.domain = domain
+ self.oct_handler = domain.pf.h.oct_handler
+ self.cell_count = cell_count
+ level_counts = self.oct_handler.count_levels(
+ self.domain.pf.max_level, self.domain.domain_id, mask)
+ level_counts[1:] = level_counts[:-1]
+ level_counts[0] = 0
+ self.level_counts = np.add.accumulate(level_counts)
diff -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c -r 32e76c66e4e435ba4945c2d62e0390b957294f24 yt/frontends/ramses/data_structures.py
--- a/yt/frontends/ramses/data_structures.py
+++ b/yt/frontends/ramses/data_structures.py
@@ -29,8 +29,6 @@
import cStringIO
from yt.funcs import *
-from yt.data_objects.grid_patch import \
- AMRGridPatch
from yt.geometry.oct_geometry_handler import \
OctreeGeometryHandler
from yt.geometry.geometry_handler import \
diff -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c -r 32e76c66e4e435ba4945c2d62e0390b957294f24 yt/geometry/oct_container.pxd
--- a/yt/geometry/oct_container.pxd
+++ b/yt/geometry/oct_container.pxd
@@ -24,6 +24,7 @@
"""
cimport numpy as np
+from fp_utils cimport *
cdef struct ParticleArrays
@@ -64,5 +65,4 @@
Oct *oct
ParticleArrays *next
np.float64_t **pos
- np.int64_t *domain_id
np.int64_t np
diff -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c -r 32e76c66e4e435ba4945c2d62e0390b957294f24 yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -23,7 +23,7 @@
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
-from libc.stdlib cimport malloc, free
+from libc.stdlib cimport malloc, free, qsort
cimport numpy as np
import numpy as np
from oct_container cimport Oct, OctAllocationContainer, OctreeContainer
@@ -311,25 +311,6 @@
count[cur.my_octs[i - cur.offset].domain - 1] += 1
return count
- def check(self, int curdom):
- cdef int dind, pi
- cdef Oct oct
- cdef OctAllocationContainer *cont = self.domains[curdom - 1]
- cdef int nbad = 0
- for pi in range(cont.n_assigned):
- oct = cont.my_octs[pi]
- for i in range(2):
- for j in range(2):
- for k in range(2):
- if oct.children[i][j][k] != NULL and \
- oct.children[i][j][k].level != oct.level + 1:
- if curdom == 61:
- print pi, oct.children[i][j][k].level,
- print oct.level
- nbad += 1
- print "DOMAIN % 3i HAS % 9i BAD OCTS (%s / %s / %s)" % (curdom, nbad,
- cont.n - cont.n_assigned, cont.n_assigned, cont.n)
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -542,10 +523,21 @@
local_filled += 1
return local_filled
+cdef public int compare_octs(void *vo1, void *vo2) nogil:
+ cdef Oct *o1 = (<Oct**> vo1)[0]
+ cdef Oct *o2 = (<Oct**> vo2)[0]
+ if o1.domain < o2.domain: return -1
+ elif o1.domain == o2.domain:
+ if o1.level < o2.level: return -1
+ if o1.level > o2.level: return 1
+ else: return 0
+ elif o1.domain > o2.domain: return 1
+
cdef class ParticleOctreeContainer(OctreeContainer):
cdef ParticleArrays *first_sd
cdef ParticleArrays *last_sd
cdef Oct** oct_list
+ cdef np.int64_t *dom_offsets
def allocate_root(self):
cdef int i, j, k
@@ -567,6 +559,8 @@
for j in range(self.nn[1]):
for k in range(self.nn[2]):
self.visit_free(self.root_mesh[i][j][k])
+ free(self.oct_list)
+ free(self.dom_offsets)
cdef void visit_free(self, Oct *o):
cdef int i, j, k
@@ -580,7 +574,6 @@
for i in range(3):
free(o.sd.pos[i])
free(o.sd.pos)
- free(o.sd.domain_id)
free(o)
@cython.boundscheck(False)
@@ -669,7 +662,7 @@
def finalize(self):
self.oct_list = <Oct**> malloc(sizeof(Oct*)*self.nocts)
- cdef i = 0
+ cdef np.int64_t i = 0
cdef ParticleArrays *c = self.first_sd
while c != NULL:
self.oct_list[i] = c.oct
@@ -678,9 +671,18 @@
free(c.pos[j])
free(c.pos)
c.pos = NULL
- # We should also include a shortening of the domain IDs here
c = c.next
i += 1
+ qsort(self.oct_list, self.nocts, sizeof(Oct*), &compare_octs)
+ cdef int cur_dom = -1
+ self.dom_offsets = <np.int64_t *>malloc(sizeof(np.int64_t) *
+ self.max_domain + 2)
+ for i in range(self.nocts):
+ self.oct_list[i].local_ind = i
+ if self.oct_list[i].domain > cur_dom:
+ cur_dom = self.oct_list[i].domain
+ self.dom_offsets[cur_dom] = i
+ self.dom_offsets[cur_dom + 1] = self.nocts
cdef Oct* allocate_oct(self):
self.nocts += 1
@@ -689,6 +691,7 @@
malloc(sizeof(ParticleArrays))
cdef int i, j, k
my_oct.ind = my_oct.domain = -1
+ my_oct.domain = -1
my_oct.local_ind = self.nocts - 1
my_oct.pos[0] = my_oct.pos[1] = my_oct.pos[2] = -1
my_oct.level = -1
@@ -705,13 +708,11 @@
self.last_sd = sd
sd.oct = my_oct
sd.next = NULL
- sd.domain_id = <np.int64_t *> malloc(sizeof(np.int64_t) * 32)
sd.pos = <np.float64_t **> malloc(sizeof(np.float64_t*) * 3)
for i in range(3):
sd.pos[i] = <np.float64_t *> malloc(sizeof(np.float64_t) * 32)
for i in range(32):
sd.pos[0][i] = sd.pos[1][i] = sd.pos[2][i] = 0.0
- sd.domain_id[i] = -1
sd.np = 0
return my_oct
@@ -723,12 +724,35 @@
c = c.next
return total
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def count_levels(self, int max_level, int domain_id,
+ np.ndarray[np.uint8_t, ndim=2, cast=True] mask):
+ cdef np.ndarray[np.int64_t, ndim=1] level_count
+ cdef Oct *o
+ cdef int oi, i
+ level_count = np.zeros(max_level+1, 'int64')
+ ndo
+ cdef np.int64_t ndo, doff
+ ndo = self.dom_offsets[domain_id + 1] \
+ - self.dom_offsets[domain_id]
+ doff = self.dom_offsets[domain_id]
+ print "Domain %s Total %s" % (domain_id, ndo)
+ for oi in range(ndo):
+ o = self.oct_list[oi + doff]
+ for i in range(8):
+ if mask[o.local_ind, i] == 0: continue
+ level_count[o.level] += 1
+ return level_count
+
def add(self, np.ndarray[np.float64_t, ndim=2] pos, np.int64_t domain_id):
cdef int no = pos.shape[0]
cdef int p, i, level
cdef np.float64_t dds[3], cp[3], pp[3]
cdef int ind[3]
self.max_domain = max(self.max_domain, domain_id)
+ cdef int mid, mad
if self.root_mesh[0][0][0] == NULL: self.allocate_root()
for p in range(no):
level = 0
@@ -740,7 +764,7 @@
cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
if cur == NULL:
raise RuntimeError
- if cur.sd.np == 32:
+ if self._check_refine(cur, cp, domain_id) == 1:
self.refine_oct(cur, cp)
while cur.sd.np < 0:
for i in range(3):
@@ -753,16 +777,25 @@
cp[i] += dds[i]/2.0
cur = cur.children[ind[0]][ind[1]][ind[2]]
level += 1
- if cur.sd.np == 32:
+ if self._check_refine(cur, cp, domain_id) == 1:
self.refine_oct(cur, cp)
# Now we copy in our particle
pi = cur.sd.np
cur.level = level
for i in range(3):
cur.sd.pos[i][pi] = pp[i]
- cur.sd.domain_id[pi] = domain_id
+ cur.domain = domain_id
cur.sd.np += 1
+ cdef int _check_refine(self, Oct *cur, np.float64_t cp[3], int domain_id):
+ if cur.children[0][0][0] != NULL:
+ return 0
+ elif cur.sd.np == 32:
+ return 1
+ elif cur.domain >= 0 and cur.domain != domain_id:
+ return 1
+ return 0
+
cdef void refine_oct(self, Oct *o, np.float64_t pos[3]):
cdef int i, j, k, m, ind[3]
cdef Oct *noct
@@ -776,7 +809,7 @@
noct.pos[2] = (o.pos[2] << 1) + k
noct.parent = o
o.children[i][j][k] = noct
- for m in range(32):
+ for m in range(o.sd.np):
for i in range(3):
if o.sd.pos[i][m] < pos[i]:
ind[i] = 0
@@ -786,12 +819,12 @@
k = noct.sd.np
for i in range(3):
noct.sd.pos[i][k] = o.sd.pos[i][m]
- noct.sd.domain_id[k] = o.sd.domain_id[k]
+ noct.domain = o.domain
noct.sd.np += 1
o.sd.np = -1
+ o.domain = -1
for i in range(3):
free(o.sd.pos[i])
- free(o.sd.domain_id)
free(o.sd.pos)
def recursively_count(self):
@@ -827,14 +860,13 @@
for oi in range(self.nocts):
m = 0
o = self.oct_list[oi]
- if o.sd.np <= 0: continue
+ if o.sd.np <= 0 or o.domain == -1: continue
for i in range(8):
if mask[oi, i] == 1:
m = 1
break
if m == 0: continue
- for i in range(o.sd.np):
- dmask[o.sd.domain_id[i]] = 1
+ dmask[o.domain] = 1
return dmask.astype("bool")
@cython.boundscheck(False)
diff -r 5f02c6f858e6fa0b99ebc3c72e56f44424cafc1c -r 32e76c66e4e435ba4945c2d62e0390b957294f24 yt/geometry/setup.py
--- a/yt/geometry/setup.py
+++ b/yt/geometry/setup.py
@@ -9,6 +9,7 @@
config = Configuration('geometry',parent_package,top_path)
config.add_extension("oct_container",
["yt/geometry/oct_container.pyx"],
+ include_dirs=["yt/utilities/lib/"],
libraries=["m"],
depends=["yt/utilities/lib/fp_utils.pxd",
"yt/geometry/oct_container.pxd",
https://bitbucket.org/yt_analysis/yt-3.0/commits/9678327f0b7d/
changeset: 9678327f0b7d
branch: yt-3.0
user: MatthewTurk
date: 2012-12-21 17:04:08
summary: Moving 'gadget' frontend to 'sph,' where it will be made more generic.
affected #: 16 files
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/gadget/__init__.py
--- a/yt/frontends/gadget/__init__.py
+++ /dev/null
@@ -1,29 +0,0 @@
-"""
-API for yt.frontends.gadget
-
-Author: Matthew Turk <matthewturk at gmail.com>
-Affiliation: UCSD
-Author: J.S. Oishi <jsoishi at gmail.com>
-Affiliation: KIPAC/SLAC/Stanford
-Author: Britton Smith <brittonsmith at gmail.com>
-Affiliation: MSU
-Homepage: http://yt-project.org/
-License:
- Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
-
- This file is part of yt.
-
- yt is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-"""
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/gadget/api.py
--- a/yt/frontends/gadget/api.py
+++ /dev/null
@@ -1,30 +0,0 @@
-"""
-API for yt.frontends.gadget
-
-Author: Matthew Turk <matthewturk at gmail.com>
-Affiliation: UCSD
-Author: J.S. Oishi <jsoishi at gmail.com>
-Affiliation: KIPAC/SLAC/Stanford
-Author: Britton Smith <brittonsmith at gmail.com>
-Affiliation: MSU
-Homepage: http://yt-project.org/
-License:
- Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
-
- This file is part of yt.
-
- yt is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-"""
-
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/gadget/data_structures.py
--- a/yt/frontends/gadget/data_structures.py
+++ /dev/null
@@ -1,58 +0,0 @@
-"""
-Data structures for a generic SPH/Gadget frontend.
-
-Author: Matthew Turk <matthewturk at gmail.com>
-Affiliation: Columbia University
-Homepage: http://yt-project.org/
-License:
- Copyright (C) 2007-2012 Matthew Turk. All Rights Reserved.
-
- This file is part of yt.
-
- yt is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-"""
-
-import h5py
-import numpy as np
-from itertools import izip
-
-from yt.funcs import *
-from yt.geometry.oct_geometry_handler import \
- OctreeGeometryHandler
-from yt.geometry.geometry_handler import \
- GeometryHandler, YTDataChunk
-from yt.data_objects.static_output import \
- StaticOutput
-
-from yt.data_objects.field_info_container import \
- FieldInfoContainer, NullFunc
-
-class GadgetDomainFile(object):
- def select(self, selector):
- pass
-
- def count(self, selector):
- pass
-
-class GadgetDomainSubset(object):
- def __init__(self, domain, mask, cell_count):
- self.mask = mask
- self.domain = domain
- self.oct_handler = domain.pf.h.oct_handler
- self.cell_count = cell_count
- level_counts = self.oct_handler.count_levels(
- self.domain.pf.max_level, self.domain.domain_id, mask)
- level_counts[1:] = level_counts[:-1]
- level_counts[0] = 0
- self.level_counts = np.add.accumulate(level_counts)
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/gadget/fields.py
--- a/yt/frontends/gadget/fields.py
+++ /dev/null
@@ -1,133 +0,0 @@
-"""
-Gadget-specific fields
-
-Author: Christopher E Moody <juxtaposicion at gmail.com>
-Affiliation: UC Santa Cruz
-Homepage: http://yt-project.org/
-License:
- Copyright (C) 2010-2011 Christopher E Moody, Matthew Turk. All Rights Reserved.
-
- This file is part of yt.
-
- yt is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-"""
-
-import numpy as np
-
-from yt.funcs import *
-from yt.data_objects.field_info_container import \
- FieldInfoContainer, \
- FieldInfo, \
- ValidateParameter, \
- ValidateDataField, \
- ValidateProperty, \
- ValidateSpatial, \
- ValidateGridType
-import yt.data_objects.universal_fields
-
-GadgetFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
-add_gadget_field = GadgetFieldInfo.add_field
-
-add_field = add_gadget_field
-
-translation_dict = {"particle_position_x" : "position_x",
- "particle_position_y" : "position_y",
- "particle_position_z" : "position_z",
- }
-
-def _generate_translation(mine, theirs):
- pfield = mine.startswith("particle")
- add_field(theirs, function=lambda a, b: b[mine], take_log=True,
- particle_type = pfield)
-
-for f,v in translation_dict.items():
- if v not in GadgetFieldInfo:
- # Note here that it's the yt field that we check for particle nature
- pfield = f.startswith("particle")
- add_field(v, function=lambda a,b: None, take_log=False,
- validators = [ValidateDataField(v)],
- particle_type = pfield)
- print "Setting up translator from %s to %s" % (v, f)
- _generate_translation(v, f)
-
-
-#for f,v in translation_dict.items():
-# add_field(f, function=lambda a,b: None, take_log=True,
-# validators = [ValidateDataField(v)],
-# units=r"\rm{cm}")
-# add_field(v, function=lambda a,b: None, take_log=True,
-# validators = [ValidateDataField(v)],
-# units=r"\rm{cm}")
-
-
-
-add_field("position_x", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("position_x")],
- particle_type = True,
- units=r"\rm{cm}")
-
-add_field("position_y", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("position_y")],
- particle_type = True,
- units=r"\rm{cm}")
-
-add_field("position_z", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("position_z")],
- particle_type = True,
- units=r"\rm{cm}")
-
-add_field("VEL", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("VEL")],
- units=r"")
-
-add_field("id", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("ID")],
- units=r"")
-
-add_field("mass", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("mass")],
- units=r"\rm{g}")
-def _particle_mass(field, data):
- return data["mass"]/just_one(data["CellVolume"])
-def _convert_particle_mass(data):
- return 1.0
-add_field("particle_mass", function=_particle_mass, take_log=True,
- convert_function=_convert_particle_mass,
- validators = [ValidateSpatial(0)],
- units=r"\mathrm{g}\/\mathrm{cm}^{-3}")
-
-add_field("U", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("U")],
- units=r"")
-
-add_field("NE", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("NE")],
- units=r"")
-
-add_field("POT", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("POT")],
- units=r"")
-
-add_field("ACCE", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("ACCE")],
- units=r"")
-
-add_field("ENDT", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("ENDT")],
- units=r"")
-
-add_field("TSTP", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("TSTP")],
- units=r"")
-
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/gadget/io.py
--- a/yt/frontends/gadget/io.py
+++ /dev/null
@@ -1,55 +0,0 @@
-"""
-Gadget-specific data-file handling function
-
-Author: Christopher E Moody <juxtaposicion at gmail.com>
-Affiliation: UC Santa Cruz
-Homepage: http://yt-project.org/
-License:
- Copyright (C) 2010-2011 Christopher E Moody, Matthew Turk. All Rights Reserved.
-
- This file is part of yt.
-
- yt is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-"""
-
-import h5py
-import numpy as np
-
-from yt.utilities.io_handler import \
- BaseIOHandler
-
-class IOHandlerGadget(BaseIOHandler):
- _data_style = 'gadget_infrastructure'
- def _read_data_set(self, grid, field):
- data = []
- fh = h5py.File(grid.filename,mode='r')
- for ptype in grid.particle_types:
- address = '/data/grid_%010i/particles/%s/%s' % (grid.id, ptype, field)
- data.append(fh[address][:])
- if len(data) > 0:
- data = np.concatenate(data)
- fh.close()
- return np.array(data)
- def _read_field_names(self,grid):
- adr = grid.Address
- fh = h5py.File(grid.filename,mode='r')
- rets = cPickle.loads(fh['/root'].attrs['fieldnames'])
- fh.close()
- return rets
-
- def _read_data_slice(self,grid, field, axis, coord):
- #how would we implement axis here?
- dat = self._read_data_set(grid,field)
- return dat[coord]
-
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/gadget/setup.py
--- a/yt/frontends/gadget/setup.py
+++ /dev/null
@@ -1,13 +0,0 @@
-#!/usr/bin/env python
-import setuptools
-import os
-import sys
-import os.path
-
-
-def configuration(parent_package='', top_path=None):
- from numpy.distutils.misc_util import Configuration
- config = Configuration('gadget', parent_package, top_path)
- config.make_config_py() # installs __config__.py
- #config.make_svn_version_py()
- return config
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/sph/gadget/__init__.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/__init__.py
@@ -0,0 +1,29 @@
+"""
+API for yt.frontends.gadget
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Author: J.S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: MSU
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/sph/gadget/api.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/api.py
@@ -0,0 +1,30 @@
+"""
+API for yt.frontends.gadget
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Author: J.S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: MSU
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/sph/gadget/data_structures.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/data_structures.py
@@ -0,0 +1,58 @@
+"""
+Data structures for a generic SPH/Gadget frontend.
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2007-2012 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import h5py
+import numpy as np
+from itertools import izip
+
+from yt.funcs import *
+from yt.geometry.oct_geometry_handler import \
+ OctreeGeometryHandler
+from yt.geometry.geometry_handler import \
+ GeometryHandler, YTDataChunk
+from yt.data_objects.static_output import \
+ StaticOutput
+
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, NullFunc
+
+class GadgetDomainFile(object):
+ def select(self, selector):
+ pass
+
+ def count(self, selector):
+ pass
+
+class GadgetDomainSubset(object):
+ def __init__(self, domain, mask, cell_count):
+ self.mask = mask
+ self.domain = domain
+ self.oct_handler = domain.pf.h.oct_handler
+ self.cell_count = cell_count
+ level_counts = self.oct_handler.count_levels(
+ self.domain.pf.max_level, self.domain.domain_id, mask)
+ level_counts[1:] = level_counts[:-1]
+ level_counts[0] = 0
+ self.level_counts = np.add.accumulate(level_counts)
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/sph/gadget/fields.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/fields.py
@@ -0,0 +1,133 @@
+"""
+Gadget-specific fields
+
+Author: Christopher E Moody <juxtaposicion at gmail.com>
+Affiliation: UC Santa Cruz
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Christopher E Moody, Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+
+from yt.funcs import *
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, \
+ FieldInfo, \
+ ValidateParameter, \
+ ValidateDataField, \
+ ValidateProperty, \
+ ValidateSpatial, \
+ ValidateGridType
+import yt.data_objects.universal_fields
+
+GadgetFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
+add_gadget_field = GadgetFieldInfo.add_field
+
+add_field = add_gadget_field
+
+translation_dict = {"particle_position_x" : "position_x",
+ "particle_position_y" : "position_y",
+ "particle_position_z" : "position_z",
+ }
+
+def _generate_translation(mine, theirs):
+ pfield = mine.startswith("particle")
+ add_field(theirs, function=lambda a, b: b[mine], take_log=True,
+ particle_type = pfield)
+
+for f,v in translation_dict.items():
+ if v not in GadgetFieldInfo:
+ # Note here that it's the yt field that we check for particle nature
+ pfield = f.startswith("particle")
+ add_field(v, function=lambda a,b: None, take_log=False,
+ validators = [ValidateDataField(v)],
+ particle_type = pfield)
+ print "Setting up translator from %s to %s" % (v, f)
+ _generate_translation(v, f)
+
+
+#for f,v in translation_dict.items():
+# add_field(f, function=lambda a,b: None, take_log=True,
+# validators = [ValidateDataField(v)],
+# units=r"\rm{cm}")
+# add_field(v, function=lambda a,b: None, take_log=True,
+# validators = [ValidateDataField(v)],
+# units=r"\rm{cm}")
+
+
+
+add_field("position_x", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("position_x")],
+ particle_type = True,
+ units=r"\rm{cm}")
+
+add_field("position_y", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("position_y")],
+ particle_type = True,
+ units=r"\rm{cm}")
+
+add_field("position_z", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("position_z")],
+ particle_type = True,
+ units=r"\rm{cm}")
+
+add_field("VEL", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("VEL")],
+ units=r"")
+
+add_field("id", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("ID")],
+ units=r"")
+
+add_field("mass", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("mass")],
+ units=r"\rm{g}")
+def _particle_mass(field, data):
+ return data["mass"]/just_one(data["CellVolume"])
+def _convert_particle_mass(data):
+ return 1.0
+add_field("particle_mass", function=_particle_mass, take_log=True,
+ convert_function=_convert_particle_mass,
+ validators = [ValidateSpatial(0)],
+ units=r"\mathrm{g}\/\mathrm{cm}^{-3}")
+
+add_field("U", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("U")],
+ units=r"")
+
+add_field("NE", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("NE")],
+ units=r"")
+
+add_field("POT", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("POT")],
+ units=r"")
+
+add_field("ACCE", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("ACCE")],
+ units=r"")
+
+add_field("ENDT", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("ENDT")],
+ units=r"")
+
+add_field("TSTP", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("TSTP")],
+ units=r"")
+
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/sph/gadget/io.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/io.py
@@ -0,0 +1,55 @@
+"""
+Gadget-specific data-file handling function
+
+Author: Christopher E Moody <juxtaposicion at gmail.com>
+Affiliation: UC Santa Cruz
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Christopher E Moody, Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import h5py
+import numpy as np
+
+from yt.utilities.io_handler import \
+ BaseIOHandler
+
+class IOHandlerGadget(BaseIOHandler):
+ _data_style = 'gadget_infrastructure'
+ def _read_data_set(self, grid, field):
+ data = []
+ fh = h5py.File(grid.filename,mode='r')
+ for ptype in grid.particle_types:
+ address = '/data/grid_%010i/particles/%s/%s' % (grid.id, ptype, field)
+ data.append(fh[address][:])
+ if len(data) > 0:
+ data = np.concatenate(data)
+ fh.close()
+ return np.array(data)
+ def _read_field_names(self,grid):
+ adr = grid.Address
+ fh = h5py.File(grid.filename,mode='r')
+ rets = cPickle.loads(fh['/root'].attrs['fieldnames'])
+ fh.close()
+ return rets
+
+ def _read_data_slice(self,grid, field, axis, coord):
+ #how would we implement axis here?
+ dat = self._read_data_set(grid,field)
+ return dat[coord]
+
diff -r 32e76c66e4e435ba4945c2d62e0390b957294f24 -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 yt/frontends/sph/gadget/setup.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/setup.py
@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+import setuptools
+import os
+import sys
+import os.path
+
+
+def configuration(parent_package='', top_path=None):
+ from numpy.distutils.misc_util import Configuration
+ config = Configuration('gadget', parent_package, top_path)
+ config.make_config_py() # installs __config__.py
+ #config.make_svn_version_py()
+ return config
https://bitbucket.org/yt_analysis/yt-3.0/commits/a91f4af64570/
changeset: a91f4af64570
branch: yt-3.0
user: MatthewTurk
date: 2012-12-21 19:07:13
summary: Initial import of the 'sph' frontend, which will include a few different types
of SPH output.
affected #: 9 files
diff -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 -r a91f4af645706634a602ef6ef36768e0ea48ce71 yt/frontends/gdf/data_structures.py
--- a/yt/frontends/gdf/data_structures.py
+++ b/yt/frontends/gdf/data_structures.py
@@ -247,7 +247,9 @@
try:
fileh = h5py.File(args[0],'r')
if "gridded_data_format" in fileh:
+ fileh.close()
return True
+ fileh.close()
except:
pass
return False
diff -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 -r a91f4af645706634a602ef6ef36768e0ea48ce71 yt/frontends/sph/__init__.py
--- /dev/null
+++ b/yt/frontends/sph/__init__.py
@@ -0,0 +1,29 @@
+"""
+API for yt.frontends.gadget
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Author: J.S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: MSU
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
diff -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 -r a91f4af645706634a602ef6ef36768e0ea48ce71 yt/frontends/sph/api.py
--- /dev/null
+++ b/yt/frontends/sph/api.py
@@ -0,0 +1,35 @@
+"""
+API for yt.frontends.sph
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Author: J.S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: MSU
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+
+from .data_structures import \
+ OWLSStaticOutput
+
+from .io import \
+ IOHandlerOWLS
diff -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 -r a91f4af645706634a602ef6ef36768e0ea48ce71 yt/frontends/sph/data_structures.py
--- /dev/null
+++ b/yt/frontends/sph/data_structures.py
@@ -0,0 +1,107 @@
+"""
+Data structures for a generic SPH/Gadget frontend.
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2007-2012 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import h5py
+import numpy as np
+import stat
+from itertools import izip
+
+from yt.funcs import *
+from yt.geometry.oct_geometry_handler import \
+ OctreeGeometryHandler
+from yt.geometry.geometry_handler import \
+ GeometryHandler, YTDataChunk
+from yt.data_objects.static_output import \
+ StaticOutput
+from yt.utilities.definitions import \
+ mpc_conversion, sec_conversion
+from .fields import \
+ OWLSFieldInfo, \
+ KnownOWLSFields
+
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, NullFunc
+
+class ParticleDomainFile(object):
+ def select(self, selector):
+ pass
+
+ def count(self, selector):
+ pass
+
+class ParticleDomainSubset(object):
+ def __init__(self, domain, mask, cell_count):
+ pass
+
+class OWLSStaticOutput(StaticOutput):
+ _fieldinfo_fallback = OWLSFieldInfo
+ _fieldinfo_known = KnownOWLSFields
+
+ def __init__(self, filename, data_style="OWLS", root_dimensions = 64):
+ self._root_dimensions = root_dimensions
+ super(OWLSStaticOutput, self).__init__(filename, data_style)
+
+ def __repr__(self):
+ return os.path.basename(self.parameter_filename).split(".")[0]
+
+ def _set_units(self):
+ pass
+
+ def _parse_parameter_file(self):
+ handle = h5py.File(self.parameter_filename)
+ hvals = {}
+ hvals.update(handle["/Header"].attrs)
+
+ self.dimensionality = 3
+ self.refine_by = 2
+ self.parameters["HydroMethod"] = "sph"
+ self.unique_identifier = \
+ int(os.stat(self.parameter_filename)[stat.ST_CTIME])
+ # Set standard values
+ self.current_time = hvals["Time_GYR"] * sec_conversion["Gyr"]
+ self.domain_left_edge = np.zeros(3, "float64")
+ self.domain_right_edge = np.ones(3, "float64") * hvals["BoxSize"]
+ self.domain_dimensions = np.ones(3, "int32") * self._root_dimensions
+ self.cosmological_simulation = 1
+ self.current_redshift = hvals["Redshift"]
+ self.omega_lambda = hvals["OmegaLambda"]
+ self.omega_matter = hvals["Omega0"]
+ self.hubble_constant = hvals["HubbleParam"]
+ self.parameters = hvals
+
+ handle.close()
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ try:
+ fileh = h5py.File(args[0],'r')
+ if "Constants" in fileh["/"].keys() and \
+ "Header" in fileh["/"].keys():
+ fileh.close()
+ return True
+ fileh.close()
+ except:
+ pass
+ return False
diff -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 -r a91f4af645706634a602ef6ef36768e0ea48ce71 yt/frontends/sph/fields.py
--- /dev/null
+++ b/yt/frontends/sph/fields.py
@@ -0,0 +1,44 @@
+"""
+OWLS-specific fields
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2012 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+
+from yt.funcs import *
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, \
+ FieldInfo, \
+ ValidateParameter, \
+ ValidateDataField, \
+ ValidateProperty, \
+ ValidateSpatial, \
+ ValidateGridType
+import yt.data_objects.universal_fields
+
+OWLSFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
+add_owls_field = OWLSFieldInfo.add_field
+
+KnownOWLSFields = FieldInfoContainer()
+add_OWLS_field = KnownOWLSFields.add_field
+
diff -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 -r a91f4af645706634a602ef6ef36768e0ea48ce71 yt/frontends/sph/io.py
--- /dev/null
+++ b/yt/frontends/sph/io.py
@@ -0,0 +1,74 @@
+"""
+Gadget-specific data-file handling function
+
+Author: Christopher E Moody <juxtaposicion at gmail.com>
+Affiliation: UC Santa Cruz
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Christopher E Moody, Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import h5py
+import numpy as np
+from yt.funcs import *
+
+from yt.utilities.io_handler import \
+ BaseIOHandler
+
+class IOHandlerOWLS(BaseIOHandler):
+ _data_style = "OWLS"
+
+ def _read_fluid_selection(self, chunks, selector, fields, size):
+ raise NotImplementedError
+
+ def _read_particle_selection_by_type(self, chunks, selector, fields):
+ rv = {}
+ # We first need a set of masks for each particle type
+ ptf = defaultdict(lambda: list)
+ psize = defaultdict(lambda: 0)
+ for ftype, fname in fields:
+ ptf[ftype].append(fname)
+ for chunk in chunks: # Will be OWLS domains
+ for subset in chunk.objs:
+ for ptype, field_list in sorted(ptf.items()):
+ f = h5py.File(subset.domain.filename, "r")
+ coords = f["/PartType%s/Coordinates" % ptype][:]
+ psize[ptype] += subset.count_particles(selector,
+ coords[:,0], coords[:,1], coords[:,2])
+ del coords
+ # Now we have all the sizes, and we can allocate
+ ind = {}
+ for field in fields:
+ rv[field] = np.empty(size, dtype="float64")
+ ind[field] = 0
+ for chunk in chunks: # Will be OWLS domains
+ for subset in chunk.objs:
+ for ptype, field_list in sorted(ptf.items()):
+ f = h5py.File(subset.domain.filename, "r")
+ g = f["/PartType%s" % ptype]
+ coords = g["Coordinates"][:]
+ mask = subset.select_particles(selector,
+ coords[:,0], coords[:,1], coords[:,2])
+ del coords
+ if mask is None: continue
+ for field in field_list:
+ data = g[field][mask,...]
+ my_ind = ind[ptype, field]
+ rv[ptype, field][my_ind:my_ind + data.size] = data
+ ind[ptype, field] += data.shape[0]
+ return rv
diff -r 9678327f0b7d467c4895746ddfc47480fa9b8c63 -r a91f4af645706634a602ef6ef36768e0ea48ce71 yt/frontends/sph/setup.py
--- /dev/null
+++ b/yt/frontends/sph/setup.py
@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+import setuptools
+import os
+import sys
+import os.path
+
+
+def configuration(parent_package='', top_path=None):
+ from numpy.distutils.misc_util import Configuration
+ config = Configuration('gadget', parent_package, top_path)
+ config.make_config_py() # installs __config__.py
+ #config.make_svn_version_py()
+ return config
https://bitbucket.org/yt_analysis/yt-3.0/commits/603c6241559a/
changeset: 603c6241559a
branch: yt-3.0
user: MatthewTurk
date: 2012-12-21 22:11:08
summary: Initial import of OWLS frontend. Reads particles, but no known particle fields
yet. No density estimator or fluid deposition.
affected #: 5 files
diff -r a91f4af645706634a602ef6ef36768e0ea48ce71 -r 603c6241559a80611ce716c3cb8d588d57175582 yt/frontends/sph/data_structures.py
--- a/yt/frontends/sph/data_structures.py
+++ b/yt/frontends/sph/data_structures.py
@@ -26,11 +26,14 @@
import h5py
import numpy as np
import stat
+import weakref
from itertools import izip
from yt.funcs import *
from yt.geometry.oct_geometry_handler import \
OctreeGeometryHandler
+from yt.geometry.oct_container import \
+ ParticleOctreeContainer
from yt.geometry.geometry_handler import \
GeometryHandler, YTDataChunk
from yt.data_objects.static_output import \
@@ -45,6 +48,13 @@
FieldInfoContainer, NullFunc
class ParticleDomainFile(object):
+ def __init__(self, pf, io, domain_filename, domain_number):
+ self.pf = pf
+ self.io = weakref.proxy(io)
+ self.domain_filename = domain_filename
+ self.domain_number = domain_number
+ self.total_particles = self.io._count_particles(domain_filename)
+
def select(self, selector):
pass
@@ -52,22 +62,122 @@
pass
class ParticleDomainSubset(object):
- def __init__(self, domain, mask, cell_count):
- pass
+ def __init__(self, domain, mask, count):
+ self.domain = domain
+ self.mask = mask
+ self.cell_count = count
+ self.oct_handler = domain.pf.h.oct_handler
+
+ def icoords(self, dobj):
+ return self.oct_handler.icoords(self.domain.domain_id, self.mask,
+ self.cell_count,
+ self.level_counts.copy())
+
+ def fcoords(self, dobj):
+ return self.oct_handler.fcoords(self.domain.domain_id, self.mask,
+ self.cell_count,
+ self.level_counts.copy())
+
+ def fwidth(self, dobj):
+ # Recall domain_dimensions is the number of cells, not octs
+ base_dx = 1.0/self.domain.pf.domain_dimensions
+ widths = np.empty((self.cell_count, 3), dtype="float64")
+ dds = (2**self.ires(dobj))
+ for i in range(3):
+ widths[:,i] = base_dx[i] / dds
+ return widths
+
+ def ires(self, dobj):
+ return self.oct_handler.ires(self.domain.domain_id, self.mask,
+ self.cell_count,
+ self.level_counts.copy())
+
+
+class ParticleGeometryHandler(OctreeGeometryHandler):
+
+ def __init__(self, pf, data_style):
+ self.data_style = data_style
+ self.parameter_file = weakref.proxy(pf)
+ # for now, the hierarchy file is the parameter file!
+ self.hierarchy_filename = self.parameter_file.parameter_filename
+ self.directory = os.path.dirname(self.hierarchy_filename)
+ self.float_type = np.float64
+ super(ParticleGeometryHandler, self).__init__(pf, data_style)
+
+ def _initialize_oct_handler(self):
+ template = self.parameter_file.domain_template
+ ndoms = self.parameter_file.domain_count
+ self.domains = [ParticleDomainFile(self.parameter_file, self.io, template % i, i)
+ for i in range(ndoms)]
+ total_particles = sum(d.total_particles for d in self.domains)
+ self.oct_handler = ParticleOctreeContainer(
+ self.parameter_file.domain_dimensions,
+ self.parameter_file.domain_left_edge,
+ self.parameter_file.domain_right_edge)
+ mylog.info("Allocating for %0.3e particles", total_particles)
+ for dom in self.domains:
+ self.io._initialize_octree(dom, self.oct_handler)
+ self.oct_handler.finalize()
+ tot = self.oct_handler.linearly_count()
+ mylog.info("Identified %0.3e octs", tot)
+
+ def _detect_fields(self):
+ # TODO: Add additional fields
+ pfl = []
+ for dom in self.domains:
+ fl = self.io._identify_fields(dom.domain_filename)
+ for f in fl:
+ if f not in pfl: pfl.append(f)
+ self.field_list = pfl
+ pf = self.parameter_file
+ pf.particle_types = tuple(set(pt for pt, pf in pfl))
+
+ def _setup_classes(self):
+ dd = self._get_data_reader_dict()
+ super(ParticleGeometryHandler, self)._setup_classes(dd)
+ self.object_types.sort()
+
+ def _identify_base_chunk(self, dobj):
+ if getattr(dobj, "_chunk_info", None) is None:
+ mask = dobj.selector.select_octs(self.oct_handler)
+ counts = self.oct_handler.count_cells(dobj.selector, mask)
+ subsets = [ParticleDomainSubset(d, mask, c)
+ for d, c in zip(self.domains, counts) if c > 0]
+ dobj._chunk_info = subsets
+ dobj.size = sum(counts)
+ dobj.shape = (dobj.size,)
+ dobj._current_chunk = list(self._chunk_all(dobj))[0]
+
+ def _chunk_all(self, dobj):
+ oobjs = getattr(dobj._current_chunk, "objs", dobj._chunk_info)
+ yield YTDataChunk(dobj, "all", oobjs, dobj.size)
+
+ def _chunk_spatial(self, dobj, ngz):
+ raise NotImplementedError
+
+ def _chunk_io(self, dobj):
+ oobjs = getattr(dobj._current_chunk, "objs", dobj._chunk_info)
+ for subset in oobjs:
+ yield YTDataChunk(dobj, "io", [subset], subset.cell_count)
class OWLSStaticOutput(StaticOutput):
+ _hierarchy_class = ParticleGeometryHandler
_fieldinfo_fallback = OWLSFieldInfo
_fieldinfo_known = KnownOWLSFields
def __init__(self, filename, data_style="OWLS", root_dimensions = 64):
self._root_dimensions = root_dimensions
+ # Set up the template for domain files
+ self.storage_filename = None
super(OWLSStaticOutput, self).__init__(filename, data_style)
def __repr__(self):
return os.path.basename(self.parameter_filename).split(".")[0]
def _set_units(self):
- pass
+ self.units = {}
+ self.time_units = {}
+ self.conversion_factors = {}
def _parse_parameter_file(self):
handle = h5py.File(self.parameter_filename)
@@ -91,6 +201,11 @@
self.hubble_constant = hvals["HubbleParam"]
self.parameters = hvals
+ prefix = self.parameter_filename.split(".", 1)[0]
+ suffix = self.parameter_filename.rsplit(".", 1)[-1]
+ self.domain_template = "%s.%%i.%s" % (prefix, suffix)
+ self.domain_count = hvals["NumFilesPerSnapshot"]
+
handle.close()
@classmethod
diff -r a91f4af645706634a602ef6ef36768e0ea48ce71 -r 603c6241559a80611ce716c3cb8d588d57175582 yt/frontends/sph/io.py
--- a/yt/frontends/sph/io.py
+++ b/yt/frontends/sph/io.py
@@ -26,6 +26,7 @@
import h5py
import numpy as np
from yt.funcs import *
+from yt.utilities.exceptions import *
from yt.utilities.io_handler import \
BaseIOHandler
@@ -36,39 +37,74 @@
def _read_fluid_selection(self, chunks, selector, fields, size):
raise NotImplementedError
- def _read_particle_selection_by_type(self, chunks, selector, fields):
+ def _read_particle_selection(self, chunks, selector, fields):
rv = {}
# We first need a set of masks for each particle type
- ptf = defaultdict(lambda: list)
+ ptf = defaultdict(list)
psize = defaultdict(lambda: 0)
+ chunks = list(chunks)
for ftype, fname in fields:
ptf[ftype].append(fname)
for chunk in chunks: # Will be OWLS domains
for subset in chunk.objs:
for ptype, field_list in sorted(ptf.items()):
- f = h5py.File(subset.domain.filename, "r")
- coords = f["/PartType%s/Coordinates" % ptype][:]
- psize[ptype] += subset.count_particles(selector,
- coords[:,0], coords[:,1], coords[:,2])
+ f = h5py.File(subset.domain.domain_filename, "r")
+ coords = f["/%s/Coordinates" % ptype][:].astype("float64")
+ psize[ptype] += selector.count_points(
+ coords[:,0], coords[:,1], coords[:,2])
del coords
+ f.close()
# Now we have all the sizes, and we can allocate
ind = {}
for field in fields:
- rv[field] = np.empty(size, dtype="float64")
+ mylog.debug("Allocating %s values for %s", psize[field[0]], field)
+ rv[field] = np.empty(psize[field[0]], dtype="float64")
ind[field] = 0
for chunk in chunks: # Will be OWLS domains
for subset in chunk.objs:
for ptype, field_list in sorted(ptf.items()):
- f = h5py.File(subset.domain.filename, "r")
- g = f["/PartType%s" % ptype]
- coords = g["Coordinates"][:]
- mask = subset.select_particles(selector,
+ f = h5py.File(subset.domain.domain_filename, "r")
+ g = f["/%s" % ptype]
+ coords = g["Coordinates"][:].astype("float64")
+ mask = selector.select_points(
coords[:,0], coords[:,1], coords[:,2])
del coords
if mask is None: continue
for field in field_list:
data = g[field][mask,...]
my_ind = ind[ptype, field]
+ mylog.debug("Filling from %s to %s with %s",
+ my_ind, my_ind+data.size, field)
rv[ptype, field][my_ind:my_ind + data.size] = data
ind[ptype, field] += data.shape[0]
+ f.close()
return rv
+
+ def _initialize_octree(self, domain, octree):
+ f = h5py.File(domain.domain_filename, "r")
+ for key in f.keys():
+ if not key.startswith("PartType"): continue
+ pos = f[key]["Coordinates"][:].astype("float64")
+ octree.add(pos, domain.domain_number)
+ f.close()
+
+ def _count_particles(self, domain_filename):
+ f = h5py.File(domain_filename, "r")
+ npart = f["/Header"].attrs["NumPart_ThisFile"].sum()
+ f.close()
+ return npart
+
+ def _identify_fields(self, domain_filename):
+ f = h5py.File(domain_filename, "r")
+ fields = []
+ for key in f.keys():
+ if not key.startswith("PartType"): continue
+ g = f[key]
+ #ptype = int(key[8:])
+ ptype = str(key)
+ for k in g.keys():
+ if not hasattr(g[k], "shape"): continue
+ # str => not unicode!
+ fields.append((ptype, str(k)))
+ f.close()
+ return fields
diff -r a91f4af645706634a602ef6ef36768e0ea48ce71 -r 603c6241559a80611ce716c3cb8d588d57175582 yt/geometry/geometry_handler.py
--- a/yt/geometry/geometry_handler.py
+++ b/yt/geometry/geometry_handler.py
@@ -59,12 +59,12 @@
mylog.debug("Setting up classes.")
self._setup_classes()
+ mylog.debug("Initializing data grid data IO")
+ self._setup_data_io()
+
mylog.debug("Setting up domain geometry.")
self._setup_geometry()
- mylog.debug("Initializing data grid data IO")
- self._setup_data_io()
-
mylog.debug("Detecting fields.")
self._detect_fields()
diff -r a91f4af645706634a602ef6ef36768e0ea48ce71 -r 603c6241559a80611ce716c3cb8d588d57175582 yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -883,3 +883,25 @@
tnp += neighbors[i].sd.np
return tnp
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def count_cells(self, SelectorObject selector,
+ np.ndarray[np.uint8_t, ndim=2, cast=True] mask):
+ cdef int i, j, k, oi
+ # pos here is CELL center, not OCT center.
+ cdef np.float64_t pos[3]
+ cdef int n = mask.shape[0]
+ cdef np.float64_t base_dx[3], dx[3]
+ cdef int eterm[3]
+ cdef np.ndarray[np.int64_t, ndim=1] count
+ count = np.zeros(self.max_domain + 1, 'int64')
+ print count.shape[0], mask.shape[0], mask.shape[1], self.nocts
+ for oi in range(n):
+ o = self.oct_list[oi]
+ if o.domain == -1: continue
+ for i in range(8):
+ count[o.domain] += mask[oi,i]
+ return count
+
+
diff -r a91f4af645706634a602ef6ef36768e0ea48ce71 -r 603c6241559a80611ce716c3cb8d588d57175582 yt/utilities/exceptions.py
--- a/yt/utilities/exceptions.py
+++ b/yt/utilities/exceptions.py
@@ -208,3 +208,10 @@
def __str__(self):
return "Enzo test output file (OutputLog) not generated for: " + \
"'%s'" % (self.testname) + ".\nTest did not complete."
+
+class YTFieldNotParseable(YTException):
+ def __init__(self, field):
+ self.field = field
+
+ def __str__(self):
+ return "Cannot identify field %s" % self.field
https://bitbucket.org/yt_analysis/yt-3.0/commits/aae18b668296/
changeset: aae18b668296
branch: yt-3.0
user: MatthewTurk
date: 2012-12-21 22:25:35
summary: Merging.
affected #: 2 files
diff -r 603c6241559a80611ce716c3cb8d588d57175582 -r aae18b668296b5abf5e795c196bad3466f187b33 yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -24,6 +24,7 @@
"""
from libc.stdlib cimport malloc, free, qsort
+from libc.math cimport floor
cimport numpy as np
import numpy as np
from oct_container cimport Oct, OctAllocationContainer, OctreeContainer
@@ -145,7 +146,7 @@
for i in range(3):
pp[i] = ppos[i] - self.DLE[i]
dds[i] = (self.DRE[i] - self.DLE[i])/self.nn[i]
- ind[i] = <np.int64_t> (pp[i]/dds[i])
+ ind[i] = <np.int64_t> (floor(pp[i]/dds[i]))
cp[i] = (ind[i] + 0.5) * dds[i]
cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
while cur.children[0][0][0] != NULL:
@@ -311,6 +312,25 @@
count[cur.my_octs[i - cur.offset].domain - 1] += 1
return count
+ def check(self, int curdom):
+ cdef int dind, pi
+ cdef Oct oct
+ cdef OctAllocationContainer *cont = self.domains[curdom - 1]
+ cdef int nbad = 0
+ for pi in range(cont.n_assigned):
+ oct = cont.my_octs[pi]
+ for i in range(2):
+ for j in range(2):
+ for k in range(2):
+ if oct.children[i][j][k] != NULL and \
+ oct.children[i][j][k].level != oct.level + 1:
+ if curdom == 61:
+ print pi, oct.children[i][j][k].level,
+ print oct.level
+ nbad += 1
+ print "DOMAIN % 3i HAS % 9i BAD OCTS (%s / %s / %s)" % (curdom, nbad,
+ cont.n - cont.n_assigned, cont.n_assigned, cont.n)
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -691,7 +711,6 @@
malloc(sizeof(ParticleArrays))
cdef int i, j, k
my_oct.ind = my_oct.domain = -1
- my_oct.domain = -1
my_oct.local_ind = self.nocts - 1
my_oct.pos[0] = my_oct.pos[1] = my_oct.pos[2] = -1
my_oct.level = -1
diff -r 603c6241559a80611ce716c3cb8d588d57175582 -r aae18b668296b5abf5e795c196bad3466f187b33 yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -758,6 +758,62 @@
if self._plot_type is None:
self._plot_type = kwargs.pop("plot_type")
PWViewer.__init__(self, *args, **kwargs)
+
+ def _setup_origin(self):
+ origin = self.origin
+ axis_index = self.data_source.axis
+ if isinstance(origin, basestring):
+ origin = tuple(origin.split('-'))[:3]
+ if 1 == len(origin):
+ origin = ('lower', 'left') + origin
+ elif 2 == len(origin) and origin[0] in set(['left','right','center']):
+ o0map = {'left': 'lower', 'right': 'upper', 'center': 'center'}
+ origin = (o0map[origin[0]],) + origin
+ elif 2 == len(origin) and origin[0] in set(['lower','upper','center']):
+ origin = (origin[0], 'center', origin[-1])
+ assert origin[-1] in ['window', 'domain']
+
+ if origin[2] == 'window':
+ xllim, xrlim = self.xlim
+ yllim, yrlim = self.ylim
+ elif origin[2] == 'domain':
+ xllim = self.pf.domain_left_edge[x_dict[axis_index]]
+ xrlim = self.pf.domain_right_edge[x_dict[axis_index]]
+ yllim = self.pf.domain_left_edge[y_dict[axis_index]]
+ yrlim = self.pf.domain_right_edge[y_dict[axis_index]]
+ else:
+ mylog.warn("origin = {0}".format(origin))
+ msg = ('origin keyword "{0}" not recognized, must declare "domain" '
+ 'or "center" as the last term in origin.').format(self.origin)
+ raise RuntimeError(msg)
+
+ if origin[0] == 'lower':
+ yc = yllim
+ elif origin[0] == 'upper':
+ yc = yrlim
+ elif origin[0] == 'center':
+ yc = (yllim + yrlim)/2.0
+ else:
+ mylog.warn("origin = {0}".format(origin))
+ msg = ('origin keyword "{0}" not recognized, must declare "lower" '
+ '"upper" or "center" as the first term in origin.')
+ msg = msg.format(self.origin)
+ raise RuntimeError(msg)
+
+ if origin[1] == 'left':
+ xc = xllim
+ elif origin[1] == 'right':
+ xc = xrlim
+ elif origin[1] == 'center':
+ xc = (xllim + xrlim)/2.0
+ else:
+ mylog.warn("origin = {0}".format(origin))
+ msg = ('origin keyword "{0}" not recognized, must declare "left" '
+ '"right" or "center" as the second term in origin.')
+ msg = msg.format(self.origin)
+ raise RuntimeError(msg)
+
+ return xc, yc
def _setup_plots(self):
if self._current_field is not None:
@@ -769,20 +825,7 @@
md = self.get_metadata(f, strip_mathml = False, return_string = False)
axis_index = self.data_source.axis
- if self.origin == 'center-window':
- xc = (self.xlim[0]+self.xlim[1])/2
- yc = (self.ylim[0]+self.ylim[1])/2
- elif self.origin == 'center-domain':
- xc = (self.pf.domain_left_edge[x_dict[axis_index]]+
- self.pf.domain_right_edge[x_dict[axis_index]])/2
- yc = (self.pf.domain_left_edge[y_dict[axis_index]]+
- self.pf.domain_right_edge[y_dict[axis_index]])/2
- elif self.origin == 'left-domain':
- xc = self.pf.domain_left_edge[x_dict[axis_index]]
- yc = self.pf.domain_left_edge[y_dict[axis_index]]
- else:
- raise RuntimeError(
- 'origin keyword: \"%(k)s\" not recognized' % {'k': self.origin})
+ xc, yc = self._setup_origin()
(unit_x, unit_y) = md['axes_unit_names']
@@ -1045,12 +1088,32 @@
Defaults to None, which automatically picks an appropriate unit.
If axes_unit is '1', 'u', or 'unitary', it will not display the
units, and only show the axes name.
- origin : string
- The location of the origin of the plot coordinate system.
- Currently, can be set to three options: 'left-domain', corresponding
- to the bottom-left hand corner of the simulation domain, 'center-domain',
- corresponding the center of the simulation domain, or 'center-window' for
- the center of the plot window.
+ origin : string or length 1, 2, or 3 sequence of strings
+ The location of the origin of the plot coordinate system. This is
+ represented by '-' separated string or a tuple of strings. In the
+ first index the y-location is given by 'lower', 'upper', or 'center'.
+ The second index is the x-location, given as 'left', 'right', or
+ 'center'. Finally, the whether the origin is applied in 'domain' space
+ or plot 'window' space is given. For example, both 'upper-right-domain'
+ and ['upper', 'right', 'domain'] both place the origin in the upper
+ right hand corner of domain space. If x or y are not given, a value is
+ inffered. For instance, 'left-domain' corresponds to the lower-left
+ hand corner of the simulation domain, 'center-domain' corresponds to the
+ center of the simulation domain, or 'center-window' for the center of
+ the plot window. Further examples:
+
+ ================================== ============================
+ format example
+ ================================== ============================
+ '{space}' 'domain'
+ '{xloc}-{space}' 'left-window'
+ '{yloc}-{space}' 'upper-domain'
+ '{yloc}-{xloc}-{space}' 'lower-right-window'
+ ('{space}',) ('window',)
+ ('{xloc}', '{space}') ('right', 'domain')
+ ('{yloc}', '{space}') ('lower', 'window')
+ ('{yloc}', '{xloc}', '{space}') ('lower', 'right', 'window')
+ ================================== ============================
fontsize : integer
The size of the fonts for the axis, colorbar, and tick labels.
@@ -1133,12 +1196,33 @@
Defaults to None, which automatically picks an appropriate unit.
If axes_unit is '1', 'u', or 'unitary', it will not display the
units, and only show the axes name.
- origin : A string
- The location of the origin of the plot coordinate system.
- Currently, can be set to three options: 'left-domain', corresponding
- to the bottom-left hand corner of the simulation domain, 'center-domain',
- corresponding the center of the simulation domain, or 'center-window' for
- the center of the plot window.
+ origin : string or length 1, 2, or 3 sequence of strings
+ The location of the origin of the plot coordinate system. This is
+ represented by '-' separated string or a tuple of strings. In the
+ first index the y-location is given by 'lower', 'upper', or 'center'.
+ The second index is the x-location, given as 'left', 'right', or
+ 'center'. Finally, the whether the origin is applied in 'domain' space
+ or plot 'window' space is given. For example, both 'upper-right-domain'
+ and ['upper', 'right', 'domain'] both place the origin in the upper
+ right hand corner of domain space. If x or y are not given, a value is
+ inffered. For instance, 'left-domain' corresponds to the lower-left
+ hand corner of the simulation domain, 'center-domain' corresponds to the
+ center of the simulation domain, or 'center-window' for the center of
+ the plot window.Further examples:
+
+ ================================== ============================
+ format example
+ ================================== ============================
+ '{space}' 'domain'
+ '{xloc}-{space}' 'left-window'
+ '{yloc}-{space}' 'upper-domain'
+ '{yloc}-{xloc}-{space}' 'lower-right-window'
+ ('{space}',) ('window',)
+ ('{xloc}', '{space}') ('right', 'domain')
+ ('{yloc}', '{space}') ('lower', 'window')
+ ('{yloc}', '{xloc}', '{space}') ('lower', 'right', 'window')
+ ================================== ============================
+
weight_field : string
The name of the weighting field. Set to None for no weight.
fontsize : integer
https://bitbucket.org/yt_analysis/yt-3.0/commits/b4b5540f4967/
changeset: b4b5540f4967
branch: yt-3.0
user: MatthewTurk
date: 2012-12-22 21:16:28
summary: Ensure vector fields are passed through correctly in SPH data.
affected #: 1 file
diff -r aae18b668296b5abf5e795c196bad3466f187b33 -r b4b5540f4967d5922a41a82bd64072629cd49c66 yt/frontends/sph/io.py
--- a/yt/frontends/sph/io.py
+++ b/yt/frontends/sph/io.py
@@ -31,6 +31,8 @@
from yt.utilities.io_handler import \
BaseIOHandler
+_vector_fields = ("Coordinates", "Velocity")
+
class IOHandlerOWLS(BaseIOHandler):
_data_style = "OWLS"
@@ -58,7 +60,11 @@
ind = {}
for field in fields:
mylog.debug("Allocating %s values for %s", psize[field[0]], field)
- rv[field] = np.empty(psize[field[0]], dtype="float64")
+ if field[1] in _vector_fields:
+ shape = (psize[field[0]], 3)
+ else:
+ shape = psize[field[0]]
+ rv[field] = np.empty(shape, dtype="float64")
ind[field] = 0
for chunk in chunks: # Will be OWLS domains
for subset in chunk.objs:
@@ -74,8 +80,8 @@
data = g[field][mask,...]
my_ind = ind[ptype, field]
mylog.debug("Filling from %s to %s with %s",
- my_ind, my_ind+data.size, field)
- rv[ptype, field][my_ind:my_ind + data.size] = data
+ my_ind, my_ind+data.shape[0], field)
+ rv[ptype, field][my_ind:my_ind + data.shape[0],...] = data
ind[ptype, field] += data.shape[0]
f.close()
return rv
https://bitbucket.org/yt_analysis/yt-3.0/commits/0bba17a0eed3/
changeset: 0bba17a0eed3
branch: yt-3.0
user: MatthewTurk
date: 2013-01-04 23:08:22
summary: Updated to version 0.4.5 of voro++.
affected #: 1 file
diff -r b4b5540f4967d5922a41a82bd64072629cd49c66 -r 0bba17a0eed3a331925b6849a1f30650e3912c63 yt/utilities/voropp.pyx
--- a/yt/utilities/voropp.pyx
+++ b/yt/utilities/voropp.pyx
@@ -31,19 +31,33 @@
cimport numpy as np
cimport cython
-cdef extern from "voro++.cc":
+cdef extern from "voro++.hh" namespace "voro":
+ cdef cppclass c_loop_all
+
+ cdef cppclass voronoicell:
+ double volume()
+
cdef cppclass container:
container(double xmin, double xmax, double ymin, double ymax,
double zmin, double zmax, int nx, int ny, int nz,
libcpp.bool xper, libcpp.bool yper, libcpp.bool zper, int alloc)
void put(int n, double x, double y, double z)
void store_cell_volumes(double *vols)
+ int compute_cell(voronoicell c, c_loop_all vl)
+ double sum_cell_volumes()
+
+ cdef cppclass c_loop_all:
+ c_loop_all(container &con)
+ int inc()
+ int start()
cdef class VoronoiVolume:
cdef container *my_con
- cdef int npart
- def __init__(self, xi, yi, zi):
- self.my_con = new container(0.0, 1.0, 0.0, 1.0, 0.0, 1.0,
+ cdef public int npart
+ def __init__(self, xi, yi, zi, left_edge, right_edge):
+ self.my_con = new container(left_edge[0], right_edge[0],
+ left_edge[1], right_edge[1],
+ left_edge[2], right_edge[2],
xi, yi, zi, False, False, False, 8)
self.npart = 0
@@ -65,5 +79,15 @@
def get_volumes(self):
cdef np.ndarray vol = np.zeros(self.npart, 'double')
cdef double *vdouble = <double *> vol.data
- self.my_con.store_cell_volumes(vdouble)
+ #self.my_con.store_cell_volumes(vdouble)
+ cdef c_loop_all *vl = new c_loop_all(deref(self.my_con))
+ cdef voronoicell c
+ if not vl.start(): return
+ cdef int i = 0
+ while 1:
+ if self.my_con.compute_cell(c, deref(vl)):
+ vol[i] = c.volume()
+ if not vl.inc(): break
+ i += 1
+ del vl
return vol
https://bitbucket.org/yt_analysis/yt-3.0/commits/db23b0022ebc/
changeset: db23b0022ebc
branch: yt-3.0
user: MatthewTurk
date: 2013-01-07 19:34:53
summary: This fixes a terribly pernicious indexing bug in the particle octree code and
adds tests for the particle octree.
affected #: 2 files
diff -r 0bba17a0eed3a331925b6849a1f30650e3912c63 -r db23b0022ebcadede5a8ae3733bb1dcc13513c2e yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -146,7 +146,7 @@
for i in range(3):
pp[i] = ppos[i] - self.DLE[i]
dds[i] = (self.DRE[i] - self.DLE[i])/self.nn[i]
- ind[i] = <np.int64_t> (floor(pp[i]/dds[i]))
+ ind[i] = <np.int64_t> ((pp[i] - self.DLE[i])/dds[i])
cp[i] = (ind[i] + 0.5) * dds[i]
cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
while cur.children[0][0][0] != NULL:
@@ -695,14 +695,16 @@
i += 1
qsort(self.oct_list, self.nocts, sizeof(Oct*), &compare_octs)
cdef int cur_dom = -1
+ # We always need at least 2, and if max_domain is 0, we need 3.
self.dom_offsets = <np.int64_t *>malloc(sizeof(np.int64_t) *
- self.max_domain + 2)
+ (self.max_domain + 3))
+ self.dom_offsets[0] = 0
for i in range(self.nocts):
self.oct_list[i].local_ind = i
if self.oct_list[i].domain > cur_dom:
cur_dom = self.oct_list[i].domain
- self.dom_offsets[cur_dom] = i
- self.dom_offsets[cur_dom + 1] = self.nocts
+ self.dom_offsets[cur_dom + 1] = i
+ self.dom_offsets[cur_dom + 2] = self.nocts
cdef Oct* allocate_oct(self):
self.nocts += 1
@@ -752,12 +754,10 @@
cdef Oct *o
cdef int oi, i
level_count = np.zeros(max_level+1, 'int64')
- ndo
cdef np.int64_t ndo, doff
- ndo = self.dom_offsets[domain_id + 1] \
- - self.dom_offsets[domain_id]
- doff = self.dom_offsets[domain_id]
- print "Domain %s Total %s" % (domain_id, ndo)
+ ndo = self.dom_offsets[domain_id + 2] \
+ - self.dom_offsets[domain_id + 1]
+ doff = self.dom_offsets[domain_id + 1]
for oi in range(ndo):
o = self.oct_list[oi + doff]
for i in range(8):
@@ -778,7 +778,7 @@
for i in range(3):
pp[i] = pos[p, i]
dds[i] = (self.DRE[i] + self.DLE[i])/self.nn[i]
- ind[i] = <np.int64_t> (pp[i]/dds[i])
+ ind[i] = <np.int64_t> ((pp[i] - self.DLE[i])/dds[i])
cp[i] = (ind[i] + 0.5) * dds[i]
cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
if cur == NULL:
diff -r 0bba17a0eed3a331925b6849a1f30650e3912c63 -r db23b0022ebcadede5a8ae3733bb1dcc13513c2e yt/geometry/tests/test_particle_octree.py
--- /dev/null
+++ b/yt/geometry/tests/test_particle_octree.py
@@ -0,0 +1,39 @@
+from yt.testing import *
+import numpy as np
+from yt.geometry.oct_container import ParticleOctreeContainer
+import time, os
+
+NPART = 32**3
+NDIM = 64
+DLE = np.array([0.0, 0.0, 0.0])
+DRE = np.array([10.0, 10.0, 10.0])
+
+def test_add_particles_random():
+ np.random.seed(int(0x4d3d3d3))
+ pos = np.random.normal(0.5, scale=0.05, size=(NPART,3)) * (DRE-DLE) + DLE
+ for i in range(3):
+ np.clip(pos[:,i], DLE[i], DRE[i], pos[:,i])
+ for ndom in [1, 2, 4, 8]:
+ octree = ParticleOctreeContainer((NDIM, NDIM, NDIM), DLE, DRE)
+ for dom in range(ndom):
+ octree.add(pos[dom::ndom,:], dom)
+ octree.finalize()
+ # This visits every oct.
+ lin_count = octree.linearly_count()
+ tc = octree.recursively_count()
+ total_count = np.zeros(len(tc), dtype="int32")
+ for i in sorted(tc):
+ total_count[i] = tc[i]
+ yield assert_equal, lin_count, total_count.sum()
+ mask = np.ones((total_count.sum(), 8), dtype="bool")
+ # This visits every cell -- including those covered by octs.
+ level_count = octree.count_levels(total_count.size-1, -1, mask)
+ for dom in range(ndom):
+ level_count += octree.count_levels(total_count.size-1, dom, mask)
+ yield assert_equal, level_count[0], NDIM**3 * 8
+ yield assert_equal, level_count, total_count * 8
+
+if __name__=="__main__":
+ for i in test_add_particles_random():
+ i[0](*i[1:])
+ time.sleep(1)
https://bitbucket.org/yt_analysis/yt-3.0/commits/d000eb4f9e52/
changeset: d000eb4f9e52
branch: yt-3.0
user: ngoldbaum
date: 2013-01-10 15:56:32
summary: Merged in MatthewTurk/yt-3.0 (pull request #11: Initial N-body data support)
affected #: 34 files
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/enzo/data_structures.py
--- a/yt/frontends/enzo/data_structures.py
+++ b/yt/frontends/enzo/data_structures.py
@@ -463,7 +463,6 @@
else:
field_list = None
field_list = self.comm.mpi_bcast(field_list)
- self.save_data(list(field_list),"/","DataFields",passthrough=True)
self.field_list = list(field_list)
def _generate_random_grids(self):
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/gadget/__init__.py
--- a/yt/frontends/gadget/__init__.py
+++ /dev/null
@@ -1,29 +0,0 @@
-"""
-API for yt.frontends.gadget
-
-Author: Matthew Turk <matthewturk at gmail.com>
-Affiliation: UCSD
-Author: J.S. Oishi <jsoishi at gmail.com>
-Affiliation: KIPAC/SLAC/Stanford
-Author: Britton Smith <brittonsmith at gmail.com>
-Affiliation: MSU
-Homepage: http://yt-project.org/
-License:
- Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
-
- This file is part of yt.
-
- yt is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-"""
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/gadget/api.py
--- a/yt/frontends/gadget/api.py
+++ /dev/null
@@ -1,41 +0,0 @@
-"""
-API for yt.frontends.gadget
-
-Author: Matthew Turk <matthewturk at gmail.com>
-Affiliation: UCSD
-Author: J.S. Oishi <jsoishi at gmail.com>
-Affiliation: KIPAC/SLAC/Stanford
-Author: Britton Smith <brittonsmith at gmail.com>
-Affiliation: MSU
-Homepage: http://yt-project.org/
-License:
- Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
-
- This file is part of yt.
-
- yt is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-"""
-
-from .data_structures import \
- GadgetGrid, \
- GadgetHierarchy, \
- GadgetStaticOutput
-
-from .fields import \
- GadgetFieldInfo, \
- add_gadget_field
-
-from .io import \
- IOHandlerGadget
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/gadget/data_structures.py
--- a/yt/frontends/gadget/data_structures.py
+++ /dev/null
@@ -1,201 +0,0 @@
-"""
-Data structures for Gadget.
-
-Author: Matthew Turk <matthewturk at gmail.com>
-Affiliation: UCSD
-Author: Chris Moody <cemoody at ucsc.edu>
-Affiliation: UCSC
-Homepage: http://yt-project.org/
-License:
- Copyright (C) 2007-2011 Matthew Turk. All Rights Reserved.
-
- This file is part of yt.
-
- yt is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-"""
-
-import h5py
-import numpy as np
-from itertools import izip
-
-from yt.funcs import *
-from yt.data_objects.grid_patch import \
- AMRGridPatch
-from yt.geometry.grid_geometry_handler import \
- GridGeometryHandler
-from yt.data_objects.static_output import \
- StaticOutput
-from yt.utilities.definitions import \
- sec_conversion
-
-from .fields import GadgetFieldInfo, KnownGadgetFields
-from yt.data_objects.field_info_container import \
- FieldInfoContainer, NullFunc
-
-class GadgetGrid(AMRGridPatch):
- _id_offset = 0
- def __init__(self, hierarchy, id, dimensions, start,
- level, parent_id, particle_count):
- AMRGridPatch.__init__(self, id, filename = hierarchy.filename,
- hierarchy = hierarchy)
- self.Parent = [] # Only one parent per grid
- self.Children = []
- self.Level = level
- self.ActiveDimensions = dimensions.copy()
- self.NumberOfParticles = particle_count
- self.start_index = start.copy().astype("int64")
- self.stop_index = self.start_index + dimensions.copy()
- self.id = id
- self._parent_id = parent_id
-
- try:
- padd = '/data/grid_%010i/particles' % id
- self.particle_types = self.hierarchy._handle[padd].keys()
- except:
- self.particle_types = ()
- self.filename = hierarchy.filename
-
- def __repr__(self):
- return 'GadgetGrid_%05i'%self.id
-
-class GadgetHierarchy(GridGeometryHandler):
- grid = GadgetGrid
-
- def __init__(self, pf, data_style='gadget_hdf5'):
- self.filename = pf.filename
- self.directory = os.path.dirname(pf.filename)
- self.data_style = data_style
- self._handle = h5py.File(pf.filename)
- GridGeometryHandler.__init__(self, pf, data_style)
- self._handle.close()
- self._handle = None
-
-
- def _initialize_data_storage(self):
- pass
-
- def _detect_fields(self):
- #this adds all the fields in
- #/particle_types/{Gas/Stars/etc.}/{position_x, etc.}
- self.field_list = []
- for ptype in self._handle['particle_types'].keys():
- for field in self._handle['particle_types'+'/'+ptype]:
- if field not in self.field_list:
- self.field_list += field,
-
- def _setup_classes(self):
- dd = self._get_data_reader_dict()
- GridGeometryHandler._setup_classes(self, dd)
- self.object_types.sort()
-
- def _count_grids(self):
- self.num_grids = len(self._handle['/grid_dimensions'])
-
- def _parse_hierarchy(self):
- f = self._handle # shortcut
- npa = np.array
- DLE = self.parameter_file.domain_left_edge
- DRE = self.parameter_file.domain_right_edge
- DW = (DRE - DLE)
-
- self.grid_levels.flat[:] = f['/grid_level'][:].astype("int32")
- LI = f['/grid_left_index'][:]
- print LI
- self.grid_dimensions[:] = f['/grid_dimensions'][:]
- self.grid_left_edge[:] = (LI * DW + DLE)
- dxs = 1.0/(2**(self.grid_levels+1)) * DW
- self.grid_right_edge[:] = self.grid_left_edge \
- + dxs *(1 + self.grid_dimensions)
- self.grid_particle_count.flat[:] = f['/grid_particle_count'][:].astype("int32")
- grid_parent_id = f['/grid_parent_id'][:]
- self.max_level = np.max(self.grid_levels)
-
- args = izip(xrange(self.num_grids), self.grid_levels.flat,
- grid_parent_id, LI,
- self.grid_dimensions, self.grid_particle_count.flat)
- self.grids = np.empty(len(args), dtype='object')
- for gi, (j,lvl,p, le, d, n) in enumerate(args):
- self.grids[gi] = self.grid(self,j,d,le,lvl,p,n)
-
- def _populate_grid_objects(self):
- for g in self.grids:
- if g._parent_id >= 0:
- parent = self.grids[g._parent_id]
- g.Parent = parent
- parent.Children.append(g)
- for g in self.grids:
- g._prepare_grid()
- g._setup_dx()
-
- def _setup_derived_fields(self):
- self.derived_field_list = []
-
-class GadgetStaticOutput(StaticOutput):
- _hierarchy_class = GadgetHierarchy
- _fieldinfo_fallback = GadgetFieldInfo
- _fieldinfo_known = KnownGadgetFields
-
- def __init__(self, filename,storage_filename=None) :
- self.storage_filename = storage_filename
- self.filename = filename
-
- StaticOutput.__init__(self, filename, 'gadget_infrastructure')
- self._set_units()
-
- def _set_units(self):
- self.units = {}
- self.time_units = {}
- self.time_units['1'] = 1
- self.units['1'] = 1.0
- self.units['cm'] = 1.0
- self.units['unitary'] = 1.0 / \
- (self.domain_right_edge - self.domain_left_edge).max()
- for unit in sec_conversion.keys():
- self.time_units[unit] = 1.0 / sec_conversion[unit]
-
- def _parse_parameter_file(self):
- fileh = h5py.File(self.filename)
- sim_param = fileh['/simulation_parameters'].attrs
- self.refine_by = sim_param['refine_by']
- self.dimensionality = sim_param['dimensionality']
- self.num_ghost_zones = sim_param['num_ghost_zones']
- self.field_ordering = sim_param['field_ordering']
- self.domain_dimensions = sim_param['domain_dimensions']
- self.current_time = sim_param['current_time']
- self.domain_left_edge = sim_param['domain_left_edge']
- self.domain_right_edge = sim_param['domain_right_edge']
- self.unique_identifier = sim_param['unique_identifier']
- self.cosmological_simulation = sim_param['cosmological_simulation']
- self.current_redshift = sim_param['current_redshift']
- self.omega_lambda = sim_param['omega_lambda']
- self.hubble_constant = sim_param['hubble_constant']
- fileh.close()
-
-
- @classmethod
- def _is_valid(self, *args, **kwargs):
- format = 'Gadget Infrastructure'
- add1 = 'griddded_data_format'
- add2 = 'data_software'
- try:
- fileh = h5py.File(args[0],'r')
- if add1 in fileh['/'].items():
- if add2 in fileh['/'+add1].attrs.keys():
- if fileh['/'+add1].attrs[add2] == format:
- fileh.close()
- return True
- fileh.close()
- except:
- pass
- return False
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/gadget/fields.py
--- a/yt/frontends/gadget/fields.py
+++ /dev/null
@@ -1,133 +0,0 @@
-"""
-Gadget-specific fields
-
-Author: Christopher E Moody <juxtaposicion at gmail.com>
-Affiliation: UC Santa Cruz
-Homepage: http://yt-project.org/
-License:
- Copyright (C) 2010-2011 Christopher E Moody, Matthew Turk. All Rights Reserved.
-
- This file is part of yt.
-
- yt is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-"""
-
-import numpy as np
-
-from yt.funcs import *
-from yt.data_objects.field_info_container import \
- FieldInfoContainer, \
- FieldInfo, \
- ValidateParameter, \
- ValidateDataField, \
- ValidateProperty, \
- ValidateSpatial, \
- ValidateGridType
-import yt.data_objects.universal_fields
-
-GadgetFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
-add_gadget_field = GadgetFieldInfo.add_field
-
-add_field = add_gadget_field
-
-translation_dict = {"particle_position_x" : "position_x",
- "particle_position_y" : "position_y",
- "particle_position_z" : "position_z",
- }
-
-def _generate_translation(mine, theirs):
- pfield = mine.startswith("particle")
- add_field(theirs, function=lambda a, b: b[mine], take_log=True,
- particle_type = pfield)
-
-for f,v in translation_dict.items():
- if v not in GadgetFieldInfo:
- # Note here that it's the yt field that we check for particle nature
- pfield = f.startswith("particle")
- add_field(v, function=lambda a,b: None, take_log=False,
- validators = [ValidateDataField(v)],
- particle_type = pfield)
- print "Setting up translator from %s to %s" % (v, f)
- _generate_translation(v, f)
-
-
-#for f,v in translation_dict.items():
-# add_field(f, function=lambda a,b: None, take_log=True,
-# validators = [ValidateDataField(v)],
-# units=r"\rm{cm}")
-# add_field(v, function=lambda a,b: None, take_log=True,
-# validators = [ValidateDataField(v)],
-# units=r"\rm{cm}")
-
-
-
-add_field("position_x", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("position_x")],
- particle_type = True,
- units=r"\rm{cm}")
-
-add_field("position_y", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("position_y")],
- particle_type = True,
- units=r"\rm{cm}")
-
-add_field("position_z", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("position_z")],
- particle_type = True,
- units=r"\rm{cm}")
-
-add_field("VEL", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("VEL")],
- units=r"")
-
-add_field("id", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("ID")],
- units=r"")
-
-add_field("mass", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("mass")],
- units=r"\rm{g}")
-def _particle_mass(field, data):
- return data["mass"]/just_one(data["CellVolume"])
-def _convert_particle_mass(data):
- return 1.0
-add_field("particle_mass", function=_particle_mass, take_log=True,
- convert_function=_convert_particle_mass,
- validators = [ValidateSpatial(0)],
- units=r"\mathrm{g}\/\mathrm{cm}^{-3}")
-
-add_field("U", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("U")],
- units=r"")
-
-add_field("NE", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("NE")],
- units=r"")
-
-add_field("POT", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("POT")],
- units=r"")
-
-add_field("ACCE", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("ACCE")],
- units=r"")
-
-add_field("ENDT", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("ENDT")],
- units=r"")
-
-add_field("TSTP", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("TSTP")],
- units=r"")
-
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/gadget/io.py
--- a/yt/frontends/gadget/io.py
+++ /dev/null
@@ -1,55 +0,0 @@
-"""
-Gadget-specific data-file handling function
-
-Author: Christopher E Moody <juxtaposicion at gmail.com>
-Affiliation: UC Santa Cruz
-Homepage: http://yt-project.org/
-License:
- Copyright (C) 2010-2011 Christopher E Moody, Matthew Turk. All Rights Reserved.
-
- This file is part of yt.
-
- yt is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-"""
-
-import h5py
-import numpy as np
-
-from yt.utilities.io_handler import \
- BaseIOHandler
-
-class IOHandlerGadget(BaseIOHandler):
- _data_style = 'gadget_infrastructure'
- def _read_data_set(self, grid, field):
- data = []
- fh = h5py.File(grid.filename,mode='r')
- for ptype in grid.particle_types:
- address = '/data/grid_%010i/particles/%s/%s' % (grid.id, ptype, field)
- data.append(fh[address][:])
- if len(data) > 0:
- data = np.concatenate(data)
- fh.close()
- return np.array(data)
- def _read_field_names(self,grid):
- adr = grid.Address
- fh = h5py.File(grid.filename,mode='r')
- rets = cPickle.loads(fh['/root'].attrs['fieldnames'])
- fh.close()
- return rets
-
- def _read_data_slice(self,grid, field, axis, coord):
- #how would we implement axis here?
- dat = self._read_data_set(grid,field)
- return dat[coord]
-
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/gadget/setup.py
--- a/yt/frontends/gadget/setup.py
+++ /dev/null
@@ -1,13 +0,0 @@
-#!/usr/bin/env python
-import setuptools
-import os
-import sys
-import os.path
-
-
-def configuration(parent_package='', top_path=None):
- from numpy.distutils.misc_util import Configuration
- config = Configuration('gadget', parent_package, top_path)
- config.make_config_py() # installs __config__.py
- #config.make_svn_version_py()
- return config
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/gdf/data_structures.py
--- a/yt/frontends/gdf/data_structures.py
+++ b/yt/frontends/gdf/data_structures.py
@@ -247,7 +247,9 @@
try:
fileh = h5py.File(args[0],'r')
if "gridded_data_format" in fileh:
+ fileh.close()
return True
+ fileh.close()
except:
pass
return False
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/ramses/data_structures.py
--- a/yt/frontends/ramses/data_structures.py
+++ b/yt/frontends/ramses/data_structures.py
@@ -29,8 +29,6 @@
import cStringIO
from yt.funcs import *
-from yt.data_objects.grid_patch import \
- AMRGridPatch
from yt.geometry.oct_geometry_handler import \
OctreeGeometryHandler
from yt.geometry.geometry_handler import \
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/__init__.py
--- /dev/null
+++ b/yt/frontends/sph/__init__.py
@@ -0,0 +1,29 @@
+"""
+API for yt.frontends.gadget
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Author: J.S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: MSU
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/api.py
--- /dev/null
+++ b/yt/frontends/sph/api.py
@@ -0,0 +1,35 @@
+"""
+API for yt.frontends.sph
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Author: J.S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: MSU
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+
+from .data_structures import \
+ OWLSStaticOutput
+
+from .io import \
+ IOHandlerOWLS
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/data_structures.py
--- /dev/null
+++ b/yt/frontends/sph/data_structures.py
@@ -0,0 +1,222 @@
+"""
+Data structures for a generic SPH/Gadget frontend.
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2007-2012 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import h5py
+import numpy as np
+import stat
+import weakref
+from itertools import izip
+
+from yt.funcs import *
+from yt.geometry.oct_geometry_handler import \
+ OctreeGeometryHandler
+from yt.geometry.oct_container import \
+ ParticleOctreeContainer
+from yt.geometry.geometry_handler import \
+ GeometryHandler, YTDataChunk
+from yt.data_objects.static_output import \
+ StaticOutput
+from yt.utilities.definitions import \
+ mpc_conversion, sec_conversion
+from .fields import \
+ OWLSFieldInfo, \
+ KnownOWLSFields
+
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, NullFunc
+
+class ParticleDomainFile(object):
+ def __init__(self, pf, io, domain_filename, domain_number):
+ self.pf = pf
+ self.io = weakref.proxy(io)
+ self.domain_filename = domain_filename
+ self.domain_number = domain_number
+ self.total_particles = self.io._count_particles(domain_filename)
+
+ def select(self, selector):
+ pass
+
+ def count(self, selector):
+ pass
+
+class ParticleDomainSubset(object):
+ def __init__(self, domain, mask, count):
+ self.domain = domain
+ self.mask = mask
+ self.cell_count = count
+ self.oct_handler = domain.pf.h.oct_handler
+
+ def icoords(self, dobj):
+ return self.oct_handler.icoords(self.domain.domain_id, self.mask,
+ self.cell_count,
+ self.level_counts.copy())
+
+ def fcoords(self, dobj):
+ return self.oct_handler.fcoords(self.domain.domain_id, self.mask,
+ self.cell_count,
+ self.level_counts.copy())
+
+ def fwidth(self, dobj):
+ # Recall domain_dimensions is the number of cells, not octs
+ base_dx = 1.0/self.domain.pf.domain_dimensions
+ widths = np.empty((self.cell_count, 3), dtype="float64")
+ dds = (2**self.ires(dobj))
+ for i in range(3):
+ widths[:,i] = base_dx[i] / dds
+ return widths
+
+ def ires(self, dobj):
+ return self.oct_handler.ires(self.domain.domain_id, self.mask,
+ self.cell_count,
+ self.level_counts.copy())
+
+
+class ParticleGeometryHandler(OctreeGeometryHandler):
+
+ def __init__(self, pf, data_style):
+ self.data_style = data_style
+ self.parameter_file = weakref.proxy(pf)
+ # for now, the hierarchy file is the parameter file!
+ self.hierarchy_filename = self.parameter_file.parameter_filename
+ self.directory = os.path.dirname(self.hierarchy_filename)
+ self.float_type = np.float64
+ super(ParticleGeometryHandler, self).__init__(pf, data_style)
+
+ def _initialize_oct_handler(self):
+ template = self.parameter_file.domain_template
+ ndoms = self.parameter_file.domain_count
+ self.domains = [ParticleDomainFile(self.parameter_file, self.io, template % i, i)
+ for i in range(ndoms)]
+ total_particles = sum(d.total_particles for d in self.domains)
+ self.oct_handler = ParticleOctreeContainer(
+ self.parameter_file.domain_dimensions,
+ self.parameter_file.domain_left_edge,
+ self.parameter_file.domain_right_edge)
+ mylog.info("Allocating for %0.3e particles", total_particles)
+ for dom in self.domains:
+ self.io._initialize_octree(dom, self.oct_handler)
+ self.oct_handler.finalize()
+ tot = self.oct_handler.linearly_count()
+ mylog.info("Identified %0.3e octs", tot)
+
+ def _detect_fields(self):
+ # TODO: Add additional fields
+ pfl = []
+ for dom in self.domains:
+ fl = self.io._identify_fields(dom.domain_filename)
+ for f in fl:
+ if f not in pfl: pfl.append(f)
+ self.field_list = pfl
+ pf = self.parameter_file
+ pf.particle_types = tuple(set(pt for pt, pf in pfl))
+
+ def _setup_classes(self):
+ dd = self._get_data_reader_dict()
+ super(ParticleGeometryHandler, self)._setup_classes(dd)
+ self.object_types.sort()
+
+ def _identify_base_chunk(self, dobj):
+ if getattr(dobj, "_chunk_info", None) is None:
+ mask = dobj.selector.select_octs(self.oct_handler)
+ counts = self.oct_handler.count_cells(dobj.selector, mask)
+ subsets = [ParticleDomainSubset(d, mask, c)
+ for d, c in zip(self.domains, counts) if c > 0]
+ dobj._chunk_info = subsets
+ dobj.size = sum(counts)
+ dobj.shape = (dobj.size,)
+ dobj._current_chunk = list(self._chunk_all(dobj))[0]
+
+ def _chunk_all(self, dobj):
+ oobjs = getattr(dobj._current_chunk, "objs", dobj._chunk_info)
+ yield YTDataChunk(dobj, "all", oobjs, dobj.size)
+
+ def _chunk_spatial(self, dobj, ngz):
+ raise NotImplementedError
+
+ def _chunk_io(self, dobj):
+ oobjs = getattr(dobj._current_chunk, "objs", dobj._chunk_info)
+ for subset in oobjs:
+ yield YTDataChunk(dobj, "io", [subset], subset.cell_count)
+
+class OWLSStaticOutput(StaticOutput):
+ _hierarchy_class = ParticleGeometryHandler
+ _fieldinfo_fallback = OWLSFieldInfo
+ _fieldinfo_known = KnownOWLSFields
+
+ def __init__(self, filename, data_style="OWLS", root_dimensions = 64):
+ self._root_dimensions = root_dimensions
+ # Set up the template for domain files
+ self.storage_filename = None
+ super(OWLSStaticOutput, self).__init__(filename, data_style)
+
+ def __repr__(self):
+ return os.path.basename(self.parameter_filename).split(".")[0]
+
+ def _set_units(self):
+ self.units = {}
+ self.time_units = {}
+ self.conversion_factors = {}
+
+ def _parse_parameter_file(self):
+ handle = h5py.File(self.parameter_filename)
+ hvals = {}
+ hvals.update(handle["/Header"].attrs)
+
+ self.dimensionality = 3
+ self.refine_by = 2
+ self.parameters["HydroMethod"] = "sph"
+ self.unique_identifier = \
+ int(os.stat(self.parameter_filename)[stat.ST_CTIME])
+ # Set standard values
+ self.current_time = hvals["Time_GYR"] * sec_conversion["Gyr"]
+ self.domain_left_edge = np.zeros(3, "float64")
+ self.domain_right_edge = np.ones(3, "float64") * hvals["BoxSize"]
+ self.domain_dimensions = np.ones(3, "int32") * self._root_dimensions
+ self.cosmological_simulation = 1
+ self.current_redshift = hvals["Redshift"]
+ self.omega_lambda = hvals["OmegaLambda"]
+ self.omega_matter = hvals["Omega0"]
+ self.hubble_constant = hvals["HubbleParam"]
+ self.parameters = hvals
+
+ prefix = self.parameter_filename.split(".", 1)[0]
+ suffix = self.parameter_filename.rsplit(".", 1)[-1]
+ self.domain_template = "%s.%%i.%s" % (prefix, suffix)
+ self.domain_count = hvals["NumFilesPerSnapshot"]
+
+ handle.close()
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ try:
+ fileh = h5py.File(args[0],'r')
+ if "Constants" in fileh["/"].keys() and \
+ "Header" in fileh["/"].keys():
+ fileh.close()
+ return True
+ fileh.close()
+ except:
+ pass
+ return False
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/fields.py
--- /dev/null
+++ b/yt/frontends/sph/fields.py
@@ -0,0 +1,44 @@
+"""
+OWLS-specific fields
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2012 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+
+from yt.funcs import *
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, \
+ FieldInfo, \
+ ValidateParameter, \
+ ValidateDataField, \
+ ValidateProperty, \
+ ValidateSpatial, \
+ ValidateGridType
+import yt.data_objects.universal_fields
+
+OWLSFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
+add_owls_field = OWLSFieldInfo.add_field
+
+KnownOWLSFields = FieldInfoContainer()
+add_OWLS_field = KnownOWLSFields.add_field
+
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/gadget/__init__.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/__init__.py
@@ -0,0 +1,29 @@
+"""
+API for yt.frontends.gadget
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Author: J.S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: MSU
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/gadget/api.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/api.py
@@ -0,0 +1,30 @@
+"""
+API for yt.frontends.gadget
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Author: J.S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: MSU
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/gadget/data_structures.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/data_structures.py
@@ -0,0 +1,58 @@
+"""
+Data structures for a generic SPH/Gadget frontend.
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2007-2012 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import h5py
+import numpy as np
+from itertools import izip
+
+from yt.funcs import *
+from yt.geometry.oct_geometry_handler import \
+ OctreeGeometryHandler
+from yt.geometry.geometry_handler import \
+ GeometryHandler, YTDataChunk
+from yt.data_objects.static_output import \
+ StaticOutput
+
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, NullFunc
+
+class GadgetDomainFile(object):
+ def select(self, selector):
+ pass
+
+ def count(self, selector):
+ pass
+
+class GadgetDomainSubset(object):
+ def __init__(self, domain, mask, cell_count):
+ self.mask = mask
+ self.domain = domain
+ self.oct_handler = domain.pf.h.oct_handler
+ self.cell_count = cell_count
+ level_counts = self.oct_handler.count_levels(
+ self.domain.pf.max_level, self.domain.domain_id, mask)
+ level_counts[1:] = level_counts[:-1]
+ level_counts[0] = 0
+ self.level_counts = np.add.accumulate(level_counts)
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/gadget/fields.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/fields.py
@@ -0,0 +1,133 @@
+"""
+Gadget-specific fields
+
+Author: Christopher E Moody <juxtaposicion at gmail.com>
+Affiliation: UC Santa Cruz
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Christopher E Moody, Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+
+from yt.funcs import *
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, \
+ FieldInfo, \
+ ValidateParameter, \
+ ValidateDataField, \
+ ValidateProperty, \
+ ValidateSpatial, \
+ ValidateGridType
+import yt.data_objects.universal_fields
+
+GadgetFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
+add_gadget_field = GadgetFieldInfo.add_field
+
+add_field = add_gadget_field
+
+translation_dict = {"particle_position_x" : "position_x",
+ "particle_position_y" : "position_y",
+ "particle_position_z" : "position_z",
+ }
+
+def _generate_translation(mine, theirs):
+ pfield = mine.startswith("particle")
+ add_field(theirs, function=lambda a, b: b[mine], take_log=True,
+ particle_type = pfield)
+
+for f,v in translation_dict.items():
+ if v not in GadgetFieldInfo:
+ # Note here that it's the yt field that we check for particle nature
+ pfield = f.startswith("particle")
+ add_field(v, function=lambda a,b: None, take_log=False,
+ validators = [ValidateDataField(v)],
+ particle_type = pfield)
+ print "Setting up translator from %s to %s" % (v, f)
+ _generate_translation(v, f)
+
+
+#for f,v in translation_dict.items():
+# add_field(f, function=lambda a,b: None, take_log=True,
+# validators = [ValidateDataField(v)],
+# units=r"\rm{cm}")
+# add_field(v, function=lambda a,b: None, take_log=True,
+# validators = [ValidateDataField(v)],
+# units=r"\rm{cm}")
+
+
+
+add_field("position_x", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("position_x")],
+ particle_type = True,
+ units=r"\rm{cm}")
+
+add_field("position_y", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("position_y")],
+ particle_type = True,
+ units=r"\rm{cm}")
+
+add_field("position_z", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("position_z")],
+ particle_type = True,
+ units=r"\rm{cm}")
+
+add_field("VEL", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("VEL")],
+ units=r"")
+
+add_field("id", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("ID")],
+ units=r"")
+
+add_field("mass", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("mass")],
+ units=r"\rm{g}")
+def _particle_mass(field, data):
+ return data["mass"]/just_one(data["CellVolume"])
+def _convert_particle_mass(data):
+ return 1.0
+add_field("particle_mass", function=_particle_mass, take_log=True,
+ convert_function=_convert_particle_mass,
+ validators = [ValidateSpatial(0)],
+ units=r"\mathrm{g}\/\mathrm{cm}^{-3}")
+
+add_field("U", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("U")],
+ units=r"")
+
+add_field("NE", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("NE")],
+ units=r"")
+
+add_field("POT", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("POT")],
+ units=r"")
+
+add_field("ACCE", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("ACCE")],
+ units=r"")
+
+add_field("ENDT", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("ENDT")],
+ units=r"")
+
+add_field("TSTP", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("TSTP")],
+ units=r"")
+
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/gadget/io.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/io.py
@@ -0,0 +1,55 @@
+"""
+Gadget-specific data-file handling function
+
+Author: Christopher E Moody <juxtaposicion at gmail.com>
+Affiliation: UC Santa Cruz
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Christopher E Moody, Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import h5py
+import numpy as np
+
+from yt.utilities.io_handler import \
+ BaseIOHandler
+
+class IOHandlerGadget(BaseIOHandler):
+ _data_style = 'gadget_infrastructure'
+ def _read_data_set(self, grid, field):
+ data = []
+ fh = h5py.File(grid.filename,mode='r')
+ for ptype in grid.particle_types:
+ address = '/data/grid_%010i/particles/%s/%s' % (grid.id, ptype, field)
+ data.append(fh[address][:])
+ if len(data) > 0:
+ data = np.concatenate(data)
+ fh.close()
+ return np.array(data)
+ def _read_field_names(self,grid):
+ adr = grid.Address
+ fh = h5py.File(grid.filename,mode='r')
+ rets = cPickle.loads(fh['/root'].attrs['fieldnames'])
+ fh.close()
+ return rets
+
+ def _read_data_slice(self,grid, field, axis, coord):
+ #how would we implement axis here?
+ dat = self._read_data_set(grid,field)
+ return dat[coord]
+
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/gadget/setup.py
--- /dev/null
+++ b/yt/frontends/sph/gadget/setup.py
@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+import setuptools
+import os
+import sys
+import os.path
+
+
+def configuration(parent_package='', top_path=None):
+ from numpy.distutils.misc_util import Configuration
+ config = Configuration('gadget', parent_package, top_path)
+ config.make_config_py() # installs __config__.py
+ #config.make_svn_version_py()
+ return config
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/io.py
--- /dev/null
+++ b/yt/frontends/sph/io.py
@@ -0,0 +1,116 @@
+"""
+Gadget-specific data-file handling function
+
+Author: Christopher E Moody <juxtaposicion at gmail.com>
+Affiliation: UC Santa Cruz
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2010-2011 Christopher E Moody, Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import h5py
+import numpy as np
+from yt.funcs import *
+from yt.utilities.exceptions import *
+
+from yt.utilities.io_handler import \
+ BaseIOHandler
+
+_vector_fields = ("Coordinates", "Velocity")
+
+class IOHandlerOWLS(BaseIOHandler):
+ _data_style = "OWLS"
+
+ def _read_fluid_selection(self, chunks, selector, fields, size):
+ raise NotImplementedError
+
+ def _read_particle_selection(self, chunks, selector, fields):
+ rv = {}
+ # We first need a set of masks for each particle type
+ ptf = defaultdict(list)
+ psize = defaultdict(lambda: 0)
+ chunks = list(chunks)
+ for ftype, fname in fields:
+ ptf[ftype].append(fname)
+ for chunk in chunks: # Will be OWLS domains
+ for subset in chunk.objs:
+ for ptype, field_list in sorted(ptf.items()):
+ f = h5py.File(subset.domain.domain_filename, "r")
+ coords = f["/%s/Coordinates" % ptype][:].astype("float64")
+ psize[ptype] += selector.count_points(
+ coords[:,0], coords[:,1], coords[:,2])
+ del coords
+ f.close()
+ # Now we have all the sizes, and we can allocate
+ ind = {}
+ for field in fields:
+ mylog.debug("Allocating %s values for %s", psize[field[0]], field)
+ if field[1] in _vector_fields:
+ shape = (psize[field[0]], 3)
+ else:
+ shape = psize[field[0]]
+ rv[field] = np.empty(shape, dtype="float64")
+ ind[field] = 0
+ for chunk in chunks: # Will be OWLS domains
+ for subset in chunk.objs:
+ for ptype, field_list in sorted(ptf.items()):
+ f = h5py.File(subset.domain.domain_filename, "r")
+ g = f["/%s" % ptype]
+ coords = g["Coordinates"][:].astype("float64")
+ mask = selector.select_points(
+ coords[:,0], coords[:,1], coords[:,2])
+ del coords
+ if mask is None: continue
+ for field in field_list:
+ data = g[field][mask,...]
+ my_ind = ind[ptype, field]
+ mylog.debug("Filling from %s to %s with %s",
+ my_ind, my_ind+data.shape[0], field)
+ rv[ptype, field][my_ind:my_ind + data.shape[0],...] = data
+ ind[ptype, field] += data.shape[0]
+ f.close()
+ return rv
+
+ def _initialize_octree(self, domain, octree):
+ f = h5py.File(domain.domain_filename, "r")
+ for key in f.keys():
+ if not key.startswith("PartType"): continue
+ pos = f[key]["Coordinates"][:].astype("float64")
+ octree.add(pos, domain.domain_number)
+ f.close()
+
+ def _count_particles(self, domain_filename):
+ f = h5py.File(domain_filename, "r")
+ npart = f["/Header"].attrs["NumPart_ThisFile"].sum()
+ f.close()
+ return npart
+
+ def _identify_fields(self, domain_filename):
+ f = h5py.File(domain_filename, "r")
+ fields = []
+ for key in f.keys():
+ if not key.startswith("PartType"): continue
+ g = f[key]
+ #ptype = int(key[8:])
+ ptype = str(key)
+ for k in g.keys():
+ if not hasattr(g[k], "shape"): continue
+ # str => not unicode!
+ fields.append((ptype, str(k)))
+ f.close()
+ return fields
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/frontends/sph/setup.py
--- /dev/null
+++ b/yt/frontends/sph/setup.py
@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+import setuptools
+import os
+import sys
+import os.path
+
+
+def configuration(parent_package='', top_path=None):
+ from numpy.distutils.misc_util import Configuration
+ config = Configuration('gadget', parent_package, top_path)
+ config.make_config_py() # installs __config__.py
+ #config.make_svn_version_py()
+ return config
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/geometry/geometry_handler.py
--- a/yt/geometry/geometry_handler.py
+++ b/yt/geometry/geometry_handler.py
@@ -59,12 +59,12 @@
mylog.debug("Setting up classes.")
self._setup_classes()
+ mylog.debug("Initializing data grid data IO")
+ self._setup_data_io()
+
mylog.debug("Setting up domain geometry.")
self._setup_geometry()
- mylog.debug("Initializing data grid data IO")
- self._setup_data_io()
-
mylog.debug("Detecting fields.")
self._detect_fields()
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/geometry/oct_container.pxd
--- a/yt/geometry/oct_container.pxd
+++ b/yt/geometry/oct_container.pxd
@@ -24,6 +24,7 @@
"""
cimport numpy as np
+from fp_utils cimport *
cdef struct ParticleArrays
@@ -64,5 +65,4 @@
Oct *oct
ParticleArrays *next
np.float64_t **pos
- np.int64_t *domain_id
np.int64_t np
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -23,7 +23,7 @@
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
-from libc.stdlib cimport malloc, free
+from libc.stdlib cimport malloc, free, qsort
from libc.math cimport floor
cimport numpy as np
import numpy as np
@@ -146,7 +146,7 @@
for i in range(3):
pp[i] = ppos[i] - self.DLE[i]
dds[i] = (self.DRE[i] - self.DLE[i])/self.nn[i]
- ind[i] = <np.int64_t> (floor(pp[i]/dds[i]))
+ ind[i] = <np.int64_t> ((pp[i] - self.DLE[i])/dds[i])
cp[i] = (ind[i] + 0.5) * dds[i]
cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
while cur.children[0][0][0] != NULL:
@@ -543,10 +543,21 @@
local_filled += 1
return local_filled
+cdef public int compare_octs(void *vo1, void *vo2) nogil:
+ cdef Oct *o1 = (<Oct**> vo1)[0]
+ cdef Oct *o2 = (<Oct**> vo2)[0]
+ if o1.domain < o2.domain: return -1
+ elif o1.domain == o2.domain:
+ if o1.level < o2.level: return -1
+ if o1.level > o2.level: return 1
+ else: return 0
+ elif o1.domain > o2.domain: return 1
+
cdef class ParticleOctreeContainer(OctreeContainer):
cdef ParticleArrays *first_sd
cdef ParticleArrays *last_sd
cdef Oct** oct_list
+ cdef np.int64_t *dom_offsets
def allocate_root(self):
cdef int i, j, k
@@ -568,6 +579,8 @@
for j in range(self.nn[1]):
for k in range(self.nn[2]):
self.visit_free(self.root_mesh[i][j][k])
+ free(self.oct_list)
+ free(self.dom_offsets)
cdef void visit_free(self, Oct *o):
cdef int i, j, k
@@ -581,7 +594,6 @@
for i in range(3):
free(o.sd.pos[i])
free(o.sd.pos)
- free(o.sd.domain_id)
free(o)
@cython.boundscheck(False)
@@ -670,7 +682,7 @@
def finalize(self):
self.oct_list = <Oct**> malloc(sizeof(Oct*)*self.nocts)
- cdef i = 0
+ cdef np.int64_t i = 0
cdef ParticleArrays *c = self.first_sd
while c != NULL:
self.oct_list[i] = c.oct
@@ -679,9 +691,20 @@
free(c.pos[j])
free(c.pos)
c.pos = NULL
- # We should also include a shortening of the domain IDs here
c = c.next
i += 1
+ qsort(self.oct_list, self.nocts, sizeof(Oct*), &compare_octs)
+ cdef int cur_dom = -1
+ # We always need at least 2, and if max_domain is 0, we need 3.
+ self.dom_offsets = <np.int64_t *>malloc(sizeof(np.int64_t) *
+ (self.max_domain + 3))
+ self.dom_offsets[0] = 0
+ for i in range(self.nocts):
+ self.oct_list[i].local_ind = i
+ if self.oct_list[i].domain > cur_dom:
+ cur_dom = self.oct_list[i].domain
+ self.dom_offsets[cur_dom + 1] = i
+ self.dom_offsets[cur_dom + 2] = self.nocts
cdef Oct* allocate_oct(self):
self.nocts += 1
@@ -706,13 +729,11 @@
self.last_sd = sd
sd.oct = my_oct
sd.next = NULL
- sd.domain_id = <np.int64_t *> malloc(sizeof(np.int64_t) * 32)
sd.pos = <np.float64_t **> malloc(sizeof(np.float64_t*) * 3)
for i in range(3):
sd.pos[i] = <np.float64_t *> malloc(sizeof(np.float64_t) * 32)
for i in range(32):
sd.pos[0][i] = sd.pos[1][i] = sd.pos[2][i] = 0.0
- sd.domain_id[i] = -1
sd.np = 0
return my_oct
@@ -724,24 +745,45 @@
c = c.next
return total
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def count_levels(self, int max_level, int domain_id,
+ np.ndarray[np.uint8_t, ndim=2, cast=True] mask):
+ cdef np.ndarray[np.int64_t, ndim=1] level_count
+ cdef Oct *o
+ cdef int oi, i
+ level_count = np.zeros(max_level+1, 'int64')
+ cdef np.int64_t ndo, doff
+ ndo = self.dom_offsets[domain_id + 2] \
+ - self.dom_offsets[domain_id + 1]
+ doff = self.dom_offsets[domain_id + 1]
+ for oi in range(ndo):
+ o = self.oct_list[oi + doff]
+ for i in range(8):
+ if mask[o.local_ind, i] == 0: continue
+ level_count[o.level] += 1
+ return level_count
+
def add(self, np.ndarray[np.float64_t, ndim=2] pos, np.int64_t domain_id):
cdef int no = pos.shape[0]
cdef int p, i, level
cdef np.float64_t dds[3], cp[3], pp[3]
cdef int ind[3]
self.max_domain = max(self.max_domain, domain_id)
+ cdef int mid, mad
if self.root_mesh[0][0][0] == NULL: self.allocate_root()
for p in range(no):
level = 0
for i in range(3):
pp[i] = pos[p, i]
dds[i] = (self.DRE[i] + self.DLE[i])/self.nn[i]
- ind[i] = <np.int64_t> (pp[i]/dds[i])
+ ind[i] = <np.int64_t> ((pp[i] - self.DLE[i])/dds[i])
cp[i] = (ind[i] + 0.5) * dds[i]
cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
if cur == NULL:
raise RuntimeError
- if cur.sd.np == 32:
+ if self._check_refine(cur, cp, domain_id) == 1:
self.refine_oct(cur, cp)
while cur.sd.np < 0:
for i in range(3):
@@ -754,16 +796,25 @@
cp[i] += dds[i]/2.0
cur = cur.children[ind[0]][ind[1]][ind[2]]
level += 1
- if cur.sd.np == 32:
+ if self._check_refine(cur, cp, domain_id) == 1:
self.refine_oct(cur, cp)
# Now we copy in our particle
pi = cur.sd.np
cur.level = level
for i in range(3):
cur.sd.pos[i][pi] = pp[i]
- cur.sd.domain_id[pi] = domain_id
+ cur.domain = domain_id
cur.sd.np += 1
+ cdef int _check_refine(self, Oct *cur, np.float64_t cp[3], int domain_id):
+ if cur.children[0][0][0] != NULL:
+ return 0
+ elif cur.sd.np == 32:
+ return 1
+ elif cur.domain >= 0 and cur.domain != domain_id:
+ return 1
+ return 0
+
cdef void refine_oct(self, Oct *o, np.float64_t pos[3]):
cdef int i, j, k, m, ind[3]
cdef Oct *noct
@@ -777,7 +828,7 @@
noct.pos[2] = (o.pos[2] << 1) + k
noct.parent = o
o.children[i][j][k] = noct
- for m in range(32):
+ for m in range(o.sd.np):
for i in range(3):
if o.sd.pos[i][m] < pos[i]:
ind[i] = 0
@@ -787,12 +838,12 @@
k = noct.sd.np
for i in range(3):
noct.sd.pos[i][k] = o.sd.pos[i][m]
- noct.sd.domain_id[k] = o.sd.domain_id[k]
+ noct.domain = o.domain
noct.sd.np += 1
o.sd.np = -1
+ o.domain = -1
for i in range(3):
free(o.sd.pos[i])
- free(o.sd.domain_id)
free(o.sd.pos)
def recursively_count(self):
@@ -828,14 +879,13 @@
for oi in range(self.nocts):
m = 0
o = self.oct_list[oi]
- if o.sd.np <= 0: continue
+ if o.sd.np <= 0 or o.domain == -1: continue
for i in range(8):
if mask[oi, i] == 1:
m = 1
break
if m == 0: continue
- for i in range(o.sd.np):
- dmask[o.sd.domain_id[i]] = 1
+ dmask[o.domain] = 1
return dmask.astype("bool")
@cython.boundscheck(False)
@@ -852,3 +902,25 @@
tnp += neighbors[i].sd.np
return tnp
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def count_cells(self, SelectorObject selector,
+ np.ndarray[np.uint8_t, ndim=2, cast=True] mask):
+ cdef int i, j, k, oi
+ # pos here is CELL center, not OCT center.
+ cdef np.float64_t pos[3]
+ cdef int n = mask.shape[0]
+ cdef np.float64_t base_dx[3], dx[3]
+ cdef int eterm[3]
+ cdef np.ndarray[np.int64_t, ndim=1] count
+ count = np.zeros(self.max_domain + 1, 'int64')
+ print count.shape[0], mask.shape[0], mask.shape[1], self.nocts
+ for oi in range(n):
+ o = self.oct_list[oi]
+ if o.domain == -1: continue
+ for i in range(8):
+ count[o.domain] += mask[oi,i]
+ return count
+
+
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/geometry/setup.py
--- a/yt/geometry/setup.py
+++ b/yt/geometry/setup.py
@@ -9,6 +9,7 @@
config = Configuration('geometry',parent_package,top_path)
config.add_extension("oct_container",
["yt/geometry/oct_container.pyx"],
+ include_dirs=["yt/utilities/lib/"],
libraries=["m"],
depends=["yt/utilities/lib/fp_utils.pxd",
"yt/geometry/oct_container.pxd",
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/geometry/tests/test_particle_octree.py
--- /dev/null
+++ b/yt/geometry/tests/test_particle_octree.py
@@ -0,0 +1,39 @@
+from yt.testing import *
+import numpy as np
+from yt.geometry.oct_container import ParticleOctreeContainer
+import time, os
+
+NPART = 32**3
+NDIM = 64
+DLE = np.array([0.0, 0.0, 0.0])
+DRE = np.array([10.0, 10.0, 10.0])
+
+def test_add_particles_random():
+ np.random.seed(int(0x4d3d3d3))
+ pos = np.random.normal(0.5, scale=0.05, size=(NPART,3)) * (DRE-DLE) + DLE
+ for i in range(3):
+ np.clip(pos[:,i], DLE[i], DRE[i], pos[:,i])
+ for ndom in [1, 2, 4, 8]:
+ octree = ParticleOctreeContainer((NDIM, NDIM, NDIM), DLE, DRE)
+ for dom in range(ndom):
+ octree.add(pos[dom::ndom,:], dom)
+ octree.finalize()
+ # This visits every oct.
+ lin_count = octree.linearly_count()
+ tc = octree.recursively_count()
+ total_count = np.zeros(len(tc), dtype="int32")
+ for i in sorted(tc):
+ total_count[i] = tc[i]
+ yield assert_equal, lin_count, total_count.sum()
+ mask = np.ones((total_count.sum(), 8), dtype="bool")
+ # This visits every cell -- including those covered by octs.
+ level_count = octree.count_levels(total_count.size-1, -1, mask)
+ for dom in range(ndom):
+ level_count += octree.count_levels(total_count.size-1, dom, mask)
+ yield assert_equal, level_count[0], NDIM**3 * 8
+ yield assert_equal, level_count, total_count * 8
+
+if __name__=="__main__":
+ for i in test_add_particles_random():
+ i[0](*i[1:])
+ time.sleep(1)
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/utilities/exceptions.py
--- a/yt/utilities/exceptions.py
+++ b/yt/utilities/exceptions.py
@@ -208,3 +208,10 @@
def __str__(self):
return "Enzo test output file (OutputLog) not generated for: " + \
"'%s'" % (self.testname) + ".\nTest did not complete."
+
+class YTFieldNotParseable(YTException):
+ def __init__(self, field):
+ self.field = field
+
+ def __str__(self):
+ return "Cannot identify field %s" % self.field
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r d000eb4f9e52fccc2b8277e97156018ae3165e4d yt/utilities/voropp.pyx
--- a/yt/utilities/voropp.pyx
+++ b/yt/utilities/voropp.pyx
@@ -31,19 +31,33 @@
cimport numpy as np
cimport cython
-cdef extern from "voro++.cc":
+cdef extern from "voro++.hh" namespace "voro":
+ cdef cppclass c_loop_all
+
+ cdef cppclass voronoicell:
+ double volume()
+
cdef cppclass container:
container(double xmin, double xmax, double ymin, double ymax,
double zmin, double zmax, int nx, int ny, int nz,
libcpp.bool xper, libcpp.bool yper, libcpp.bool zper, int alloc)
void put(int n, double x, double y, double z)
void store_cell_volumes(double *vols)
+ int compute_cell(voronoicell c, c_loop_all vl)
+ double sum_cell_volumes()
+
+ cdef cppclass c_loop_all:
+ c_loop_all(container &con)
+ int inc()
+ int start()
cdef class VoronoiVolume:
cdef container *my_con
- cdef int npart
- def __init__(self, xi, yi, zi):
- self.my_con = new container(0.0, 1.0, 0.0, 1.0, 0.0, 1.0,
+ cdef public int npart
+ def __init__(self, xi, yi, zi, left_edge, right_edge):
+ self.my_con = new container(left_edge[0], right_edge[0],
+ left_edge[1], right_edge[1],
+ left_edge[2], right_edge[2],
xi, yi, zi, False, False, False, 8)
self.npart = 0
@@ -65,5 +79,15 @@
def get_volumes(self):
cdef np.ndarray vol = np.zeros(self.npart, 'double')
cdef double *vdouble = <double *> vol.data
- self.my_con.store_cell_volumes(vdouble)
+ #self.my_con.store_cell_volumes(vdouble)
+ cdef c_loop_all *vl = new c_loop_all(deref(self.my_con))
+ cdef voronoicell c
+ if not vl.start(): return
+ cdef int i = 0
+ while 1:
+ if self.my_con.compute_cell(c, deref(vl)):
+ vol[i] = c.volume()
+ if not vl.inc(): break
+ i += 1
+ del vl
return vol
Repository URL: https://bitbucket.org/yt_analysis/yt-3.0/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list