[yt-svn] commit/yt-3.0: 14 new changesets
Bitbucket
commits-noreply at bitbucket.org
Thu Jan 10 07:09:02 PST 2013
14 new commits in yt-3.0:
https://bitbucket.org/yt_analysis/yt-3.0/commits/4cd73bea3369/
changeset: 4cd73bea3369
branch: yt-3.0
user: MatthewTurk
date: 2013-01-08 20:54:20
summary: Allowing select_grids to use min_level and max_level attributes.
In principle this should be portable between selector objects, but the first
method that will use it will be the covering grids.
affected #: 3 files
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r 4cd73bea3369035e2a1a6291ee1664cae795e089 yt/geometry/grid_geometry_handler.py
--- a/yt/geometry/grid_geometry_handler.py
+++ b/yt/geometry/grid_geometry_handler.py
@@ -197,7 +197,8 @@
dobj._chunk_info[0] = dobj
elif getattr(dobj, "_grids", None) is None:
gi = dobj.selector.select_grids(self.grid_left_edge,
- self.grid_right_edge)
+ self.grid_right_edge,
+ self.grid_levels)
grids = list(sorted(self.grids[gi], key = lambda g: g.filename))
dobj._chunk_info = np.empty(len(grids), dtype='object')
for i, g in enumerate(grids):
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r 4cd73bea3369035e2a1a6291ee1664cae795e089 yt/geometry/selection_routines.pxd
--- a/yt/geometry/selection_routines.pxd
+++ b/yt/geometry/selection_routines.pxd
@@ -28,12 +28,16 @@
cdef struct Oct
cdef class SelectorObject:
+ cdef public np.int32_t min_level
+ cdef public np.int32_t max_level
+
cdef void recursively_select_octs(self, Oct *root,
np.float64_t pos[3], np.float64_t dds[3],
np.ndarray[np.uint8_t, ndim=2] mask,
int level = ?)
cdef int select_grid(self, np.float64_t left_edge[3],
- np.float64_t right_edge[3]) nogil
+ np.float64_t right_edge[3],
+ np.int32_t level) nogil
cdef int select_cell(self, np.float64_t pos[3], np.float64_t dds[3],
int eterm[3]) nogil
cdef void set_bounds(self,
diff -r aaaca172321776fadf448c18ce8753dee31e3999 -r 4cd73bea3369035e2a1a6291ee1664cae795e089 yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -118,12 +118,17 @@
cdef class SelectorObject:
+ def __cinit__(self, dobj):
+ self.min_level = getattr(dobj, "min_level", 0)
+ self.max_level = getattr(dobj, "max_level", 99)
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def select_grids(self,
np.ndarray[np.float64_t, ndim=2] left_edges,
- np.ndarray[np.float64_t, ndim=2] right_edges):
+ np.ndarray[np.float64_t, ndim=2] right_edges,
+ np.ndarray[np.int32_t, ndim=2] levels):
cdef int i, n
cdef int ng = left_edges.shape[0]
cdef np.ndarray[np.uint8_t, ndim=1] gridi = np.zeros(ng, dtype='uint8')
@@ -135,7 +140,7 @@
for i in range(3):
LE[i] = left_edges[n, i]
RE[i] = right_edges[n, i]
- gridi[n] = self.select_grid(LE, RE)
+ gridi[n] = self.select_grid(LE, RE, levels[n,0])
return gridi.astype("bool")
@cython.boundscheck(False)
@@ -158,7 +163,7 @@
if octree.root_mesh[i][j][k] == NULL: continue
self.recursively_select_octs(
octree.root_mesh[i][j][k],
- pos, dds, mask)
+ pos, dds, mask, 0) # last is level
pos[2] += dds[2]
pos[1] += dds[2]
pos[0] += dds[2]
@@ -182,7 +187,7 @@
LE[i] = pos[i] - dds[i]/2.0
RE[i] = pos[i] + dds[i]/2.0
#print LE[0], RE[0], LE[1], RE[1], LE[2], RE[2]
- res = self.select_grid(LE, RE)
+ res = self.select_grid(LE, RE, level)
cdef int eterm[3]
eterm[0] = eterm[1] = eterm[2] = 0
if res == 0:
@@ -192,6 +197,15 @@
# Now we visit all our children. We subtract off sdds for the first
# pass because we center it on the first cell.
spos[0] = pos[0] - sdds[0]/2.0
+ cdef int next_level, this_level
+ # next_level: an int that says whether or not we can progress to children
+ # this_level: an int that says whether or not we can select from this
+ # level
+ next_level = this_level = 1
+ if level == self.max_level:
+ next_level = 0
+ if level < self.min_level or level > self.max_level:
+ this_level = 0
for i in range(2):
spos[1] = pos[1] - sdds[1]/2.0
for j in range(2):
@@ -199,11 +213,11 @@
for k in range(2):
ii = ((k*2)+j)*2+i
ch = root.children[i][j][k]
- if ch != NULL:
+ if this_level == 1 and ch != NULL:
mask[root.local_ind, ii] = 0
self.recursively_select_octs(
ch, spos, sdds, mask, level + 1)
- else:
+ elif this_level == 1:
mask[root.local_ind, ii] = \
self.select_cell(spos, sdds, eterm)
spos[2] += sdds[2]
@@ -211,7 +225,8 @@
spos[0] += sdds[0]
cdef int select_grid(self, np.float64_t left_edge[3],
- np.float64_t right_edge[3]) nogil:
+ np.float64_t right_edge[3],
+ np.int32_t level) nogil:
return 0
cdef int select_cell(self, np.float64_t pos[3], np.float64_t dds[3],
@@ -235,6 +250,14 @@
dds[i] = odds[i]
cdef int count = 0
cdef int check = 1
+ # Check for the level bounds
+ cdef np.int32_t level = gobj.Level
+ if level < self.min_level or level > self.max_level:
+ return count
+ # We set this to 1 if we ignore child_mask
+ cdef int this_level = 0
+ if level == self.max_level:
+ this_level = 1
cdef int eterm[3]
self.set_bounds(<np.float64_t *> left_edge.data,
<np.float64_t *> right_edge.data,
@@ -250,7 +273,7 @@
pos[2] = left_edge[2] + dds[2] * 0.5
for k in range(ind[2][0], ind[2][1]):
eterm[2] = 0
- if child_mask[i,j,k] == 1:
+ if child_mask[i,j,k] == 1 or this_level == 1:
count += self.select_cell(pos, dds, eterm)
if eterm[2] == 1: break
pos[2] += dds[1]
@@ -262,7 +285,8 @@
for i in range(ind[0][0], ind[0][1]):
for j in range(ind[1][0], ind[1][1]):
for k in range(ind[2][0], ind[2][1]):
- count += child_mask[i,j,k]
+ if child_mask[i,j,k] == 1 or this_level == 1:
+ count += 1
return count
@cython.boundscheck(False)
@@ -290,6 +314,14 @@
<np.float64_t *> right_edge.data,
dds, ind, &check)
cdef int temp
+ # Check for the level bounds
+ cdef np.int32_t level = gobj.Level
+ if level < self.min_level or level > self.max_level:
+ return mask
+ # We set this to 1 if we ignore child_mask
+ cdef int this_level = 0
+ if level == self.max_level:
+ this_level = 1
with nogil:
if check == 1:
pos[0] = left_edge[0] + dds[0] * 0.5
@@ -301,7 +333,7 @@
pos[2] = left_edge[2] + dds[2] * 0.5
for k in range(ind[2][0], ind[2][1]):
eterm[2] = 0
- if child_mask[i,j,k] == 1:
+ if child_mask[i,j,k] == 1 or this_level == 1:
mask[i,j,k] = self.select_cell(pos, dds, eterm)
total += mask[i,j,k]
if eterm[2] == 1: break
@@ -314,7 +346,8 @@
for i in range(ind[0][0], ind[0][1]):
for j in range(ind[1][0], ind[1][1]):
for k in range(ind[2][0], ind[2][1]):
- mask[i,j,k] = child_mask[i,j,k]
+ if child_mask[i,j,k] == 1 or this_level == 1:
+ mask[i,j,k] = 1
total += mask[i,j,k]
if total == 0: return None
return mask.astype("bool")
@@ -381,7 +414,8 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int select_grid(self, np.float64_t left_edge[3],
- np.float64_t right_edge[3]) nogil:
+ np.float64_t right_edge[3],
+ np.int32_t level) nogil:
cdef np.float64_t box_center, relcenter, closest, dist, edge
return 1
cdef int id
@@ -432,7 +466,9 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int select_grid(self, np.float64_t left_edge[3],
- np.float64_t right_edge[3]) nogil:
+ np.float64_t right_edge[3],
+ np.int32_t level) nogil:
+ if level < self.min_level or level > self.max_level: return 0
for i in range(3):
if left_edge[i] >= self.right_edge[i]: return 0
if right_edge[i] <= self.left_edge[i]: return 0
@@ -495,7 +531,8 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int select_grid(self, np.float64_t left_edge[3],
- np.float64_t right_edge[3]) nogil:
+ np.float64_t right_edge[3],
+ np.int32_t level) nogil:
cdef np.float64_t *arr[2]
cdef np.float64_t pos[3], H, D, R2, temp
cdef int i, j, k, n
@@ -562,7 +599,8 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int select_grid(self, np.float64_t left_edge[3],
- np.float64_t right_edge[3]) nogil:
+ np.float64_t right_edge[3],
+ np.int32_t level) nogil:
cdef int i, j, k, n
cdef np.float64_t *arr[2]
cdef np.float64_t pos[3]
@@ -629,7 +667,8 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int select_grid(self, np.float64_t left_edge[3],
- np.float64_t right_edge[3]) nogil:
+ np.float64_t right_edge[3],
+ np.int32_t level) nogil:
if right_edge[self.axis] > self.coord \
and left_edge[self.axis] <= self.coord:
return 1
@@ -666,7 +705,8 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int select_grid(self, np.float64_t left_edge[3],
- np.float64_t right_edge[3]) nogil:
+ np.float64_t right_edge[3],
+ np.int32_t level) nogil:
if ( (self.px >= left_edge[self.px_ax])
and (self.px < right_edge[self.px_ax])
and (self.py >= left_edge[self.py_ax])
@@ -744,7 +784,8 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int select_grid(self, np.float64_t left_edge[3],
- np.float64_t right_edge[3]) nogil:
+ np.float64_t right_edge[3],
+ np.int32_t level) nogil:
cdef int i, ax
cdef int i1, i2
cdef np.float64_t vs[3], t, v[3]
@@ -885,7 +926,8 @@
@cython.cdivision(True)
def select_grids(self,
np.ndarray[np.float64_t, ndim=2] left_edges,
- np.ndarray[np.float64_t, ndim=2] right_edges):
+ np.ndarray[np.float64_t, ndim=2] right_edges,
+ np.ndarray[np.int32_t, ndim=2] levels):
cdef int i, n
cdef int ng = left_edges.shape[0]
cdef np.ndarray[np.uint8_t, ndim=1] gridi = np.zeros(ng, dtype='uint8')
@@ -933,7 +975,8 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int select_grid(self, np.float64_t left_edge[3],
- np.float64_t right_edge[3]) nogil:
+ np.float64_t right_edge[3],
+ np.int32_t level) nogil:
# This is the sphere selection
cdef np.float64_t radius2, box_center, relcenter, closest, dist, edge
return 1
@@ -996,7 +1039,8 @@
@cython.cdivision(True)
def select_grids(self,
np.ndarray[np.float64_t, ndim=2] left_edges,
- np.ndarray[np.float64_t, ndim=2] right_edges):
+ np.ndarray[np.float64_t, ndim=2] right_edges,
+ np.ndarray[np.int32_t, ndim=2] levels):
cdef int ng = left_edges.shape[0]
cdef np.ndarray[np.uint8_t, ndim=1] gridi = np.zeros(ng, dtype='uint8')
gridi[self.ind] = 1
https://bitbucket.org/yt_analysis/yt-3.0/commits/5a01ae4eb720/
changeset: 5a01ae4eb720
branch: yt-3.0
user: MatthewTurk
date: 2013-01-08 22:46:36
summary: Adding a new fill_region routine (and tests) for flat position arrays.
affected #: 2 files
diff -r 4cd73bea3369035e2a1a6291ee1664cae795e089 -r 5a01ae4eb7201ff937ec80f8176944dfa8c8b33d yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -450,3 +450,39 @@
if ofield[i, j, k] == 1:
nfield[ni, nj, nk] = 1
return nfield.astype("bool")
+
+def fill_region(input_fields, output_fields,
+ np.int32_t output_level,
+ np.ndarray[np.int64_t, ndim=1] left_index,
+ np.ndarray[np.int64_t, ndim=2] ipos,
+ np.ndarray[np.int64_t, ndim=1] ires,
+ np.int64_t refine_by = 2
+ ):
+ cdef int i, n, skip
+ cdef np.int64_t iind[3], oind[3], dim[3], oi, oj, ok, rf
+ cdef np.ndarray[np.float64_t, ndim=3] ofield
+ cdef np.ndarray[np.float64_t, ndim=1] ifield
+ nf = len(input_fields)
+ for i in range(3):
+ dim[i] = output_fields[0].shape[i]
+ for n in range(nf):
+ ofield = output_fields[n]
+ ifield = input_fields[n]
+ for i in range(ipos.shape[0]):
+ for n in range(3):
+ iind[n] = ipos[i, n] - left_index[n]
+ rf = refine_by**(output_level - ires[i])
+ for oi in range(rf):
+ oind[0] = oi + iind[0] * rf
+ if oind[0] < left_index[0] or oind[0] >= left_index[0] + dim[0]:
+ continue
+ for oj in range(rf):
+ oind[1] = oj + iind[1] * rf
+ if oind[1] < left_index[1] or oind[1] >= left_index[1] + dim[1]:
+ continue
+ for ok in range(rf):
+ oind[2] = ok + iind[2] * rf
+ if oind[2] < left_index[2] or oind[2] >= left_index[2] + dim[2]:
+ continue
+ ofield[oind[0], oind[1], oind[2]] = \
+ ifield[i]
diff -r 4cd73bea3369035e2a1a6291ee1664cae795e089 -r 5a01ae4eb7201ff937ec80f8176944dfa8c8b33d yt/utilities/lib/tests/test_fill_region.py
--- /dev/null
+++ b/yt/utilities/lib/tests/test_fill_region.py
@@ -0,0 +1,29 @@
+from yt.testing import *
+from yt.utilities.lib import fill_region
+
+NDIM = 32
+
+def test_fill_region():
+ for level in range(2):
+ rf = 2**level
+ output_fields = [np.zeros((NDIM*rf,NDIM*rf,NDIM*rf), "float64")
+ for i in range(3)]
+ input_fields = [np.empty(NDIM**3, "float64")
+ for i in range(3)]
+ v = np.mgrid[0.0:1.0:NDIM*1j, 0.0:1.0:NDIM*1j, 0.0:1.0:NDIM*1j]
+ input_fields[0][:] = v[0].ravel()
+ input_fields[1][:] = v[1].ravel()
+ input_fields[2][:] = v[2].ravel()
+ left_index = np.zeros(3, "int64")
+ ipos = np.empty((NDIM**3, 3), dtype="int64")
+ ind = np.indices((NDIM,NDIM,NDIM))
+ ipos[:,0] = ind[0].ravel()
+ ipos[:,1] = ind[1].ravel()
+ ipos[:,2] = ind[2].ravel()
+ ires = np.zeros(NDIM*NDIM*NDIM, "int64")
+ fill_region(input_fields, output_fields, level,
+ left_index, ipos, ires)
+ for r in range(level + 1):
+ for o, i in zip(output_fields, v):
+ assert_equal( o[r::rf,r::rf,r::rf], i)
+
https://bitbucket.org/yt_analysis/yt-3.0/commits/8f40060e1c99/
changeset: 8f40060e1c99
branch: yt-3.0
user: MatthewTurk
date: 2013-01-09 05:11:52
summary: Covering grid works for many use cases.
Removed overlap projection (may return eventually) and the fixed res
projection.
Re-enabled covering grid tests. They do not all pass yet, as the covering grid
still needs to have NeedsOriginalGrid and generate_field implementation work.
affected #: 3 files
diff -r 5a01ae4eb7201ff937ec80f8176944dfa8c8b33d -r 8f40060e1c99d1fa45a328278beb33fd2e7963b5 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -42,7 +42,7 @@
from .field_info_container import \
NeedsOriginalGrid
from yt.utilities.lib import \
- QuadTree, ghost_zone_interpolate
+ QuadTree, ghost_zone_interpolate, fill_region
from yt.utilities.data_point_utilities import CombineGrids,\
DataCubeRefine, DataCubeReplace, FillRegion, FillBuffer
from yt.utilities.definitions import axis_names, x_dict, y_dict
@@ -364,382 +364,6 @@
pw = self._get_pw(fields, center, width, origin, axes_unit, 'Projection')
return pw
-class YTOverlapProjBase(YTSelectionContainer2D):
- _top_node = "/Projections"
- _key_fields = YTSelectionContainer2D._key_fields + ['weight_field']
- _type_name = "overlap_proj"
- _con_args = ('axis', 'weight_field')
- def __init__(self, field, axis, weight_field = None,
- max_level = None, center = None, pf = None,
- source=None, node_name = None, field_cuts = None,
- preload_style=None, serialize=True, **kwargs):
- """
- This is a data object corresponding to a line integral through the
- simulation domain.
-
- This object is typically accessed through the `proj` object that
- hangs off of hierarchy objects. AMRProj is a projection of a `field`
- along an `axis`. The field can have an associated `weight_field`, in
- which case the values are multiplied by a weight before being summed,
- and then divided by the sum of that weight; the two fundamental modes
- of operating are direct line integral (no weighting) and average along
- a line of sight (weighting.) Note also that lines of sight are
- integrated at every projected finest-level cell
-
- Parameters
- ----------
- field : string
- This is the field which will be "projected" along the axis. If
- multiple are specified (in a list) they will all be projected in
- the first pass.
- axis : int or axis
- The axis along which to slice. Can be 0, 1, or 2 for x, y, z.
- weight_field : string
- If supplied, the field being projected will be multiplied by this
- weight value before being integrated, and at the conclusion of the
- projection the resultant values will be divided by the projected
- `weight_field`.
- max_level : int
- If supplied, only cells at or below this level will be projected.
- center : array_like, optional
- The 'center' supplied to fields that use it. Note that this does
- not have to have `coord` as one value. Strictly optional.
- source : `yt.data_objects.api.YTDataContainer`, optional
- If specified, this will be the data source used for selecting
- regions to project.
- node_name: string, optional
- The node in the .yt file to find or store this slice at. Should
- probably not be used.
- field_cuts : list of strings, optional
- If supplied, each of these strings will be evaluated to cut a
- region of a grid out. They can be of the form "grid['Temperature']
- > 100" for instance.
- preload_style : string
- Either 'level', 'all', or None (default). Defines how grids are
- loaded -- either level by level, or all at once. Only applicable
- during parallel runs.
- serialize : bool, optional
- Whether we should store this projection in the .yt file or not.
- kwargs : dict of items
- Any additional values are passed as field parameters that can be
- accessed by generated fields.
-
- Examples
- --------
-
- >>> pf = load("RedshiftOutput0005")
- >>> proj = pf.h.proj("Density", "x")
- >>> print proj["Density"]
- """
- YTSelectionContainer2D.__init__(self, axis, field, pf, node_name = None, **kwargs)
- self.weight_field = weight_field
- self._field_cuts = field_cuts
- self.serialize = serialize
- self._set_center(center)
- if center is not None: self.set_field_parameter('center',center)
- self._node_name = node_name
- self._initialize_source(source)
- self._grids = self.source._grids
- if max_level == None:
- max_level = self.hierarchy.max_level
- if self.source is not None:
- max_level = min(max_level, self.source.grid_levels.max())
- self._max_level = max_level
- self._weight = weight_field
- self.preload_style = preload_style
- self.func = np.sum # for the future
- self.__retval_coords = {}
- self.__retval_fields = {}
- self.__retval_coarse = {}
- self.__overlap_masks = {}
- self._deserialize(node_name)
- self._refresh_data()
- if self._okay_to_serialize and self.serialize: self._serialize(node_name=self._node_name)
-
- def _convert_field_name(self, field):
- if field == "weight_field": return "weight_field_%s" % self._weight
- if field in self._key_fields: return field
- return "%s_%s" % (field, self._weight)
-
- def _initialize_source(self, source = None):
- if source is None:
- check, source = self.partition_hierarchy_2d(self.axis)
- self._check_region = check
- #self._okay_to_serialize = (not check)
- else:
- self._distributed = False
- self._okay_to_serialize = False
- self._check_region = True
- self.source = source
- if self._field_cuts is not None:
- # Override if field cuts are around; we don't want to serialize!
- self._check_region = True
- self._okay_to_serialize = False
- if self._node_name is not None:
- self._node_name = "%s/%s" % (self._top_node,self._node_name)
- self._okay_to_serialize = True
-
- def __calculate_overlap(self, level):
- s = self.source
- mylog.info("Generating overlap masks for level %s", level)
- i = 0
- pbar = get_pbar("Reading and masking grids ", len(s._grids))
- mylog.debug("Examining level %s", level)
- grids = s.select_grid_indices(level)
- RE = s.grid_right_edge[grids]
- LE = s.grid_left_edge[grids]
- for grid in s._grids[grids]:
- pbar.update(i)
- self.__overlap_masks[grid.id] = \
- grid._generate_overlap_masks(self.axis, LE, RE)
- i += 1
- pbar.finish()
- mylog.info("Finished calculating overlap.")
-
- def _get_dls(self, grid, fields):
- # Place holder for a time when maybe we will not be doing just
- # a single dx for every field.
- dls = []
- convs = []
- for field in fields + [self._weight]:
- if field is None: continue
- dls.append(just_one(grid['d%s' % axis_names[self.axis]]))
- convs.append(self.pf.units[self.pf.field_info[field].projection_conversion])
- return np.array(dls), np.array(convs)
-
- def __project_level(self, level, fields):
- grids_to_project = self.source.select_grids(level)
- dls, convs = self._get_dls(grids_to_project[0], fields)
- zero_out = (level != self._max_level)
- pbar = get_pbar('Projecting level % 2i / % 2i ' \
- % (level, self._max_level), len(grids_to_project))
- for pi, grid in enumerate(grids_to_project):
- g_coords, g_fields = self._project_grid(grid, fields, zero_out)
- self.__retval_coords[grid.id] = g_coords
- self.__retval_fields[grid.id] = g_fields
- for fi in range(len(fields)): g_fields[fi] *= dls[fi]
- if self._weight is not None: g_coords[3] *= dls[-1]
- pbar.update(pi)
- grid.clear_data()
- pbar.finish()
- self.__combine_grids_on_level(level) # In-place
- if level > 0 and level <= self._max_level:
- self.__refine_to_level(level) # In-place
- coord_data = []
- field_data = []
- for grid in grids_to_project:
- coarse = self.__retval_coords[grid.id][2]==0 # Where childmask = 0
- fine = ~coarse
- coord_data.append([pi[fine] for pi in self.__retval_coords[grid.id]])
- field_data.append([pi[fine] for pi in self.__retval_fields[grid.id]])
- self.__retval_coords[grid.id] = [pi[coarse] for pi in self.__retval_coords[grid.id]]
- self.__retval_fields[grid.id] = [pi[coarse] for pi in self.__retval_fields[grid.id]]
- coord_data = np.concatenate(coord_data, axis=1)
- field_data = np.concatenate(field_data, axis=1)
- if self._weight is not None:
- field_data = field_data / coord_data[3,:].reshape((1,coord_data.shape[1]))
- else:
- field_data *= convs[...,np.newaxis]
- mylog.info("Level %s done: %s final", \
- level, coord_data.shape[1])
- pdx = grids_to_project[0].dds[x_dict[self.axis]] # this is our dl
- pdy = grids_to_project[0].dds[y_dict[self.axis]] # this is our dl
- return coord_data, pdx, pdy, field_data
-
- def __combine_grids_on_level(self, level):
- grids = self.source.select_grids(level)
- grids_i = self.source.select_grid_indices(level)
- pbar = get_pbar('Combining level % 2i / % 2i ' \
- % (level, self._max_level), len(grids))
- # We have an N^2 check, so we try to be as quick as possible
- # and to skip as many as possible
- for pi, grid1 in enumerate(grids):
- pbar.update(pi)
- if self.__retval_coords[grid1.id][0].shape[0] == 0: continue
- for grid2 in self.source._grids[grids_i][self.__overlap_masks[grid1.id]]:
- if self.__retval_coords[grid2.id][0].shape[0] == 0 \
- or grid1.id == grid2.id:
- continue
- args = [] # First is source, then destination
- args += self.__retval_coords[grid2.id] + [self.__retval_fields[grid2.id]]
- args += self.__retval_coords[grid1.id] + [self.__retval_fields[grid1.id]]
- args.append(1) # Refinement factor
- args.append(np.ones(args[0].shape, dtype='int64'))
- kk = CombineGrids(*args)
- goodI = args[-1].astype('bool')
- self.__retval_coords[grid2.id] = \
- [coords[goodI] for coords in self.__retval_coords[grid2.id]]
- self.__retval_fields[grid2.id] = \
- [fields[goodI] for fields in self.__retval_fields[grid2.id]]
- pbar.finish()
-
- def __refine_to_level(self, level):
- grids = self.source.select_grids(level)
- grids_up = self.source.select_grid_indices(level - 1)
- pbar = get_pbar('Refining to level % 2i / % 2i ' \
- % (level, self._max_level), len(grids))
- for pi, grid1 in enumerate(grids):
- pbar.update(pi)
- for parent in ensure_list(grid1.Parent):
- if parent.id not in self.__overlap_masks: continue
- for grid2 in self.source._grids[grids_up][self.__overlap_masks[parent.id]]:
- if self.__retval_coords[grid2.id][0].shape[0] == 0: continue
- args = []
- args += self.__retval_coords[grid2.id] + [self.__retval_fields[grid2.id]]
- args += self.__retval_coords[grid1.id] + [self.__retval_fields[grid1.id]]
- # Refinement factor, which is same in all directions. Note
- # that this complicated rounding is because sometimes
- # epsilon differences in dds between the grids causes this
- # to round to up or down from the expected value.
- args.append(int(np.rint(grid2.dds / grid1.dds)[0]))
- args.append(np.ones(args[0].shape, dtype='int64'))
- kk = CombineGrids(*args)
- goodI = args[-1].astype('bool')
- self.__retval_coords[grid2.id] = \
- [coords[goodI] for coords in self.__retval_coords[grid2.id]]
- self.__retval_fields[grid2.id] = \
- [fields[goodI] for fields in self.__retval_fields[grid2.id]]
- for grid1 in self.source.select_grids(level-1):
- if not self._check_region and self.__retval_coords[grid1.id][0].size != 0:
- mylog.error("Something messed up, and %s still has %s points of data",
- grid1, self.__retval_coords[grid1.id][0].size)
- mylog.error("Please contact the yt-users mailing list.")
- raise ValueError(grid1, self.__retval_coords[grid1.id])
- pbar.finish()
-
- def get_data(self, fields):
- fields = ensure_list(fields)
- self._obtain_fields(fields, self._node_name)
- fields = [f for f in fields if f not in self.field_data]
- if len(fields) == 0: return
- coord_data = []
- field_data = []
- pdxs = []
- pdys = []
- # We do this here, but I am not convinced it should be done here
- # It is probably faster, as it consolidates IO, but if we did it in
- # _project_level, then it would be more memory conservative
- if self.preload_style == 'all':
- print "Preloading %s grids and getting %s" % (
- len(self.source._grids), self.get_dependencies(fields))
- self.comm.preload(self.source._grids,
- self.get_dependencies(fields), self.hierarchy.io)
- for level in range(0, self._max_level+1):
- if self.preload_style == 'level':
- self.comm.preload(self.source.select_grids(level),
- self.get_dependencies(fields), self.hierarchy.io)
- self.__calculate_overlap(level)
- my_coords, my_pdx, my_pdy, my_fields = \
- self.__project_level(level, fields)
- coord_data.append(my_coords)
- field_data.append(my_fields)
- pdxs.append(my_pdx * np.ones(my_coords.shape[1], dtype='float64'))
- pdys.append(my_pdx * np.ones(my_coords.shape[1], dtype='float64'))
- if self._check_region and False:
- check=self.__cleanup_level(level - 1)
- if len(check) > 0: all_data.append(check)
- # Now, we should clean up after ourselves...
- for grid in self.source.select_grids(level - 1):
- del self.__retval_coords[grid.id]
- del self.__retval_fields[grid.id]
- del self.__overlap_masks[grid.id]
- mylog.debug("End of projecting level level %s, memory usage %0.3e",
- level, get_memory_usage()/1024.)
- coord_data = np.concatenate(coord_data, axis=1)
- field_data = np.concatenate(field_data, axis=1)
- pdxs = np.concatenate(pdxs, axis=1)
- pdys = np.concatenate(pdys, axis=1)
- # We now convert to half-widths and center-points
- data = {}
- data['pdx'] = pdxs; del pdxs
- data['pdy'] = pdys; del pdys
- ox = self.pf.domain_left_edge[x_dict[self.axis]]
- oy = self.pf.domain_left_edge[y_dict[self.axis]]
- data['px'] = (coord_data[0,:]+0.5) * data['pdx'] + ox
- data['py'] = (coord_data[1,:]+0.5) * data['pdx'] + oy
- data['weight_field'] = coord_data[3,:].copy()
- del coord_data
- data['pdx'] *= 0.5
- data['pdy'] *= 0.5
- data['fields'] = field_data
- # Now we run the finalizer, which is ignored if we don't need it
- data = self.comm.par_combine_object(data, datatype='dict', op='cat')
- field_data = np.vsplit(data.pop('fields'), len(fields))
- for fi, field in enumerate(fields):
- self[field] = field_data[fi].ravel()
- if self.serialize: self._store_fields(field, self._node_name)
- for i in data.keys(): self[i] = data.pop(i)
- mylog.info("Projection completed")
-
- def add_fields(self, fields, weight = "CellMassMsun"):
- pass
-
- def _project_grid(self, grid, fields, zero_out):
- # We split this next bit into two sections to try to limit the IO load
- # on the system. This way, we perserve grid state (@restore_grid_state
- # in _get_data_from_grid *and* we attempt not to load weight data
- # independently of the standard field data.
- if self._weight is None:
- weight_data = np.ones(grid.ActiveDimensions, dtype='float64')
- if zero_out: weight_data[grid.child_indices] = 0
- masked_data = [fd.astype('float64') * weight_data
- for fd in self._get_data_from_grid(grid, fields)]
- else:
- fields_to_get = list(set(fields + [self._weight]))
- field_data = dict(zip(
- fields_to_get, self._get_data_from_grid(grid, fields_to_get)))
- weight_data = field_data[self._weight].copy().astype('float64')
- if zero_out: weight_data[grid.child_indices] = 0
- masked_data = [field_data[field].copy().astype('float64') * weight_data
- for field in fields]
- del field_data
- # if we zero it out here, then we only have to zero out the weight!
- full_proj = [self.func(field, axis=self.axis) for field in masked_data]
- weight_proj = self.func(weight_data, axis=self.axis)
- if (self._check_region and not self.source._is_fully_enclosed(grid)) or self._field_cuts is not None:
- used_data = self._get_points_in_region(grid).astype('bool')
- used_points = np.where(np.logical_or.reduce(used_data, self.axis))
- else:
- used_data = np.array([1.0], dtype='bool')
- used_points = slice(None)
- if zero_out:
- subgrid_mask = np.logical_and.reduce(
- np.logical_or(grid.child_mask,
- ~used_data),
- self.axis).astype('int64')
- else:
- subgrid_mask = np.ones(full_proj[0].shape, dtype='int64')
- xind, yind = [arr[used_points].ravel() for arr in np.indices(full_proj[0].shape)]
- start_index = grid.get_global_startindex()
- xpoints = (xind + (start_index[x_dict[self.axis]])).astype('int64')
- ypoints = (yind + (start_index[y_dict[self.axis]])).astype('int64')
- return ([xpoints, ypoints,
- subgrid_mask[used_points].ravel(),
- weight_proj[used_points].ravel()],
- [data[used_points].ravel() for data in full_proj])
-
- def _get_points_in_region(self, grid):
- pointI = self.source._get_point_indices(grid, use_child_mask=False)
- point_mask = np.zeros(grid.ActiveDimensions)
- point_mask[pointI] = 1.0
- if self._field_cuts is not None:
- for cut in self._field_cuts:
- point_mask *= eval(cut)
- return point_mask
-
- def _get_data_from_grid(self, grid, fields):
- fields = ensure_list(fields)
- if self._check_region:
- bad_points = self._get_points_in_region(grid)
- else:
- bad_points = 1.0
- return [grid[field] * bad_points for field in fields]
-
- def _gen_node_name(self):
- return "%s/%s" % \
- (self._top_node, self.axis)
-
-
class YTCoveringGridBase(YTSelectionContainer3D):
_spatial = True
_type_name = "covering_grid"
@@ -779,115 +403,26 @@
self.global_startindex = np.rint((self.left_edge-self.pf.domain_left_edge)/self.dds).astype('int64')
self.domain_width = np.rint((self.pf.domain_right_edge -
self.pf.domain_left_edge)/self.dds).astype('int64')
- self.get_data(fields)
+ self._data_source = self.pf.h.region(
+ self.center, self.left_edge, self.right_edge)
+ self._data_source.min_level = 0
+ self._data_source.max_level = self.level
- def _get_list_of_grids(self, buffer = 0.0):
- if self._grids is not None: return
- if np.any(self.left_edge - buffer < self.pf.domain_left_edge) or \
- np.any(self.right_edge + buffer > self.pf.domain_right_edge):
- grids,ind = self.pf.hierarchy.get_periodic_box_grids_below_level(
- self.left_edge - buffer,
- self.right_edge + buffer, self.level)
- else:
- grids,ind = self.pf.hierarchy.get_box_grids_below_level(
- self.left_edge - buffer,
- self.right_edge + buffer, self.level)
- sort_ind = np.argsort(self.pf.h.grid_levels.ravel()[ind])
- self._grids = self.pf.hierarchy.grids[ind][(sort_ind,)][::-1]
-
- def _refresh_data(self):
- YTSelectionContainer3D._refresh_data(self)
- self['dx'] = self.dds[0] * np.ones(self.ActiveDimensions, dtype='float64')
- self['dy'] = self.dds[1] * np.ones(self.ActiveDimensions, dtype='float64')
- self['dz'] = self.dds[2] * np.ones(self.ActiveDimensions, dtype='float64')
-
- def get_data(self, fields):
- if self._grids is None:
- self._get_list_of_grids()
- fields = ensure_list(fields)
- obtain_fields = []
- for field in fields:
- if self.field_data.has_key(field): continue
- if field not in self.hierarchy.field_list:
- try:
- #print "Generating", field
- self._generate_field(field)
- continue
- except NeedsOriginalGrid, ngt_exception:
- pass
- obtain_fields.append(field)
- self[field] = np.zeros(self.ActiveDimensions, dtype='float64') -999
- if len(obtain_fields) == 0: return
- mylog.debug("Getting fields %s from %s possible grids",
- obtain_fields, len(self._grids))
- if self._use_pbar: pbar = \
- get_pbar('Searching grids for values ', len(self._grids))
- count = self.ActiveDimensions.prod()
- for i, grid in enumerate(self._grids):
- if self._use_pbar: pbar.update(i)
- count -= self._get_data_from_grid(grid, obtain_fields)
- if count <= 0: break
- if self._use_pbar: pbar.finish()
- if count > 0 or np.any(self[obtain_fields[0]] == -999):
- # and self.dx < self.hierarchy.grids[0].dx:
- n_bad = np.where(self[obtain_fields[0]]==-999)[0].size
- mylog.error("Covering problem: %s cells are uncovered", n_bad)
- raise KeyError(n_bad)
-
- def _generate_field(self, field):
- if self.pf.field_info.has_key(field):
- # First we check the validator; this might even raise!
- self.pf.field_info[field].check_available(self)
- self[field] = self.pf.field_info[field](self)
- else: # Can't find the field, try as it might
- raise KeyError(field)
-
- def flush_data(self, fields=None):
- """
- Any modifications made to the data in this object are pushed back
- to the originating grids, except the cells where those grids are both
- below the current level `and` have child cells.
- """
- self._get_list_of_grids()
- # We don't generate coordinates here.
- for grid in self._grids:
- self._flush_data_to_grid(grid, ensure_list(fields))
-
- def _get_data_from_grid(self, grid, fields):
- ll = int(grid.Level == self.level)
- ref_ratio = self.pf.refine_by**(self.level - grid.Level)
- g_fields = [grid[field].astype("float64") for field in fields]
- c_fields = [self[field] for field in fields]
- count = FillRegion(ref_ratio,
- grid.get_global_startindex(), self.global_startindex,
- c_fields, g_fields,
- self.ActiveDimensions, grid.ActiveDimensions,
- grid.child_mask, self.domain_width, ll, 0)
- return count
-
- def _flush_data_to_grid(self, grid, fields):
- ll = int(grid.Level == self.level)
- ref_ratio = self.pf.refine_by**(self.level - grid.Level)
- g_fields = []
- for field in fields:
- if not grid.has_key(field): grid[field] = \
- np.zeros(grid.ActiveDimensions, dtype=self[field].dtype)
- g_fields.append(grid[field])
- c_fields = [self[field] for field in fields]
- FillRegion(ref_ratio,
- grid.get_global_startindex(), self.global_startindex,
- c_fields, g_fields,
- self.ActiveDimensions, grid.ActiveDimensions,
- grid.child_mask, self.domain_width, ll, 1)
-
- @property
- def LeftEdge(self):
- return self.left_edge
-
- @property
- def RightEdge(self):
- return self.right_edge
-
+ def get_data(self, fields = None):
+ fields = self._determine_fields(ensure_list(fields))
+ # We need a new tree for every single set of fields we add
+ if len(fields) == 0: return
+ # TODO:
+ # * Implement NeedsOriginalGrid catch
+ # * Generate some fields inside the data container
+ output_fields = [np.zeros(self.ActiveDimensions, dtype="float64")
+ for field in fields]
+ for chunk in self._data_source.chunks(fields, "io"):
+ input_fields = [chunk[field] for field in fields]
+ fill_region(input_fields, output_fields, self.level,
+ self.global_startindex, chunk.icoords, chunk.ires)
+ for name, v in zip(fields, output_fields):
+ self[name] = v
class YTSmoothedCoveringGridBase(YTCoveringGridBase):
_type_name = "smoothed_covering_grid"
@@ -1041,132 +576,3 @@
def flush_data(self, *args, **kwargs):
raise KeyError("Can't do this")
-
-class YTFixedResProjectionBase(YTSelectionContainer2D):
- _top_node = "/Projections"
- _type_name = "fixed_res_proj"
- _con_args = ('axis', 'field', 'weight_field')
- def __init__(self, axis, level, left_edge, dims, pf=None,
- field_parameters = None):
- """
- This is a data structure that projects grids, but only to fixed (rather
- than variable) resolution.
-
- This object is typically accessed through the `fixed_res_proj` object
- that hangs off of hierarchy objects. This projection mechanism is much
- simpler than the standard, variable-resolution projection. Rather than
- attempt to identify the highest-resolution element along every possible
- line of sight, this data structure simply deposits every cell into one
- of a fixed number of bins. It is suitable for inline analysis, and it
- should scale nicely.
-
- Parameters
- ----------
- axis : int
- The axis along which to project. Can be 0, 1, or 2 for x, y, z.
- level : int
- This is the level to which values will be projected. Note that
- the pixel size in the projection will be identical to a cell at
- this level of refinement in the simulation.
- left_edge : array of ints
- The left edge, in level-local integer coordinates, of the
- projection
- dims : array of ints
- The dimensions of the projection (which, in concert with the
- left_edge, serves to define its right edge.)
- field_parameters : dict of items
- Any additional values are passed as field parameters that can be
- accessed by generated fields.
-
- Examples
- --------
-
- >>> pf = load("RedshiftOutput0005")
- >>> fproj = pf.h.fixed_res_proj(1, [0, 0, 0], [64, 64, 64], ["Density"])
- >>> print fproj["Density"]
- """
- YTSelectionContainer2D.__init__(self, axis, pf, field_parameters)
- self.left_edge = np.array(left_edge)
- self.level = level
- self.dds = self.pf.h.select_grids(self.level)[0].dds.copy()
- self.dims = np.array([dims]*2)
- self.ActiveDimensions = np.array([dims]*3, dtype='int32')
- self.right_edge = self.left_edge + self.ActiveDimensions*self.dds
- self.global_startindex = np.rint((self.left_edge - self.pf.domain_left_edge)
- /self.dds).astype('int64')
- self._dls = {}
- self.domain_width = np.rint((self.pf.domain_right_edge -
- self.pf.domain_left_edge)/self.dds).astype('int64')
- self._refresh_data()
-
- def _get_list_of_grids(self):
- if self._grids is not None: return
- if np.any(self.left_edge < self.pf.domain_left_edge) or \
- np.any(self.right_edge > self.pf.domain_right_edge):
- grids,ind = self.pf.hierarchy.get_periodic_box_grids(
- self.left_edge, self.right_edge)
- else:
- grids,ind = self.pf.hierarchy.get_box_grids(
- self.left_edge, self.right_edge)
- level_ind = (self.pf.hierarchy.grid_levels.ravel()[ind] <= self.level)
- sort_ind = np.argsort(self.pf.h.grid_levels.ravel()[ind][level_ind])
- self._grids = self.pf.hierarchy.grids[ind][level_ind][(sort_ind,)][::-1]
-
- def _generate_coords(self):
- xax = x_dict[self.axis]
- yax = y_dict[self.axis]
- ci = self.left_edge + self.dds*0.5
- cf = self.left_edge + self.dds*(self.ActiveDimensions-0.5)
- cx = np.mgrid[ci[xax]:cf[xax]:self.ActiveDimensions[xax]*1j]
- cy = np.mgrid[ci[yax]:cf[yax]:self.ActiveDimensions[yax]*1j]
- blank = np.ones( (self.ActiveDimensions[xax],
- self.ActiveDimensions[yax]), dtype='float64')
- self['px'] = cx[None,:] * blank
- self['py'] = cx[:,None] * blank
- self['pdx'] = self.dds[xax]
- self['pdy'] = self.dds[yax]
-
- #@time_execution
- def get_data(self, fields):
- """
- Iterates over the list of fields and generates/reads them all.
- """
- self._get_list_of_grids()
- if not self.has_key('pdx'):
- self._generate_coords()
- fields_to_get = ensure_list(fields)
- if len(fields_to_get) == 0: return
- temp_data = {}
- for field in fields_to_get:
- self[field] = np.zeros(self.dims, dtype='float64')
- dls = self.__setup_dls(fields_to_get)
- for i,grid in enumerate(self._get_grids()):
- mylog.debug("Getting fields from %s", i)
- self._get_data_from_grid(grid, fields_to_get, dls)
- mylog.info("IO completed; summing")
- for field in fields_to_get:
- self[field] = self.comm.mpi_allreduce(self[field], op='sum')
- conv = self.pf.units[self.pf.field_info[field].projection_conversion]
- self[field] *= conv
-
- def __setup_dls(self, fields):
- dls = {}
- for level in range(self.level+1):
- dls[level] = []
- grid = self.select_grids(level)[0]
- for field in fields:
- if field is None: continue
- dls[level].append(float(just_one(grid['d%s' % axis_names[self.axis]])))
- return dls
-
- #@restore_grid_state
- def _get_data_from_grid(self, grid, fields, dls):
- g_fields = [grid[field].astype("float64") for field in fields]
- c_fields = [self[field] for field in fields]
- ref_ratio = self.pf.refine_by**(self.level - grid.Level)
- FillBuffer(ref_ratio,
- grid.get_global_startindex(), self.global_startindex,
- c_fields, g_fields,
- self.ActiveDimensions, grid.ActiveDimensions,
- grid.child_mask, self.domain_width, dls[grid.Level],
- self.axis)
diff -r 5a01ae4eb7201ff937ec80f8176944dfa8c8b33d -r 8f40060e1c99d1fa45a328278beb33fd2e7963b5 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -496,11 +496,11 @@
old_size, self.size = self.size, chunk.data_size
old_chunk, self._current_chunk = self._current_chunk, chunk
old_locked, self._locked = self._locked, False
- self.shape = (self.size,)
+ #self.shape = (self.size,)
yield
self.field_data = old_field_data
self.size = old_size
- self.shape = (old_size,)
+ #self.shape = (old_size,)
self._current_chunk = old_chunk
self._locked = old_locked
diff -r 5a01ae4eb7201ff937ec80f8176944dfa8c8b33d -r 8f40060e1c99d1fa45a328278beb33fd2e7963b5 yt/data_objects/tests/test_covering_grid.py
--- a/yt/data_objects/tests/test_covering_grid.py
+++ b/yt/data_objects/tests/test_covering_grid.py
@@ -8,7 +8,6 @@
def test_covering_grid():
# We decompose in different ways
- return
for level in [0, 1, 2]:
for nprocs in [1, 2, 4, 8]:
pf = fake_random_pf(16, nprocs = nprocs)
https://bitbucket.org/yt_analysis/yt-3.0/commits/4b45c2ba3de4/
changeset: 4b45c2ba3de4
branch: yt-3.0
user: MatthewTurk
date: 2013-01-09 12:41:50
summary: Refactored base SelectionDataContainer. Tests for covering grid now pass.
affected #: 2 files
diff -r 8f40060e1c99d1fa45a328278beb33fd2e7963b5 -r 4b45c2ba3de435b2194dedec0f22f42cbea0971c yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -415,6 +415,12 @@
# TODO:
# * Implement NeedsOriginalGrid catch
# * Generate some fields inside the data container
+ fields_to_get = self._identify_dependencies(fields)
+ fill, gen = self.pf.h._split_fields(fields_to_get)
+ if len(fill) > 0: self._fill_fields(fill)
+ if len(gen) > 0: self._generate_fields(gen)
+
+ def _fill_fields(self, fields):
output_fields = [np.zeros(self.ActiveDimensions, dtype="float64")
for field in fields]
for chunk in self._data_source.chunks(fields, "io"):
@@ -423,6 +429,7 @@
self.global_startindex, chunk.icoords, chunk.ires)
for name, v in zip(fields, output_fields):
self[name] = v
+
class YTSmoothedCoveringGridBase(YTCoveringGridBase):
_type_name = "smoothed_covering_grid"
diff -r 8f40060e1c99d1fa45a328278beb33fd2e7963b5 -r 4b45c2ba3de435b2194dedec0f22f42cbea0971c yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -429,6 +429,19 @@
# NOTE: we yield before releasing the context
yield self
+ def _identify_dependencies(self, fields_to_get):
+ inspected = 0
+ fields_to_get = fields_to_get[:]
+ for ftype, field in itertools.cycle(fields_to_get):
+ if inspected >= len(fields_to_get): break
+ inspected += 1
+ if field not in self.pf.field_dependencies: continue
+ fd = self.pf.field_dependencies[field]
+ requested = self._determine_fields(list(set(fd.requested)))
+ deps = [d for d in requested if d not in fields_to_get]
+ fields_to_get += deps
+ return fields_to_get
+
def get_data(self, fields=None):
if self._current_chunk is None:
self.hierarchy._identify_base_chunk(self)
@@ -441,15 +454,7 @@
elif self._locked == True:
raise GenerationInProgress(fields)
# At this point, we want to figure out *all* our dependencies.
- inspected = 0
- for ftype, field in itertools.cycle(fields_to_get):
- if inspected >= len(fields_to_get): break
- inspected += 1
- if field not in self.pf.field_dependencies: continue
- fd = self.pf.field_dependencies[field]
- requested = self._determine_fields(list(set(fd.requested)))
- deps = [d for d in requested if d not in fields_to_get]
- fields_to_get += deps
+ fields_to_get = self._identify_dependencies(fields_to_get)
# We now split up into readers for the types of fields
fluids, particles = [], []
for ftype, fname in fields_to_get:
@@ -469,6 +474,9 @@
particles, self, self._current_chunk)
self.field_data.update(read_particles)
fields_to_generate = gen_fluids + gen_particles
+ self._generate_fields(fields_to_generate)
+
+ def _generate_fields(self, fields_to_generate):
index = 0
with self._field_lock():
while any(f not in self.field_data for f in fields_to_generate):
https://bitbucket.org/yt_analysis/yt-3.0/commits/a2a284dadb7d/
changeset: a2a284dadb7d
branch: yt-3.0
user: MatthewTurk
date: 2013-01-09 13:16:56
summary: Supporting NeedsOriginalGrid. Adding test for same.
affected #: 2 files
diff -r 4b45c2ba3de435b2194dedec0f22f42cbea0971c -r a2a284dadb7d4343fc3ced4d64ec24750ac74cbf yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -416,10 +416,21 @@
# * Implement NeedsOriginalGrid catch
# * Generate some fields inside the data container
fields_to_get = self._identify_dependencies(fields)
- fill, gen = self.pf.h._split_fields(fields_to_get)
+ fill, gen = self._split_fields(fields_to_get)
if len(fill) > 0: self._fill_fields(fill)
if len(gen) > 0: self._generate_fields(gen)
+ def _split_fields(self, fields_to_get):
+ fill, gen = self.pf.h._split_fields(fields_to_get)
+ for field in gen:
+ finfo = self.pf._get_field_info(*field)
+ try:
+ finfo.check_available(self)
+ except NeedsOriginalGrid:
+ fill.append(field)
+ gen = [f for f in gen if f not in fill]
+ return fill, gen
+
def _fill_fields(self, fields):
output_fields = [np.zeros(self.ActiveDimensions, dtype="float64")
for field in fields]
@@ -470,7 +481,7 @@
# interpolation but are not directly inside our bounds
buffer = ((self.pf.domain_right_edge - self.pf.domain_left_edge)
/ self.pf.domain_dimensions).max()
- self._base_region = self.pf.geometry.region(
+ self._base_region = self.pf.h.region(
self.center,
self.left_edge - buffer,
self.right_edge + buffer)
diff -r 4b45c2ba3de435b2194dedec0f22f42cbea0971c -r a2a284dadb7d4343fc3ced4d64ec24750ac74cbf yt/data_objects/tests/test_covering_grid.py
--- a/yt/data_objects/tests/test_covering_grid.py
+++ b/yt/data_objects/tests/test_covering_grid.py
@@ -16,6 +16,7 @@
dn * pf.domain_dimensions)
yield assert_equal, cg["Ones"].max(), 1.0
yield assert_equal, cg["Ones"].min(), 1.0
+ yield assert_equal, cg["GridLevel"], 0
yield assert_equal, cg["CellVolume"].sum(), pf.domain_width.prod()
for g in pf.h.grids:
di = g.get_global_startindex()
https://bitbucket.org/yt_analysis/yt-3.0/commits/2d6776fc7122/
changeset: 2d6776fc7122
branch: yt-3.0
user: MatthewTurk
date: 2013-01-09 13:48:20
summary: Covering grid: Don't regenerate existing fields.
affected #: 2 files
diff -r a2a284dadb7d4343fc3ced4d64ec24750ac74cbf -r 2d6776fc71225d7dfd6096d82a39fdb03bf4ddc0 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -410,12 +410,9 @@
def get_data(self, fields = None):
fields = self._determine_fields(ensure_list(fields))
- # We need a new tree for every single set of fields we add
- if len(fields) == 0: return
- # TODO:
- # * Implement NeedsOriginalGrid catch
- # * Generate some fields inside the data container
- fields_to_get = self._identify_dependencies(fields)
+ fields_to_get = [f for f in fields if f not in self.field_data]
+ fields_to_get = self._identify_dependencies(fields_to_get)
+ if len(fields_to_get) == 0: return
fill, gen = self._split_fields(fields_to_get)
if len(fill) > 0: self._fill_fields(fill)
if len(gen) > 0: self._generate_fields(gen)
@@ -440,7 +437,6 @@
self.global_startindex, chunk.icoords, chunk.ires)
for name, v in zip(fields, output_fields):
self[name] = v
-
class YTSmoothedCoveringGridBase(YTCoveringGridBase):
_type_name = "smoothed_covering_grid"
diff -r a2a284dadb7d4343fc3ced4d64ec24750ac74cbf -r 2d6776fc71225d7dfd6096d82a39fdb03bf4ddc0 yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -451,6 +451,9 @@
nfield[ni, nj, nk] = 1
return nfield.astype("bool")
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
def fill_region(input_fields, output_fields,
np.int32_t output_level,
np.ndarray[np.int64_t, ndim=1] left_index,
https://bitbucket.org/yt_analysis/yt-3.0/commits/17d6b8d16dfb/
changeset: 17d6b8d16dfb
branch: yt-3.0
user: MatthewTurk
date: 2013-01-09 17:32:33
summary: Smoothed covering grid tests pass.
affected #: 2 files
diff -r 2d6776fc71225d7dfd6096d82a39fdb03bf4ddc0 -r 17d6b8d16dfbe8fea5a7bf0dba646e5afcb90d61 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -403,6 +403,9 @@
self.global_startindex = np.rint((self.left_edge-self.pf.domain_left_edge)/self.dds).astype('int64')
self.domain_width = np.rint((self.pf.domain_right_edge -
self.pf.domain_left_edge)/self.dds).astype('int64')
+ self._setup_data_source()
+
+ def _setup_data_source(self):
self._data_source = self.pf.h.region(
self.center, self.left_edge, self.right_edge)
self._data_source.min_level = 0
@@ -473,120 +476,77 @@
self.global_endindex = None
YTCoveringGridBase.__init__(self, *args, **kwargs)
self._final_start_index = self.global_startindex
+
+ def _setup_data_source(self, level = 0):
# We need a buffer region to allow for zones that contribute to the
# interpolation but are not directly inside our bounds
buffer = ((self.pf.domain_right_edge - self.pf.domain_left_edge)
/ self.pf.domain_dimensions).max()
- self._base_region = self.pf.h.region(
+ self._data_source = self.pf.h.region(
self.center,
self.left_edge - buffer,
self.right_edge + buffer)
+ self._data_source.min_level = 0
+ self._data_source.max_level = level
- def get_data(self, field):
- self._get_list_of_grids()
- # We don't generate coordinates here.
- if field == None:
- fields = self.fields[:]
- else:
- fields = ensure_list(field)
- fields_to_get = []
- for field in fields:
- if self.field_data.has_key(field): continue
- if field not in self.hierarchy.field_list:
- try:
- #print "Generating", field
- self._generate_field(field)
- continue
- except NeedsOriginalGrid, ngt_exception:
- pass
- fields_to_get.append(field)
- if len(fields_to_get) == 0: return
- # Note that, thanks to some trickery, we have different dimensions
- # on the field than one might think from looking at the dx and the
- # L/R edges.
- # We jump-start our task here
- mylog.debug("Getting fields %s from %s possible grids",
- fields_to_get, len(self._grids))
- self._update_level_state(0, fields_to_get)
- if self._use_pbar: pbar = \
- get_pbar('Searching grids for values ', len(self._grids))
- # The grids are assumed to be pre-sorted
- last_level = 0
- for gi, grid in enumerate(self._grids):
- if self._use_pbar: pbar.update(gi)
- if grid.Level > last_level and grid.Level <= self.level:
- mylog.debug("Updating level state to %s", last_level + 1)
- self._update_level_state(last_level + 1)
- self._refine(1, fields_to_get)
- last_level = grid.Level
- self._get_data_from_grid(grid, fields_to_get)
- while last_level < self.level:
- mylog.debug("Grid-free refinement %s to %s", last_level, last_level + 1)
- self._update_level_state(last_level + 1)
- self._refine(1, fields_to_get)
- last_level += 1
- if self.level > 0:
- for field in fields_to_get:
- self[field] = self[field][1:-1,1:-1,1:-1]
- if np.any(self[field] == -999):
- # and self.dx < self.hierarchy.grids[0].dx:
- n_bad = (self[field]==-999).sum()
- mylog.error("Covering problem: %s cells are uncovered", n_bad)
- raise KeyError(n_bad)
- if self._use_pbar: pbar.finish()
+ def _fill_fields(self, fields):
+ self._current_level = -1
+ output_fields = self._initialize_fields(fields)
+ for level in range(self.level + 1):
+ self._update_level_state(level)
+ self._setup_data_source(level)
+ for chunk in self._data_source.chunks(fields, "io"):
+ input_fields = [chunk[field] for field in fields]
+ fill_region(input_fields, output_fields, level,
+ self.global_startindex, chunk.icoords, chunk.ires)
+ if self.level > level:
+ output_fields = self._refine(1, output_fields)
+ for name, v in zip(fields, output_fields):
+ self[name] = v
def _update_level_state(self, level, fields = None):
+ if self._current_level == level: return
dx = self._base_dx / self.pf.refine_by**level
- self.field_data['cdx'] = dx[0]
- self.field_data['cdy'] = dx[1]
- self.field_data['cdz'] = dx[2]
+ self.cdx = dx
LL = self.left_edge - self.pf.domain_left_edge
self._old_global_startindex = self.global_startindex
self.global_startindex = np.rint(LL / dx).astype('int64') - 1
self.domain_width = np.rint((self.pf.domain_right_edge -
self.pf.domain_left_edge)/dx).astype('int64')
- if level == 0 and self.level > 0:
+ self._current_level = level
+
+ def _initialize_fields(self, fields):
+ field_return = []
+ self._update_level_state(0)
+ dx = self._base_dx
+ LL = self.left_edge - self.pf.domain_left_edge
+ if self.level > 0:
# We use one grid cell at LEAST, plus one buffer on all sides
idims = np.rint((self.right_edge-self.left_edge)/dx).astype('int64') + 2
- fields = ensure_list(fields)
for field in fields:
- self.field_data[field] = np.zeros(idims,dtype='float64')-999
+ field_return.append(np.zeros(idims,dtype='float64')-999)
self._cur_dims = idims.astype("int32")
- elif level == 0 and self.level == 0:
+ elif self.level == 0:
DLE = self.pf.domain_left_edge
- self.global_startindex = np.array(np.floor(LL/ dx), dtype='int64')
+ self.global_startindex = np.array(np.floor(LL/dx), dtype='int64')
idims = np.rint((self.ActiveDimensions*self.dds)/dx).astype('int64')
fields = ensure_list(fields)
for field in fields:
- self.field_data[field] = np.zeros(idims,dtype='float64')-999
+ field_return.append(np.zeros(idims,dtype='float64')-999)
self._cur_dims = idims.astype("int32")
+ return field_return
def _refine(self, dlevel, fields):
rf = float(self.pf.refine_by**dlevel)
input_left = (self._old_global_startindex + 0.5) * rf
- dx = np.fromiter((self['cd%s' % ax] for ax in 'xyz'), count=3, dtype='float64')
- output_dims = np.rint((self.ActiveDimensions*self.dds)/dx+0.5).astype('int32') + 2
+ output_dims = np.rint((self.ActiveDimensions*self.dds)/self.cdx+0.5).astype('int32') + 2
self._cur_dims = output_dims
-
- for field in fields:
+ return_fields = []
+ for input_field in fields:
output_field = np.zeros(output_dims, dtype="float64")
output_left = self.global_startindex + 0.5
- ghost_zone_interpolate(rf, self[field], input_left,
+ ghost_zone_interpolate(rf, input_field, input_left,
output_field, output_left)
- self.field_data[field] = output_field
-
- @restore_field_information_state
- def _get_data_from_grid(self, grid, fields):
- fields = ensure_list(fields)
- g_fields = [grid[field].astype("float64") for field in fields]
- c_fields = [self.field_data[field] for field in fields]
- count = FillRegion(1,
- grid.get_global_startindex(), self.global_startindex,
- c_fields, g_fields,
- self._cur_dims, grid.ActiveDimensions,
- grid.child_mask, self.domain_width, 1, 0)
- return count
-
- def flush_data(self, *args, **kwargs):
- raise KeyError("Can't do this")
+ return_fields.append(output_field)
+ return return_fields
diff -r 2d6776fc71225d7dfd6096d82a39fdb03bf4ddc0 -r 17d6b8d16dfbe8fea5a7bf0dba646e5afcb90d61 yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -461,7 +461,7 @@
np.ndarray[np.int64_t, ndim=1] ires,
np.int64_t refine_by = 2
):
- cdef int i, n, skip
+ cdef int i, n
cdef np.int64_t iind[3], oind[3], dim[3], oi, oj, ok, rf
cdef np.ndarray[np.float64_t, ndim=3] ofield
cdef np.ndarray[np.float64_t, ndim=1] ifield
https://bitbucket.org/yt_analysis/yt-3.0/commits/bcb74451defa/
changeset: bcb74451defa
branch: yt-3.0
user: MatthewTurk
date: 2013-01-09 17:43:47
summary: Help the other tests fail in a clearer way.
affected #: 3 files
diff -r 17d6b8d16dfbe8fea5a7bf0dba646e5afcb90d61 -r bcb74451defa2bbdb6ed0a344e6e76037328bd82 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -443,6 +443,7 @@
class YTSmoothedCoveringGridBase(YTCoveringGridBase):
_type_name = "smoothed_covering_grid"
+ filename = None
@wraps(YTCoveringGridBase.__init__)
def __init__(self, *args, **kwargs):
"""A 3D region with all data extracted and interpolated to a
diff -r 17d6b8d16dfbe8fea5a7bf0dba646e5afcb90d61 -r bcb74451defa2bbdb6ed0a344e6e76037328bd82 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -415,7 +415,8 @@
if self._selector is not None: return self._selector
sclass = getattr(yt.geometry.selection_routines,
"%s_selector" % self._type_name, None)
- if sclass is None: raise NotImplementedError
+ if sclass is None:
+ raise YTDataSelectorNotImplemented(self._type_name)
self._selector = sclass(self)
return self._selector
diff -r 17d6b8d16dfbe8fea5a7bf0dba646e5afcb90d61 -r bcb74451defa2bbdb6ed0a344e6e76037328bd82 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 YTDataSelectorNotImplemented(YTException):
+ def __init__(self, class_name):
+ self.class_name = class_name
+
+ def __str__(self):
+ return "Data selector '%s' not implemented." % (self.class_name)
https://bitbucket.org/yt_analysis/yt-3.0/commits/9d3ccd144c92/
changeset: 9d3ccd144c92
branch: yt-3.0
user: MatthewTurk
date: 2013-01-09 19:04:41
summary: Refinement mechanism for smoothed grids now well-behaved.
Grids themselves fail to correclty fill.
affected #: 2 files
diff -r bcb74451defa2bbdb6ed0a344e6e76037328bd82 -r 9d3ccd144c92e06c606cdb48459219f0821ad882 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -487,22 +487,24 @@
self.center,
self.left_edge - buffer,
self.right_edge + buffer)
- self._data_source.min_level = 0
+ self._data_source.min_level = level
self._data_source.max_level = level
def _fill_fields(self, fields):
self._current_level = -1
output_fields = self._initialize_fields(fields)
+ tot = 0
for level in range(self.level + 1):
- self._update_level_state(level)
self._setup_data_source(level)
for chunk in self._data_source.chunks(fields, "io"):
input_fields = [chunk[field] for field in fields]
- fill_region(input_fields, output_fields, level,
+ tot += fill_region(input_fields, output_fields, level,
self.global_startindex, chunk.icoords, chunk.ires)
if self.level > level:
+ self._update_level_state(level + 1)
output_fields = self._refine(1, output_fields)
for name, v in zip(fields, output_fields):
+ if self.level > 0: v = v[1:-1,1:-1,1:-1]
self[name] = v
def _update_level_state(self, level, fields = None):
diff -r bcb74451defa2bbdb6ed0a344e6e76037328bd82 -r 9d3ccd144c92e06c606cdb48459219f0821ad882 yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -462,6 +462,7 @@
np.int64_t refine_by = 2
):
cdef int i, n
+ cdef np.int64_t tot
cdef np.int64_t iind[3], oind[3], dim[3], oi, oj, ok, rf
cdef np.ndarray[np.float64_t, ndim=3] ofield
cdef np.ndarray[np.float64_t, ndim=1] ifield
@@ -469,6 +470,7 @@
for i in range(3):
dim[i] = output_fields[0].shape[i]
for n in range(nf):
+ tot = 0
ofield = output_fields[n]
ifield = input_fields[n]
for i in range(ipos.shape[0]):
@@ -489,3 +491,5 @@
continue
ofield[oind[0], oind[1], oind[2]] = \
ifield[i]
+ tot += 1
+ return tot
https://bitbucket.org/yt_analysis/yt-3.0/commits/9d78e820465f/
changeset: 9d78e820465f
branch: yt-3.0
user: MatthewTurk
date: 2013-01-09 20:25:26
summary: Initial creation of a LevelState object for refining smoothed covering grids.
affected #: 1 file
diff -r 9d3ccd144c92e06c606cdb48459219f0821ad882 -r 9d78e820465fce856dc7b3e70da9227a8dbfb2c1 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -441,6 +441,15 @@
for name, v in zip(fields, output_fields):
self[name] = v
+class LevelState(object):
+ current_dx = None
+ current_dims = None
+ current_level = None
+ global_startindex = None
+ old_global_startindex = None
+ domain_iwidth = None
+ fields = None
+
class YTSmoothedCoveringGridBase(YTCoveringGridBase):
_type_name = "smoothed_covering_grid"
filename = None
@@ -491,65 +500,66 @@
self._data_source.max_level = level
def _fill_fields(self, fields):
- self._current_level = -1
- output_fields = self._initialize_fields(fields)
+ ls = self._initialize_level_state(fields)
tot = 0
for level in range(self.level + 1):
- self._setup_data_source(level)
for chunk in self._data_source.chunks(fields, "io"):
input_fields = [chunk[field] for field in fields]
- tot += fill_region(input_fields, output_fields, level,
- self.global_startindex, chunk.icoords, chunk.ires)
- if self.level > level:
- self._update_level_state(level + 1)
- output_fields = self._refine(1, output_fields)
- for name, v in zip(fields, output_fields):
+ tot += fill_region(input_fields, ls.fields, level,
+ ls.global_startindex, chunk.icoords,
+ chunk.ires)
+ if tot == 0:
+ raise RuntimeError
+ self._update_level_state(ls)
+ for name, v in zip(fields, ls.fields):
if self.level > 0: v = v[1:-1,1:-1,1:-1]
self[name] = v
+ print name, tot, v.min(), v.max()
- def _update_level_state(self, level, fields = None):
- if self._current_level == level: return
- dx = self._base_dx / self.pf.refine_by**level
- self.cdx = dx
+ def _initialize_level_state(self, fields):
+ ls = LevelState()
+ ls.current_dx = self._base_dx
+ ls.current_level = 0
LL = self.left_edge - self.pf.domain_left_edge
- self._old_global_startindex = self.global_startindex
- self.global_startindex = np.rint(LL / dx).astype('int64') - 1
- self.domain_width = np.rint((self.pf.domain_right_edge -
- self.pf.domain_left_edge)/dx).astype('int64')
- self._current_level = level
-
- def _initialize_fields(self, fields):
- field_return = []
- self._update_level_state(0)
- dx = self._base_dx
- LL = self.left_edge - self.pf.domain_left_edge
+ ls.global_startindex = np.rint(LL / ls.current_dx).astype('int64') - 1
+ ls.domain_iwidth = np.rint((self.pf.domain_right_edge -
+ self.pf.domain_left_edge)/ls.current_dx).astype('int64')
if self.level > 0:
# We use one grid cell at LEAST, plus one buffer on all sides
- idims = np.rint((self.right_edge-self.left_edge)/dx).astype('int64') + 2
- for field in fields:
- field_return.append(np.zeros(idims,dtype='float64')-999)
- self._cur_dims = idims.astype("int32")
+ width = self.right_edge-self.left_edge
+ idims = np.rint(width/ls.current_dx).astype('int64') + 2
elif self.level == 0:
- DLE = self.pf.domain_left_edge
- self.global_startindex = np.array(np.floor(LL/dx), dtype='int64')
- idims = np.rint((self.ActiveDimensions*self.dds)/dx).astype('int64')
- fields = ensure_list(fields)
- for field in fields:
- field_return.append(np.zeros(idims,dtype='float64')-999)
- self._cur_dims = idims.astype("int32")
- return field_return
+ ls.global_startindex = np.array(np.floor(LL/ls.current_dx), dtype='int64')
+ idims = np.rint((self.ActiveDimensions*self.dds)/ls.current_dx).astype('int64')
+ ls.current_dims = idims.astype("int32")
+ ls.fields = [np.zeros(idims, dtype="float64")-999 for field in fields]
+ self._setup_data_source(0)
+ return ls
- def _refine(self, dlevel, fields):
- rf = float(self.pf.refine_by**dlevel)
+ def _update_level_state(self, level_state):
+ ls = level_state
+ if ls.current_level >= self.level: return
+ ls.current_level += 1
+ self._setup_data_source(ls.current_level)
+ ls.current_dx = self._base_dx / self.pf.refine_by**ls.current_level
+ LL = self.left_edge - self.pf.domain_left_edge
+ ls.old_global_startindex = ls.global_startindex
+ ls.global_startindex = np.rint(LL / ls.current_dx).astype('int64') - 1
+ ls.domain_iwidth = np.rint(self.pf.domain_width/ls.current_dx).astype('int64')
+ self._refine(ls)
- input_left = (self._old_global_startindex + 0.5) * rf
- output_dims = np.rint((self.ActiveDimensions*self.dds)/self.cdx+0.5).astype('int32') + 2
- self._cur_dims = output_dims
- return_fields = []
- for input_field in fields:
+ def _refine(self, level_state):
+ rf = float(self.pf.refine_by)
+ input_left = (level_state.old_global_startindex + 0.5) * rf
+ width = (self.ActiveDimensions*self.dds)
+ output_dims = np.rint(width/level_state.current_dx+0.5).astype("int32") + 2
+ level_state.current_dims = output_dims
+ new_fields = []
+ for input_field in level_state.fields:
output_field = np.zeros(output_dims, dtype="float64")
output_left = self.global_startindex + 0.5
ghost_zone_interpolate(rf, input_field, input_left,
output_field, output_left)
- return_fields.append(output_field)
- return return_fields
+ new_fields.append(output_field)
+ level_state.fields = new_fields
+
https://bitbucket.org/yt_analysis/yt-3.0/commits/6fc7073f4cec/
changeset: 6fc7073f4cec
branch: yt-3.0
user: MatthewTurk
date: 2013-01-09 23:09:30
summary: Fixing issues in RegionSelector and fill_region.
This fixes smoothed covering grids, which now provisionally work. Answers for
tested datasets are identical to previous answers. Early termination in
RegionSelector objects has been disabled.
affected #: 3 files
diff -r 9d78e820465fce856dc7b3e70da9227a8dbfb2c1 -r 6fc7073f4cec696f52a8b6d049bb8c6a2ffac4b8 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -449,6 +449,7 @@
old_global_startindex = None
domain_iwidth = None
fields = None
+ data_source = None
class YTSmoothedCoveringGridBase(YTCoveringGridBase):
_type_name = "smoothed_covering_grid"
@@ -487,34 +488,31 @@
YTCoveringGridBase.__init__(self, *args, **kwargs)
self._final_start_index = self.global_startindex
- def _setup_data_source(self, level = 0):
+ def _setup_data_source(self, level_state = None):
+ if level_state is None: return
# We need a buffer region to allow for zones that contribute to the
# interpolation but are not directly inside our bounds
- buffer = ((self.pf.domain_right_edge - self.pf.domain_left_edge)
- / self.pf.domain_dimensions).max()
- self._data_source = self.pf.h.region(
+ level_state.data_source = self.pf.h.region(
self.center,
- self.left_edge - buffer,
- self.right_edge + buffer)
- self._data_source.min_level = level
- self._data_source.max_level = level
+ self.left_edge - level_state.current_dx,
+ self.right_edge + level_state.current_dx)
+ level_state.data_source.min_level = level_state.current_level
+ level_state.data_source.max_level = level_state.current_level
def _fill_fields(self, fields):
ls = self._initialize_level_state(fields)
- tot = 0
for level in range(self.level + 1):
- for chunk in self._data_source.chunks(fields, "io"):
+ tot = 0
+ for chunk in ls.data_source.chunks(fields, "io"):
+ chunk[fields[0]]
input_fields = [chunk[field] for field in fields]
- tot += fill_region(input_fields, ls.fields, level,
+ tot += fill_region(input_fields, ls.fields, ls.current_level,
ls.global_startindex, chunk.icoords,
chunk.ires)
- if tot == 0:
- raise RuntimeError
self._update_level_state(ls)
for name, v in zip(fields, ls.fields):
if self.level > 0: v = v[1:-1,1:-1,1:-1]
self[name] = v
- print name, tot, v.min(), v.max()
def _initialize_level_state(self, fields):
ls = LevelState()
@@ -533,22 +531,19 @@
idims = np.rint((self.ActiveDimensions*self.dds)/ls.current_dx).astype('int64')
ls.current_dims = idims.astype("int32")
ls.fields = [np.zeros(idims, dtype="float64")-999 for field in fields]
- self._setup_data_source(0)
+ self._setup_data_source(ls)
return ls
def _update_level_state(self, level_state):
ls = level_state
if ls.current_level >= self.level: return
ls.current_level += 1
- self._setup_data_source(ls.current_level)
ls.current_dx = self._base_dx / self.pf.refine_by**ls.current_level
+ self._setup_data_source(ls)
LL = self.left_edge - self.pf.domain_left_edge
ls.old_global_startindex = ls.global_startindex
ls.global_startindex = np.rint(LL / ls.current_dx).astype('int64') - 1
ls.domain_iwidth = np.rint(self.pf.domain_width/ls.current_dx).astype('int64')
- self._refine(ls)
-
- def _refine(self, level_state):
rf = float(self.pf.refine_by)
input_left = (level_state.old_global_startindex + 0.5) * rf
width = (self.ActiveDimensions*self.dds)
@@ -562,4 +557,3 @@
output_field, output_left)
new_fields.append(output_field)
level_state.fields = new_fields
-
diff -r 9d78e820465fce856dc7b3e70da9227a8dbfb2c1 -r 6fc7073f4cec696f52a8b6d049bb8c6a2ffac4b8 yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -496,18 +496,24 @@
np.float64_t dds[3], int ind[3][2], int *check):
cdef int temp, i, all_inside
# Left pos is left_edge + 0.5 * dds
- all_inside = 1
+ # This means the grid is fully within the region
+ # Note vice versa!
+ all_inside = 1
for i in range(3):
+ # If the region starts after the grid does...
if self.left_edge[i] > left_edge[i]:
temp = <int> ((self.left_edge[i] - left_edge[i])/dds[i]) - 1
- ind[i][0] = iclip(temp, ind[i][0], ind[i][1] - 1)
+ #ind[i][0] = iclip(temp, ind[i][0], ind[i][1])
all_inside = 0
+ # If the region ends before the grid does...
if self.right_edge[i] < right_edge[i]:
- temp = <int> ((self.right_edge[i] - right_edge[i])/dds[i]) + 1
- ind[i][1] = iclip(temp, ind[i][0] + 1, ind[i][1])
+ temp = <int> ((self.right_edge[i] - left_edge[i])/dds[i]) + 1
+ #ind[i][1] = iclip(temp, ind[i][0], ind[i][1])
all_inside = 0
- check[0] = all_inside
-
+ if all_inside == 0:
+ check[0] = 1
+ else:
+ check[0] = 0
region_selector = RegionSelector
diff -r 9d78e820465fce856dc7b3e70da9227a8dbfb2c1 -r 6fc7073f4cec696f52a8b6d049bb8c6a2ffac4b8 yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -474,20 +474,20 @@
ofield = output_fields[n]
ifield = input_fields[n]
for i in range(ipos.shape[0]):
+ rf = refine_by**(output_level - ires[i])
for n in range(3):
- iind[n] = ipos[i, n] - left_index[n]
- rf = refine_by**(output_level - ires[i])
+ iind[n] = ipos[i, n] * rf - left_index[n]
for oi in range(rf):
- oind[0] = oi + iind[0] * rf
- if oind[0] < left_index[0] or oind[0] >= left_index[0] + dim[0]:
+ oind[0] = oi + iind[0]
+ if oind[0] < 0 or oind[0] >= dim[0]:
continue
for oj in range(rf):
- oind[1] = oj + iind[1] * rf
- if oind[1] < left_index[1] or oind[1] >= left_index[1] + dim[1]:
+ oind[1] = oj + iind[1]
+ if oind[1] < 0 or oind[1] >= dim[1]:
continue
for ok in range(rf):
- oind[2] = ok + iind[2] * rf
- if oind[2] < left_index[2] or oind[2] >= left_index[2] + dim[2]:
+ oind[2] = ok + iind[2]
+ if oind[2] < 0 or oind[2] >= dim[2]:
continue
ofield[oind[0], oind[1], oind[2]] = \
ifield[i]
https://bitbucket.org/yt_analysis/yt-3.0/commits/f2e039a1f6d4/
changeset: f2e039a1f6d4
branch: yt-3.0
user: MatthewTurk
date: 2013-01-10 01:49:02
summary: Fixing ValidateSpatial fields.
* Smoothed covering grids now (if appropriate) have links to _base_grid.
* Spatial chunking works slightly differently for ngz > 0.
* Fixed field_parameters for covering and smoothed covering grids.
* Re-enabled a bunch of tests; VorticityStretchingY_4 still fails for some reason.
affected #: 4 files
diff -r 6fc7073f4cec696f52a8b6d049bb8c6a2ffac4b8 -r f2e039a1f6d40342ce922b6544ec6a0fac891c83 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -368,8 +368,10 @@
_spatial = True
_type_name = "covering_grid"
_con_args = ('level', 'left_edge', 'ActiveDimensions')
+ _base_grid = None
def __init__(self, level, left_edge, dims, fields = None,
- pf = None, num_ghost_zones = 0, use_pbar = True, **kwargs):
+ pf = None, num_ghost_zones = 0, use_pbar = True,
+ field_parameters = None):
"""A 3D region with all data extracted to a single, specified
resolution.
@@ -390,8 +392,9 @@
dims=[128, 128, 128])
"""
- YTSelectionContainer3D.__init__(self, center=kwargs.pop("center", None),
- pf=pf, **kwargs)
+ YTSelectionContainer3D.__init__(self,
+ field_parameters.get("center", None),
+ pf, field_parameters)
self.left_edge = np.array(left_edge)
self.level = level
rdx = self.pf.domain_dimensions*self.pf.refine_by**level
diff -r 6fc7073f4cec696f52a8b6d049bb8c6a2ffac4b8 -r f2e039a1f6d40342ce922b6544ec6a0fac891c83 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -234,21 +234,33 @@
try:
self.pf.field_info[fname].check_available(gen_obj)
except NeedsGridType as ngt_exception:
- rv = np.empty(self.size, dtype="float64")
- ind = 0
- ngz = ngt_exception.ghost_zones
+ rv = self._generate_spatial_fluid(field, ngt_exception.ghost_zones)
+ else:
+ rv = self.pf.field_info[fname](gen_obj)
+ return rv
+
+ def _generate_spatial_fluid(self, field, ngz):
+ rv = np.empty(self.size, dtype="float64")
+ ind = 0
+ if ngz == 0:
for io_chunk in self.chunks([], "io"):
- for i,chunk in enumerate(self.chunks(field, "spatial", ngz = ngz)):
+ for i,chunk in enumerate(self.chunks(field, "spatial", ngz = 0)):
mask = self._current_chunk.objs[0].select(self.selector)
if mask is None: continue
- data = self[field]
- if ngz > 0:
- data = data[ngz:-ngz, ngz:-ngz, ngz:-ngz]
- data = data[mask]
+ data = self[field][mask]
rv[ind:ind+data.size] = data
ind += data.size
else:
- rv = self.pf.field_info[fname](gen_obj)
+ chunks = self.hierarchy._chunk(self, "spatial", ngz = ngz)
+ for i, chunk in enumerate(chunks):
+ with self._chunked_read(chunk):
+ gz = self._current_chunk.objs[0]
+ wogz = gz._base_grid
+ mask = wogz.select(self.selector)
+ if mask is None: continue
+ data = gz[field][ngz:-ngz, ngz:-ngz, ngz:-ngz][mask]
+ rv[ind:ind+data.size] = data
+ ind += data.size
return rv
def _generate_particle_field(self, field):
diff -r 6fc7073f4cec696f52a8b6d049bb8c6a2ffac4b8 -r f2e039a1f6d40342ce922b6544ec6a0fac891c83 yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -368,13 +368,18 @@
'use_pbar':False, 'fields':fields}
# This should update the arguments to set the field parameters to be
# those of this grid.
- kwargs.update(self.field_parameters)
+ field_parameters = {}
+ field_parameters.update(self.field_parameters)
if smoothed:
cube = self.hierarchy.smoothed_covering_grid(
- level, new_left_edge, **kwargs)
+ level, new_left_edge,
+ field_parameters = field_parameters,
+ **kwargs)
else:
- cube = self.hierarchy.covering_grid(level, new_left_edge, **kwargs)
-
+ cube = self.hierarchy.covering_grid(level, new_left_edge,
+ field_parameters = field_parameters,
+ **kwargs)
+ cube._base_grid = self
return cube
def get_vertex_centered_data(self, field, smoothed=True, no_ghost=False):
diff -r 6fc7073f4cec696f52a8b6d049bb8c6a2ffac4b8 -r f2e039a1f6d40342ce922b6544ec6a0fac891c83 yt/data_objects/tests/test_fields.py
--- a/yt/data_objects/tests/test_fields.py
+++ b/yt/data_objects/tests/test_fields.py
@@ -86,10 +86,6 @@
if field.startswith("particle"): continue
if field.startswith("CIC"): continue
if field.startswith("WeakLensingConvergence"): continue
- skip = False
- for v in FieldInfo[field].validators:
- if getattr(v, "ghost_zones", 0) > 0: skip = True
- if skip: continue
if FieldInfo[field].particle_type: continue
for nproc in [1, 4, 8]:
yield TestFieldAccess(field, nproc)
https://bitbucket.org/yt_analysis/yt-3.0/commits/736226b9a242/
changeset: 736226b9a242
branch: yt-3.0
user: MatthewTurk
date: 2013-01-10 02:20:31
summary: Fixed up field_parameters for covering grids.
Only remaining test failures are for not-implemented objects.
affected #: 1 file
diff -r f2e039a1f6d40342ce922b6544ec6a0fac891c83 -r 736226b9a242de0a34b14085f3d65c4f5db46769 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -392,9 +392,12 @@
dims=[128, 128, 128])
"""
+ if field_parameters is None:
+ center = None
+ else:
+ center = field_parameters.get("center", None)
YTSelectionContainer3D.__init__(self,
- field_parameters.get("center", None),
- pf, field_parameters)
+ center, pf, field_parameters)
self.left_edge = np.array(left_edge)
self.level = level
rdx = self.pf.domain_dimensions*self.pf.refine_by**level
https://bitbucket.org/yt_analysis/yt-3.0/commits/789bd796e365/
changeset: 789bd796e365
branch: yt-3.0
user: MatthewTurk
date: 2013-01-10 16:07:31
summary: Merging
affected #: 34 files
diff -r 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 yt/utilities/exceptions.py
--- a/yt/utilities/exceptions.py
+++ b/yt/utilities/exceptions.py
@@ -209,6 +209,13 @@
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
+
class YTDataSelectorNotImplemented(YTException):
def __init__(self, class_name):
self.class_name = class_name
diff -r 736226b9a242de0a34b14085f3d65c4f5db46769 -r 789bd796e365b1d277569ac4b0c3c5e0630f4504 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