[yt-svn] commit/yt: 13 new changesets
Bitbucket
commits-noreply at bitbucket.org
Tue Nov 20 12:56:25 PST 2012
13 new commits in yt:
https://bitbucket.org/yt_analysis/yt/changeset/d00afdd643c6/
changeset: d00afdd643c6
branch: yt
user: jzuhone
date: 2012-11-17 15:52:33
summary: Depth-first search for leaf grids given a set of positions.
affected #: 3 files
diff -r 928bf03996bb8a0704b89e5bede9ace9a57aafa2 -r d00afdd643c65062f357dbbe31fbfbd44c8b2818 yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -242,6 +242,23 @@
cond = np.logical_and(cond, self.LeftEdge[y] <= RE[:,y])
return cond
+ def is_in_grid(self, x, y, z) :
+ """
+ Generate a mask that shows which points in *x*, *y*, and *z*
+ fall within this grid's boundaries.
+ """
+ xcond = np.logical_and(x >= self.LeftEdge[0],
+ x < self.RightEdge[0])
+ ycond = np.logical_and(y >= self.LeftEdge[1],
+ y < self.RightEdge[1])
+ zcond = np.logical_and(z >= self.LeftEdge[2],
+ z < self.RightEdge[2])
+
+ cond = np.logical_and(xcond, ycond)
+ cond = np.logical_and(zcond, cond)
+
+ return cond
+
def __repr__(self):
return "AMRGridPatch_%04i" % (self.id)
diff -r 928bf03996bb8a0704b89e5bede9ace9a57aafa2 -r d00afdd643c65062f357dbbe31fbfbd44c8b2818 yt/data_objects/object_finding_mixin.py
--- a/yt/data_objects/object_finding_mixin.py
+++ b/yt/data_objects/object_finding_mixin.py
@@ -29,8 +29,10 @@
from yt.utilities.lib import \
get_box_grids_level, \
get_box_grids_below_level
+from yt.utilities.particle_utils import \
+ MatchPointsToGrids
-class ObjectFindingMixin(object):
+class ObjectFindingMixin(object) :
def find_ray_grids(self, coord, axis):
"""
@@ -110,6 +112,14 @@
ind = np.where(mask == 1)
return self.grids[ind], ind
+ def find_points(self, x, y, z) :
+ """
+ Returns the (objects, indices) of leaf grids containing a number of (x,y,z) points
+ """
+ find_pts = MatchPointsToGrids(self,x,y,z)
+ ind = find_pts()
+ return self.grids[ind], ind
+
def find_field_value_at_point(self, fields, coord):
r"""Find the value of fields at a point.
diff -r 928bf03996bb8a0704b89e5bede9ace9a57aafa2 -r d00afdd643c65062f357dbbe31fbfbd44c8b2818 yt/utilities/particle_utils.py
--- /dev/null
+++ b/yt/utilities/particle_utils.py
@@ -0,0 +1,95 @@
+"""
+Utilities for particles.
+
+Author: John ZuHone <jzuhone at gmail.com>
+Affiliation: NASA/Goddard Space Flight Center
+Homepage: http://yt-project.org/
+License:
+Copyright (C) 2012 John ZuHone. 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 *
+
+class MatchPointsToGrids(object) :
+
+ def __init__(self, h, x, y, z) :
+
+ self.num_points = len(x)
+
+ self.x = x
+ self.y = y
+ self.z = z
+
+ self.point_grids = -np.ones((self.num_points)).astype("int32")
+ self.all_idxs = np.arange((self.num_points))
+
+ self.root_grids = h.select_grids(0)
+
+ self.counter = 0
+
+ def __call__(self) :
+ """
+ Loops over every root grid and follows down each's subtree
+ until it finds leaf grids, then assigns that grid id to the points
+ in that grid
+ """
+
+ self.pbar = get_pbar("Finding Grids", self.num_points)
+
+ for root_grid in self.root_grids :
+
+ self.check_positions(root_grid)
+
+ self.pbar.finish()
+
+ return self.point_grids
+
+ def check_positions(self, grid) :
+ """
+ Recursive function to traverse up and down the hierarchy
+ """
+
+ if len(self.all_idxs) == 0 : return # We've run out of points
+
+ if len(grid.Children) > 0 :
+
+ # This grid has children
+
+ for child in grid.Children :
+
+ self.check_positions(child)
+
+ else :
+
+ # We hit a leaf grid
+
+ local_idxs = grid.is_in_grid(self.x[self.all_idxs], self.y[self.all_idxs],
+ self.z[self.all_idxs])
+
+ if np.any(local_idxs) :
+
+ self.point_grids[self.all_idxs[local_idxs]] = grid.id-grid._id_offset
+ self.all_idxs = self.all_idxs[np.logical_not(local_idxs)]
+
+ self.counter += local_idxs.sum()
+
+ self.pbar.update(self.counter)
+
+ return
+
https://bitbucket.org/yt_analysis/yt/changeset/d2343b2cc852/
changeset: d2343b2cc852
branch: yt
user: jzuhone
date: 2012-11-20 03:49:28
summary: Adding Cythonized GridTree for assigning points to grids
affected #: 3 files
diff -r d00afdd643c65062f357dbbe31fbfbd44c8b2818 -r d2343b2cc85275300a1f5ee29b06a0a861280f6b yt/utilities/lib/GridTree.pyx
--- /dev/null
+++ b/yt/utilities/lib/GridTree.pyx
@@ -0,0 +1,278 @@
+"""
+Matching points on the grid to specific grids
+
+Author: John ZuHone <jzuhone at gmail.com>
+Affiliation: NASA/Goddard Space Flight Center
+Homepage: http://yt-project.org/
+License:
+Copyright (C) 2012 John ZuHone. 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
+cimport numpy as np
+cimport cython
+
+from libc.stdlib cimport malloc, free
+
+cdef struct GridTreeNode :
+ int num_children
+ int level
+ int index
+ np.float64_t left_edge[3]
+ np.float64_t right_edge[3]
+ GridTreeNode **children
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef GridTreeNode Grid_initialize(np.ndarray[np.float64_t, ndim=1] le,
+ np.ndarray[np.float64_t, ndim=1] re,
+ int num_children, int level, int index) :
+
+ cdef GridTreeNode node
+ cdef int i
+
+ node.index = index
+ node.level = level
+ for i in range(3) :
+ node.left_edge[i] = le[i]
+ node.right_edge[i] = re[i]
+ node.num_children = num_children
+
+ if num_children > 0:
+ node.children = <GridTreeNode **> malloc(sizeof(GridTreeNode *) *
+ num_children)
+ for i in range(num_children) :
+ node.children[i] = NULL
+ else :
+ node.children = NULL
+
+ return node
+
+cdef class GridTree :
+
+ cdef GridTreeNode *grids
+ cdef GridTreeNode *root_grids
+ cdef int num_grids
+ cdef int num_root_grids
+ cdef int num_leaf_grids
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def __cinit__(self, int num_grids,
+ np.ndarray[np.float64_t, ndim=2] left_edge,
+ np.ndarray[np.float64_t, ndim=2] right_edge,
+ np.ndarray[np.int64_t, ndim=1] parent_ind,
+ np.ndarray[np.int64_t, ndim=1] level,
+ np.ndarray[np.int64_t, ndim=1] num_children) :
+
+ cdef int i, j, k
+ cdef np.ndarray[np.int64_t, ndim=1] child_ptr
+
+ child_ptr = np.zeros(num_grids, dtype='int64')
+
+ self.num_grids = num_grids
+ self.num_root_grids = 0
+ self.num_leaf_grids = 0
+
+ self.grids = <GridTreeNode *> malloc(sizeof(GridTreeNode) *
+ num_grids)
+
+ for i in range(num_grids) :
+
+ self.grids[i] = Grid_initialize(left_edge[i,:],
+ right_edge[i,:],
+ num_children[i],
+ level[i], i)
+ if level[i] == 0 :
+ self.num_root_grids += 1
+
+ if num_children[i] == 0 : self.num_leaf_grids += 1
+
+ self.root_grids = <GridTreeNode *> malloc(sizeof(GridTreeNode) *
+ self.num_root_grids)
+
+ k = 0
+
+ for i in range(num_grids) :
+
+ j = parent_ind[i]
+
+ if j >= 0:
+
+ self.grids[j].children[child_ptr[j]] = &self.grids[i]
+
+ child_ptr[j] += 1
+
+ else :
+
+ self.root_grids[k] = self.grids[i]
+
+ k = k + 1
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def return_tree_info(self) :
+
+ cdef int i, j
+
+ levels = []
+ indices = []
+ nchild = []
+ children = []
+
+ for i in range(self.num_grids) :
+
+ childs = []
+
+ levels.append(self.grids[i].level)
+ indices.append(self.grids[i].index)
+ nchild.append(self.grids[i].num_children)
+ for j in range(self.grids[i].num_children) :
+ childs.append(self.grids[i].children[j].index)
+ children.append(childs)
+
+ return indices, levels, nchild, children
+
+cdef class MatchPointsToGrids :
+
+ cdef int num_points
+ cdef int num_idxs_left
+ cdef int num_grids_walked
+ cdef np.float64_t * xp
+ cdef np.float64_t * yp
+ cdef np.float64_t * zp
+ cdef GridTree tree
+ cdef np.int64_t * all_idxs
+ cdef np.int64_t * point_grids
+ cdef np.int64_t * in_grid
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def __cinit__(self, GridTree tree,
+ int num_points,
+ np.ndarray[np.float64_t, ndim=1] x,
+ np.ndarray[np.float64_t, ndim=1] y,
+ np.ndarray[np.float64_t, ndim=1] z) :
+
+ cdef int i
+
+ self.num_points = num_points
+ self.num_idxs_left = num_points
+
+ self.xp = <np.float64_t *> malloc(sizeof(np.float64_t) *
+ num_points)
+ self.yp = <np.float64_t *> malloc(sizeof(np.float64_t) *
+ num_points)
+ self.zp = <np.float64_t *> malloc(sizeof(np.float64_t) *
+ num_points)
+ self.all_idxs = <np.int64_t *> malloc(sizeof(np.int64_t) *
+ num_points)
+ self.in_grid = <np.int64_t *> malloc(sizeof(np.int64_t) *
+ num_points)
+ self.point_grids = <np.int64_t *> malloc(sizeof(np.int64_t) *
+ num_points)
+
+ for i in range(num_points) :
+ self.xp[i] = x[i]
+ self.yp[i] = y[i]
+ self.zp[i] = z[i]
+ self.all_idxs[i] = i
+ self.in_grid[i] = 0
+ self.point_grids[i] = -1
+
+ self.tree = tree
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ def find_points_in_tree(self) :
+
+ cdef np.ndarray[np.int64_t, ndim=1] pt_grids
+ cdef int i
+
+ pt_grids = np.zeros(self.num_points, dtype='int64')
+
+ for i in range(self.tree.num_root_grids) :
+
+ self.check_positions(&self.tree.root_grids[i])
+
+ for i in range(self.num_points) :
+ pt_grids[i] = self.point_grids[i]
+
+ return pt_grids
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ cdef void check_positions(self, GridTreeNode * grid) :
+
+ cdef int i
+
+ if self.num_idxs_left == 0 : return
+
+ if grid.num_children > 0 :
+
+ for i in range(grid.num_children) :
+
+ self.check_positions(grid.children[i])
+
+ else :
+
+ self.is_in_grid(grid)
+
+ for i in range(self.num_points) :
+
+ if self.in_grid[i] == 1 :
+
+ #print self.num_idxs_left
+
+ self.point_grids[i] = grid.index
+ self.all_idxs[i] = -1
+ self.num_idxs_left -= 1
+
+ return
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef void is_in_grid(self, GridTreeNode * grid) :
+
+ cdef int i
+ cdef np.uint8_t xcond, ycond, zcond, cond
+
+ for i in range(self.num_points) :
+
+ if self.all_idxs[i] >= 0 :
+
+ xcond = self.xp[i] >= grid.left_edge[0] and self.xp[i] < grid.right_edge[0]
+ ycond = self.yp[i] >= grid.left_edge[1] and self.yp[i] < grid.right_edge[1]
+ zcond = self.zp[i] >= grid.left_edge[2] and self.zp[i] < grid.right_edge[2]
+
+ cond = xcond and ycond
+ cond = cond and zcond
+
+ if cond :
+ self.in_grid[i] = 1
+ else :
+ self.in_grid[i] = 0
+
+ else :
+
+ self.in_grid[i] = 0
diff -r d00afdd643c65062f357dbbe31fbfbd44c8b2818 -r d2343b2cc85275300a1f5ee29b06a0a861280f6b yt/utilities/lib/__init__.py
--- a/yt/utilities/lib/__init__.py
+++ b/yt/utilities/lib/__init__.py
@@ -38,3 +38,4 @@
from .RayIntegrators import *
from .grid_traversal import *
from .marching_cubes import *
+from .GridTree import *
diff -r d00afdd643c65062f357dbbe31fbfbd44c8b2818 -r d2343b2cc85275300a1f5ee29b06a0a861280f6b yt/utilities/lib/setup.py
--- a/yt/utilities/lib/setup.py
+++ b/yt/utilities/lib/setup.py
@@ -204,6 +204,10 @@
"yt/utilities/lib/field_interpolation_tables.pxd",
]
)
+ config.add_extension("GridTree",
+ ["yt/utilities/lib/GridTree.pyx"],
+ libraries=["m"], depends=["yt/utilities/lib/fp_utils.pxd"])
+
if os.environ.get("GPERFTOOLS", "no").upper() != "NO":
gpd = os.environ["GPERFTOOLS"]
idir = os.path.join(gpd, "include")
https://bitbucket.org/yt_analysis/yt/changeset/3283e02e6f55/
changeset: 3283e02e6f55
branch: yt
user: jzuhone
date: 2012-11-20 05:30:05
summary: Update to the GridTree. We should check each particle one-by-one and go down the the tree instead of checking each grid against lists of particles.
affected #: 1 file
diff -r d2343b2cc85275300a1f5ee29b06a0a861280f6b -r 3283e02e6f557a33094cb58cbc3246b0dca75b59 yt/utilities/lib/GridTree.pyx
--- a/yt/utilities/lib/GridTree.pyx
+++ b/yt/utilities/lib/GridTree.pyx
@@ -154,15 +154,11 @@
cdef class MatchPointsToGrids :
cdef int num_points
- cdef int num_idxs_left
- cdef int num_grids_walked
cdef np.float64_t * xp
cdef np.float64_t * yp
cdef np.float64_t * zp
cdef GridTree tree
- cdef np.int64_t * all_idxs
cdef np.int64_t * point_grids
- cdef np.int64_t * in_grid
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -176,7 +172,6 @@
cdef int i
self.num_points = num_points
- self.num_idxs_left = num_points
self.xp = <np.float64_t *> malloc(sizeof(np.float64_t) *
num_points)
@@ -184,10 +179,6 @@
num_points)
self.zp = <np.float64_t *> malloc(sizeof(np.float64_t) *
num_points)
- self.all_idxs = <np.int64_t *> malloc(sizeof(np.int64_t) *
- num_points)
- self.in_grid = <np.int64_t *> malloc(sizeof(np.int64_t) *
- num_points)
self.point_grids = <np.int64_t *> malloc(sizeof(np.int64_t) *
num_points)
@@ -195,8 +186,6 @@
self.xp[i] = x[i]
self.yp[i] = y[i]
self.zp[i] = z[i]
- self.all_idxs[i] = i
- self.in_grid[i] = 0
self.point_grids[i] = -1
self.tree = tree
@@ -206,13 +195,16 @@
def find_points_in_tree(self) :
cdef np.ndarray[np.int64_t, ndim=1] pt_grids
- cdef int i
+ cdef int i, j
pt_grids = np.zeros(self.num_points, dtype='int64')
-
- for i in range(self.tree.num_root_grids) :
- self.check_positions(&self.tree.root_grids[i])
+ for i in range(self.num_points) :
+
+ for j in range(self.tree.num_root_grids) :
+
+ self.check_position(i, self.xp[i], self.yp[i], self.zp[i],
+ &self.tree.root_grids[j])
for i in range(self.num_points) :
pt_grids[i] = self.point_grids[i]
@@ -221,58 +213,50 @@
@cython.boundscheck(False)
@cython.wraparound(False)
- cdef void check_positions(self, GridTreeNode * grid) :
+ cdef void check_position(self,
+ np.int64_t pt_index,
+ np.float64_t x,
+ np.float64_t y,
+ np.float64_t z,
+ GridTreeNode * grid) :
cdef int i
-
- if self.num_idxs_left == 0 : return
+ cdef np.uint8_t in_grid
+
+ in_grid = self.is_in_grid(x, y, z, grid)
- if grid.num_children > 0 :
+ if in_grid :
- for i in range(grid.num_children) :
-
- self.check_positions(grid.children[i])
-
+ if grid.num_children > 0 :
+
+ for i in range(grid.num_children) :
+
+ self.check_position(pt_index, x, y, z, grid.children[i])
+
+ else :
+
+ self.point_grids[pt_index] = grid.index
+
else :
- self.is_in_grid(grid)
-
- for i in range(self.num_points) :
-
- if self.in_grid[i] == 1 :
-
- #print self.num_idxs_left
-
- self.point_grids[i] = grid.index
- self.all_idxs[i] = -1
- self.num_idxs_left -= 1
-
return
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef void is_in_grid(self, GridTreeNode * grid) :
+ cdef np.uint8_t is_in_grid(self,
+ np.float64_t x,
+ np.float64_t y,
+ np.float64_t z,
+ GridTreeNode * grid) :
- cdef int i
cdef np.uint8_t xcond, ycond, zcond, cond
-
- for i in range(self.num_points) :
+
+ xcond = x >= grid.left_edge[0] and x < grid.right_edge[0]
+ ycond = y >= grid.left_edge[1] and y < grid.right_edge[1]
+ zcond = z >= grid.left_edge[2] and z < grid.right_edge[2]
+
+ cond = xcond and ycond
+ cond = cond and zcond
- if self.all_idxs[i] >= 0 :
-
- xcond = self.xp[i] >= grid.left_edge[0] and self.xp[i] < grid.right_edge[0]
- ycond = self.yp[i] >= grid.left_edge[1] and self.yp[i] < grid.right_edge[1]
- zcond = self.zp[i] >= grid.left_edge[2] and self.zp[i] < grid.right_edge[2]
-
- cond = xcond and ycond
- cond = cond and zcond
-
- if cond :
- self.in_grid[i] = 1
- else :
- self.in_grid[i] = 0
-
- else :
-
- self.in_grid[i] = 0
+ return cond
https://bitbucket.org/yt_analysis/yt/changeset/cc4fffdea959/
changeset: cc4fffdea959
branch: yt
user: jzuhone
date: 2012-11-20 05:48:29
summary: Making the Cythonized GridTree the way we find points
affected #: 2 files
diff -r 3283e02e6f557a33094cb58cbc3246b0dca75b59 -r cc4fffdea95956df5c31b4329f989542b124520d yt/data_objects/object_finding_mixin.py
--- a/yt/data_objects/object_finding_mixin.py
+++ b/yt/data_objects/object_finding_mixin.py
@@ -29,8 +29,9 @@
from yt.utilities.lib import \
get_box_grids_level, \
get_box_grids_below_level
-from yt.utilities.particle_utils import \
- MatchPointsToGrids
+from yt.utilities.lib import \
+ MatchPointsToGrids, \
+ GridTree
class ObjectFindingMixin(object) :
@@ -116,8 +117,10 @@
"""
Returns the (objects, indices) of leaf grids containing a number of (x,y,z) points
"""
- find_pts = MatchPointsToGrids(self,x,y,z)
- ind = find_pts()
+ num_points = len(x)
+ grid_tree = self.get_grid_tree()
+ pts = MatchPointsToGrids(grid_tree,num_points,x,y,z)
+ ind = pts.find_points_in_tree()
return self.grids[ind], ind
def find_field_value_at_point(self, fields, coord):
@@ -249,3 +252,24 @@
mask[gi] = True
return self.grids[mask], np.where(mask)
+ def get_grid_tree(self) :
+
+ left_edge = np.zeros((self.num_grids, 3))
+ right_edge = np.zeros((self.num_grids, 3))
+ level = np.zeros((self.num_grids), dtype='int64')
+ parent_ind = np.zeros((self.num_grids), dtype='int64')
+ num_children = np.zeros((self.num_grids), dtype='int64')
+
+ for i, grid in enumerate(self.grids) :
+
+ left_edge[i,:] = grid.LeftEdge
+ right_edge[i,:] = grid.RightEdge
+ level[i] = grid.Level
+ if grid.Parent is None :
+ parent_ind[i] = -1
+ else :
+ parent_ind[i] = grid.Parent.id - grid.Parent._id_offset
+ num_children[i] = np.int64(len(grid.Children))
+
+ return GridTree(self.num_grids, left_edge, right_edge, parent_ind,
+ level, num_children)
diff -r 3283e02e6f557a33094cb58cbc3246b0dca75b59 -r cc4fffdea95956df5c31b4329f989542b124520d yt/utilities/particle_utils.py
--- a/yt/utilities/particle_utils.py
+++ /dev/null
@@ -1,95 +0,0 @@
-"""
-Utilities for particles.
-
-Author: John ZuHone <jzuhone at gmail.com>
-Affiliation: NASA/Goddard Space Flight Center
-Homepage: http://yt-project.org/
-License:
-Copyright (C) 2012 John ZuHone. 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 *
-
-class MatchPointsToGrids(object) :
-
- def __init__(self, h, x, y, z) :
-
- self.num_points = len(x)
-
- self.x = x
- self.y = y
- self.z = z
-
- self.point_grids = -np.ones((self.num_points)).astype("int32")
- self.all_idxs = np.arange((self.num_points))
-
- self.root_grids = h.select_grids(0)
-
- self.counter = 0
-
- def __call__(self) :
- """
- Loops over every root grid and follows down each's subtree
- until it finds leaf grids, then assigns that grid id to the points
- in that grid
- """
-
- self.pbar = get_pbar("Finding Grids", self.num_points)
-
- for root_grid in self.root_grids :
-
- self.check_positions(root_grid)
-
- self.pbar.finish()
-
- return self.point_grids
-
- def check_positions(self, grid) :
- """
- Recursive function to traverse up and down the hierarchy
- """
-
- if len(self.all_idxs) == 0 : return # We've run out of points
-
- if len(grid.Children) > 0 :
-
- # This grid has children
-
- for child in grid.Children :
-
- self.check_positions(child)
-
- else :
-
- # We hit a leaf grid
-
- local_idxs = grid.is_in_grid(self.x[self.all_idxs], self.y[self.all_idxs],
- self.z[self.all_idxs])
-
- if np.any(local_idxs) :
-
- self.point_grids[self.all_idxs[local_idxs]] = grid.id-grid._id_offset
- self.all_idxs = self.all_idxs[np.logical_not(local_idxs)]
-
- self.counter += local_idxs.sum()
-
- self.pbar.update(self.counter)
-
- return
-
https://bitbucket.org/yt_analysis/yt/changeset/2f8e236ca52d/
changeset: 2f8e236ca52d
branch: yt
user: jzuhone
date: 2012-11-20 17:36:55
summary: Generalization to hierarchies which are patch-based.
affected #: 1 file
diff -r cc4fffdea95956df5c31b4329f989542b124520d -r 2f8e236ca52d6a4f0f6bfb4ede8276e6b35d5055 yt/utilities/lib/GridTree.pyx
--- a/yt/utilities/lib/GridTree.pyx
+++ b/yt/utilities/lib/GridTree.pyx
@@ -196,15 +196,19 @@
cdef np.ndarray[np.int64_t, ndim=1] pt_grids
cdef int i, j
-
+ cdef np.uint8_t in_grid
+
pt_grids = np.zeros(self.num_points, dtype='int64')
for i in range(self.num_points) :
+ in_grid = 0
+
for j in range(self.tree.num_root_grids) :
- self.check_position(i, self.xp[i], self.yp[i], self.zp[i],
- &self.tree.root_grids[j])
+ if not in_grid :
+ in_grid = self.check_position(i, self.xp[i], self.yp[i], self.zp[i],
+ &self.tree.root_grids[j])
for i in range(self.num_points) :
pt_grids[i] = self.point_grids[i]
@@ -213,12 +217,12 @@
@cython.boundscheck(False)
@cython.wraparound(False)
- cdef void check_position(self,
- np.int64_t pt_index,
- np.float64_t x,
- np.float64_t y,
- np.float64_t z,
- GridTreeNode * grid) :
+ cdef np.uint8_t check_position(self,
+ np.int64_t pt_index,
+ np.float64_t x,
+ np.float64_t y,
+ np.float64_t z,
+ GridTreeNode * grid) :
cdef int i
cdef np.uint8_t in_grid
@@ -229,17 +233,22 @@
if grid.num_children > 0 :
+ in_grid = 0
+
for i in range(grid.num_children) :
- self.check_position(pt_index, x, y, z, grid.children[i])
+ if not in_grid :
+ in_grid = self.check_position(pt_index, x, y, z, grid.children[i])
+
+ if not in_grid :
+ self.point_grids[pt_index] = grid.index
+
else :
self.point_grids[pt_index] = grid.index
- else :
-
- return
+ return in_grid
@cython.boundscheck(False)
@cython.wraparound(False)
https://bitbucket.org/yt_analysis/yt/changeset/7e31170e0544/
changeset: 7e31170e0544
branch: yt
user: jzuhone
date: 2012-11-20 18:55:46
summary: Verified that this works for patch-based refinement schemes.
affected #: 3 files
diff -r 2f8e236ca52d6a4f0f6bfb4ede8276e6b35d5055 -r 7e31170e054436aeda35a75bb243f61f81ff1acb yt/utilities/initial_conditions.py
--- a/yt/utilities/initial_conditions.py
+++ b/yt/utilities/initial_conditions.py
@@ -67,6 +67,47 @@
val = ((r[ind] - cr2) / (r2 - cr2))**0.5 * (outer_val - inner_val)
grid[field][ind] = val + inner_val
+class BetaModelSphere(FluidOperator):
+ def __init__(self, beta, core_radius, radius, center, fields):
+ self.radius = radius
+ self.center = center
+ self.fields = fields
+ self.core_radius = core_radius
+ self.beta = beta
+
+ def __call__(self, grid, sub_select = None):
+ r = np.zeros(grid.ActiveDimensions, dtype="float64")
+ r2 = self.radius**2
+ cr2 = self.core_radius**2
+ for i, ax in enumerate("xyz"):
+ np.add(r, (grid[ax] - self.center[i])**2.0, r)
+ ind = (r <= r2)
+ if sub_select is not None:
+ ind &= sub_select
+ for field, core_val in self.fields.iteritems() :
+ val = core_val*(1.+r[ind]/cr2)**(-1.5*self.beta)
+ grid[field][ind] = val
+
+class NFWModelSphere(FluidOperator):
+ def __init__(self, scale_radius, radius, center, fields):
+ self.radius = radius
+ self.center = center
+ self.fields = fields
+ self.scale_radius = scale_radius
+
+ def __call__(self, grid, sub_select = None):
+ r = np.zeros(grid.ActiveDimensions, dtype="float64")
+ for i, ax in enumerate("xyz"):
+ np.add(r, (grid[ax] - self.center[i])**2.0, r)
+ np,sqrt(r,r)
+ ind = (r <= self.radius)
+ r /= scale.radius
+ if sub_select is not None:
+ ind &= sub_select
+ for field, scale_val in self.fields.iteritems() :
+ val = scale_val/(r[ind]*(1.+r[ind])**2)
+ grid[field][ind] = val
+
class RandomFluctuation(FluidOperator):
def __init__(self, fields):
self.fields = fields
diff -r 2f8e236ca52d6a4f0f6bfb4ede8276e6b35d5055 -r 7e31170e054436aeda35a75bb243f61f81ff1acb yt/utilities/lib/CICDeposit.pyx
--- a/yt/utilities/lib/CICDeposit.pyx
+++ b/yt/utilities/lib/CICDeposit.pyx
@@ -116,6 +116,66 @@
@cython.boundscheck(False)
@cython.wraparound(False)
+def CICSample_3(np.ndarray[np.float64_t, ndim=1] posx,
+ np.ndarray[np.float64_t, ndim=1] posy,
+ np.ndarray[np.float64_t, ndim=1] posz,
+ np.ndarray[np.float64_t, ndim=1] sample,
+ np.int64_t npositions,
+ np.ndarray[np.float64_t, ndim=3] field,
+ np.ndarray[np.float64_t, ndim=1] leftEdge,
+ np.ndarray[np.int32_t, ndim=1] gridDimension,
+ np.float64_t cellSize):
+
+ cdef int i1, j1, k1, n
+ cdef double xpos, ypos, zpos
+ cdef double fact, edge0, edge1, edge2
+ cdef double le0, le1, le2
+ cdef float dx, dy, dz, dx2, dy2, dz2
+
+ # Note the accounting for the one layer of ghost cells!
+
+ edge0 = (<float> gridDimension[0]+2) - 0.5001
+ edge1 = (<float> gridDimension[1]+2) - 0.5001
+ edge2 = (<float> gridDimension[2]+2) - 0.5001
+ fact = 1.0 / cellSize
+
+ # Note the accounting for the one layer of ghost cells!
+
+ le0 = leftEdge[0] - cellSize
+ le1 = leftEdge[1] - cellSize
+ le2 = leftEdge[2] - cellSize
+
+ for n in range(npositions):
+
+ # Compute the position of the central cell
+ xpos = fmin(fmax((posx[n] - le0)*fact, 0.5001), edge0)
+ ypos = fmin(fmax((posy[n] - le1)*fact, 0.5001), edge1)
+ zpos = fmin(fmax((posz[n] - le2)*fact, 0.5001), edge2)
+
+ i1 = <int> (xpos + 0.5)
+ j1 = <int> (ypos + 0.5)
+ k1 = <int> (zpos + 0.5)
+
+ # Compute the weights
+ dx = (<float> i1) + 0.5 - xpos
+ dy = (<float> j1) + 0.5 - ypos
+ dz = (<float> k1) + 0.5 - zpos
+ dx2 = 1.0 - dx
+ dy2 = 1.0 - dy
+ dz2 = 1.0 - dz
+
+ # Interpolate from field onto the particle
+ sample[n] = (field[i1-1,j1-1,k1-1] * dx * dy * dz +
+ field[i1 ,j1-1,k1-1] * dx2 * dy * dz +
+ field[i1-1,j1 ,k1-1] * dx * dy2 * dz +
+ field[i1 ,j1 ,k1-1] * dx2 * dy2 * dz +
+ field[i1-1,j1-1,k1 ] * dx * dy * dz2 +
+ field[i1 ,j1-1,k1 ] * dx2 * dy * dz2 +
+ field[i1-1,j1 ,k1 ] * dx * dy2 * dz2 +
+ field[i1 ,j1 ,k1 ] * dx2 * dy2 * dz2)
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
@cython.cdivision(True)
def assign_particles_to_cells(np.ndarray[np.int32_t, ndim=1] levels, #for cells
np.ndarray[np.float32_t, ndim=2] left_edges, #many cells
diff -r 2f8e236ca52d6a4f0f6bfb4ede8276e6b35d5055 -r 7e31170e054436aeda35a75bb243f61f81ff1acb yt/utilities/lib/GridTree.pyx
--- a/yt/utilities/lib/GridTree.pyx
+++ b/yt/utilities/lib/GridTree.pyx
@@ -243,11 +243,13 @@
if not in_grid :
self.point_grids[pt_index] = grid.index
-
+ in_grid = 1
+
else :
self.point_grids[pt_index] = grid.index
-
+ in_grid = 1
+
return in_grid
@cython.boundscheck(False)
https://bitbucket.org/yt_analysis/yt/changeset/6d6b74fd852f/
changeset: 6d6b74fd852f
branch: yt
user: jzuhone
date: 2012-11-20 15:58:08
summary: Merging
affected #: 3 files
diff -r cc4fffdea95956df5c31b4329f989542b124520d -r 6d6b74fd852f78325179288050bfdbaf5c58b48b yt/frontends/enzo/answer_testing_support.py
--- a/yt/frontends/enzo/answer_testing_support.py
+++ b/yt/frontends/enzo/answer_testing_support.py
@@ -35,7 +35,8 @@
GridValuesTest, \
ProjectionValuesTest, \
ParentageRelationshipsTest, \
- temp_cwd
+ temp_cwd, \
+ AssertWrapper
def requires_outputlog(path = ".", prefix = ""):
def ffalse(func):
@@ -94,9 +95,10 @@
for xmin, xmax in zip(self.left_edges, self.right_edges):
mask = (position >= xmin)*(position <= xmax)
exact_field = np.interp(position[mask], exact['pos'], exact[k])
+ myname = "ShockTubeTest_%s" % k
# yield test vs analytical solution
- yield assert_allclose, field[mask], exact_field, \
- self.rtol, self.atol
+ yield AssertWrapper(myname, assert_allclose, field[mask],
+ exact_field, self.rtol, self.atol)
def get_analytical_solution(self):
# Reads in from file
diff -r cc4fffdea95956df5c31b4329f989542b124520d -r 6d6b74fd852f78325179288050bfdbaf5c58b48b yt/frontends/orion/fields.py
--- a/yt/frontends/orion/fields.py
+++ b/yt/frontends/orion/fields.py
@@ -126,7 +126,7 @@
data["x-velocity"]**2.0
+ data["y-velocity"]**2.0
+ data["z-velocity"]**2.0 )
-add_orion_field("ThermalEnergy", function=_ThermalEnergy,
+add_field("ThermalEnergy", function=_ThermalEnergy,
units=r"\rm{ergs}/\rm{cm^3}")
def _Pressure(field,data):
@@ -134,11 +134,11 @@
NB: this will need to be modified for radiation
"""
return (data.pf["Gamma"] - 1.0)*data["ThermalEnergy"]
-add_orion_field("Pressure", function=_Pressure, units=r"\rm{dyne}/\rm{cm}^{2}")
+add_field("Pressure", function=_Pressure, units=r"\rm{dyne}/\rm{cm}^{2}")
def _Temperature(field,data):
return (data.pf["Gamma"]-1.0)*data.pf["mu"]*mh*data["ThermalEnergy"]/(kboltz*data["Density"])
-add_orion_field("Temperature",function=_Temperature,units=r"\rm{Kelvin}",take_log=False)
+add_field("Temperature",function=_Temperature,units=r"\rm{Kelvin}",take_log=False)
# particle fields
diff -r cc4fffdea95956df5c31b4329f989542b124520d -r 6d6b74fd852f78325179288050bfdbaf5c58b48b yt/utilities/answer_testing/framework.py
--- a/yt/utilities/answer_testing/framework.py
+++ b/yt/utilities/answer_testing/framework.py
@@ -112,7 +112,7 @@
(os.path.realpath(options.output_dir),
options.store_name)
if not os.path.isdir(name_dir_path):
- os.mkdir(name_dir_path)
+ os.makedirs(name_dir_path)
options.store_name= "%s/%s" % \
(name_dir_path, options.store_name)
else:
@@ -554,3 +554,18 @@
yield PixelizedProjectionValuesTest(
pf_fn, axis, field, weight_field,
ds)
+
+class AssertWrapper(object):
+ """
+ Used to wrap a numpy testing assertion, in order to provide a useful name
+ for a given assertion test.
+ """
+ def __init__(self, description, *args):
+ # The key here is to add a description attribute, which nose will pick
+ # up.
+ self.args = args
+ self.description = description
+
+ def __call__(self):
+ self.args[0](*self.args[1:])
+
https://bitbucket.org/yt_analysis/yt/changeset/852a6d0b3ac3/
changeset: 852a6d0b3ac3
branch: yt
user: jzuhone
date: 2012-11-20 17:37:27
summary: Merging
affected #: 1 file
diff -r 6d6b74fd852f78325179288050bfdbaf5c58b48b -r 852a6d0b3ac32e319d5ac9eb4b0e9f328bf2cd86 yt/utilities/lib/GridTree.pyx
--- a/yt/utilities/lib/GridTree.pyx
+++ b/yt/utilities/lib/GridTree.pyx
@@ -196,15 +196,19 @@
cdef np.ndarray[np.int64_t, ndim=1] pt_grids
cdef int i, j
-
+ cdef np.uint8_t in_grid
+
pt_grids = np.zeros(self.num_points, dtype='int64')
for i in range(self.num_points) :
+ in_grid = 0
+
for j in range(self.tree.num_root_grids) :
- self.check_position(i, self.xp[i], self.yp[i], self.zp[i],
- &self.tree.root_grids[j])
+ if not in_grid :
+ in_grid = self.check_position(i, self.xp[i], self.yp[i], self.zp[i],
+ &self.tree.root_grids[j])
for i in range(self.num_points) :
pt_grids[i] = self.point_grids[i]
@@ -213,12 +217,12 @@
@cython.boundscheck(False)
@cython.wraparound(False)
- cdef void check_position(self,
- np.int64_t pt_index,
- np.float64_t x,
- np.float64_t y,
- np.float64_t z,
- GridTreeNode * grid) :
+ cdef np.uint8_t check_position(self,
+ np.int64_t pt_index,
+ np.float64_t x,
+ np.float64_t y,
+ np.float64_t z,
+ GridTreeNode * grid) :
cdef int i
cdef np.uint8_t in_grid
@@ -229,17 +233,22 @@
if grid.num_children > 0 :
+ in_grid = 0
+
for i in range(grid.num_children) :
- self.check_position(pt_index, x, y, z, grid.children[i])
+ if not in_grid :
+ in_grid = self.check_position(pt_index, x, y, z, grid.children[i])
+
+ if not in_grid :
+ self.point_grids[pt_index] = grid.index
+
else :
self.point_grids[pt_index] = grid.index
- else :
-
- return
+ return in_grid
@cython.boundscheck(False)
@cython.wraparound(False)
https://bitbucket.org/yt_analysis/yt/changeset/967446e53154/
changeset: 967446e53154
branch: yt
user: jzuhone
date: 2012-11-20 18:58:21
summary: Merging
affected #: 3 files
diff -r 7e31170e054436aeda35a75bb243f61f81ff1acb -r 967446e53154a31cdd178f8b45e740dca2bc8668 yt/frontends/enzo/answer_testing_support.py
--- a/yt/frontends/enzo/answer_testing_support.py
+++ b/yt/frontends/enzo/answer_testing_support.py
@@ -35,7 +35,8 @@
GridValuesTest, \
ProjectionValuesTest, \
ParentageRelationshipsTest, \
- temp_cwd
+ temp_cwd, \
+ AssertWrapper
def requires_outputlog(path = ".", prefix = ""):
def ffalse(func):
@@ -94,9 +95,10 @@
for xmin, xmax in zip(self.left_edges, self.right_edges):
mask = (position >= xmin)*(position <= xmax)
exact_field = np.interp(position[mask], exact['pos'], exact[k])
+ myname = "ShockTubeTest_%s" % k
# yield test vs analytical solution
- yield assert_allclose, field[mask], exact_field, \
- self.rtol, self.atol
+ yield AssertWrapper(myname, assert_allclose, field[mask],
+ exact_field, self.rtol, self.atol)
def get_analytical_solution(self):
# Reads in from file
diff -r 7e31170e054436aeda35a75bb243f61f81ff1acb -r 967446e53154a31cdd178f8b45e740dca2bc8668 yt/frontends/orion/fields.py
--- a/yt/frontends/orion/fields.py
+++ b/yt/frontends/orion/fields.py
@@ -126,7 +126,7 @@
data["x-velocity"]**2.0
+ data["y-velocity"]**2.0
+ data["z-velocity"]**2.0 )
-add_orion_field("ThermalEnergy", function=_ThermalEnergy,
+add_field("ThermalEnergy", function=_ThermalEnergy,
units=r"\rm{ergs}/\rm{cm^3}")
def _Pressure(field,data):
@@ -134,11 +134,11 @@
NB: this will need to be modified for radiation
"""
return (data.pf["Gamma"] - 1.0)*data["ThermalEnergy"]
-add_orion_field("Pressure", function=_Pressure, units=r"\rm{dyne}/\rm{cm}^{2}")
+add_field("Pressure", function=_Pressure, units=r"\rm{dyne}/\rm{cm}^{2}")
def _Temperature(field,data):
return (data.pf["Gamma"]-1.0)*data.pf["mu"]*mh*data["ThermalEnergy"]/(kboltz*data["Density"])
-add_orion_field("Temperature",function=_Temperature,units=r"\rm{Kelvin}",take_log=False)
+add_field("Temperature",function=_Temperature,units=r"\rm{Kelvin}",take_log=False)
# particle fields
diff -r 7e31170e054436aeda35a75bb243f61f81ff1acb -r 967446e53154a31cdd178f8b45e740dca2bc8668 yt/utilities/answer_testing/framework.py
--- a/yt/utilities/answer_testing/framework.py
+++ b/yt/utilities/answer_testing/framework.py
@@ -112,7 +112,7 @@
(os.path.realpath(options.output_dir),
options.store_name)
if not os.path.isdir(name_dir_path):
- os.mkdir(name_dir_path)
+ os.makedirs(name_dir_path)
options.store_name= "%s/%s" % \
(name_dir_path, options.store_name)
else:
@@ -554,3 +554,18 @@
yield PixelizedProjectionValuesTest(
pf_fn, axis, field, weight_field,
ds)
+
+class AssertWrapper(object):
+ """
+ Used to wrap a numpy testing assertion, in order to provide a useful name
+ for a given assertion test.
+ """
+ def __init__(self, description, *args):
+ # The key here is to add a description attribute, which nose will pick
+ # up.
+ self.args = args
+ self.description = description
+
+ def __call__(self):
+ self.args[0](*self.args[1:])
+
https://bitbucket.org/yt_analysis/yt/changeset/82c32faef934/
changeset: 82c32faef934
branch: yt
user: jzuhone
date: 2012-11-20 19:04:37
summary: Undoing a change that wasn't ready yet
affected #: 2 files
diff -r 967446e53154a31cdd178f8b45e740dca2bc8668 -r 82c32faef93404f2457578e97a1e6d9bd2675181 yt/utilities/initial_conditions.py
--- a/yt/utilities/initial_conditions.py
+++ b/yt/utilities/initial_conditions.py
@@ -67,47 +67,6 @@
val = ((r[ind] - cr2) / (r2 - cr2))**0.5 * (outer_val - inner_val)
grid[field][ind] = val + inner_val
-class BetaModelSphere(FluidOperator):
- def __init__(self, beta, core_radius, radius, center, fields):
- self.radius = radius
- self.center = center
- self.fields = fields
- self.core_radius = core_radius
- self.beta = beta
-
- def __call__(self, grid, sub_select = None):
- r = np.zeros(grid.ActiveDimensions, dtype="float64")
- r2 = self.radius**2
- cr2 = self.core_radius**2
- for i, ax in enumerate("xyz"):
- np.add(r, (grid[ax] - self.center[i])**2.0, r)
- ind = (r <= r2)
- if sub_select is not None:
- ind &= sub_select
- for field, core_val in self.fields.iteritems() :
- val = core_val*(1.+r[ind]/cr2)**(-1.5*self.beta)
- grid[field][ind] = val
-
-class NFWModelSphere(FluidOperator):
- def __init__(self, scale_radius, radius, center, fields):
- self.radius = radius
- self.center = center
- self.fields = fields
- self.scale_radius = scale_radius
-
- def __call__(self, grid, sub_select = None):
- r = np.zeros(grid.ActiveDimensions, dtype="float64")
- for i, ax in enumerate("xyz"):
- np.add(r, (grid[ax] - self.center[i])**2.0, r)
- np,sqrt(r,r)
- ind = (r <= self.radius)
- r /= scale.radius
- if sub_select is not None:
- ind &= sub_select
- for field, scale_val in self.fields.iteritems() :
- val = scale_val/(r[ind]*(1.+r[ind])**2)
- grid[field][ind] = val
-
class RandomFluctuation(FluidOperator):
def __init__(self, fields):
self.fields = fields
diff -r 967446e53154a31cdd178f8b45e740dca2bc8668 -r 82c32faef93404f2457578e97a1e6d9bd2675181 yt/utilities/lib/CICDeposit.pyx
--- a/yt/utilities/lib/CICDeposit.pyx
+++ b/yt/utilities/lib/CICDeposit.pyx
@@ -116,66 +116,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
-def CICSample_3(np.ndarray[np.float64_t, ndim=1] posx,
- np.ndarray[np.float64_t, ndim=1] posy,
- np.ndarray[np.float64_t, ndim=1] posz,
- np.ndarray[np.float64_t, ndim=1] sample,
- np.int64_t npositions,
- np.ndarray[np.float64_t, ndim=3] field,
- np.ndarray[np.float64_t, ndim=1] leftEdge,
- np.ndarray[np.int32_t, ndim=1] gridDimension,
- np.float64_t cellSize):
-
- cdef int i1, j1, k1, n
- cdef double xpos, ypos, zpos
- cdef double fact, edge0, edge1, edge2
- cdef double le0, le1, le2
- cdef float dx, dy, dz, dx2, dy2, dz2
-
- # Note the accounting for the one layer of ghost cells!
-
- edge0 = (<float> gridDimension[0]+2) - 0.5001
- edge1 = (<float> gridDimension[1]+2) - 0.5001
- edge2 = (<float> gridDimension[2]+2) - 0.5001
- fact = 1.0 / cellSize
-
- # Note the accounting for the one layer of ghost cells!
-
- le0 = leftEdge[0] - cellSize
- le1 = leftEdge[1] - cellSize
- le2 = leftEdge[2] - cellSize
-
- for n in range(npositions):
-
- # Compute the position of the central cell
- xpos = fmin(fmax((posx[n] - le0)*fact, 0.5001), edge0)
- ypos = fmin(fmax((posy[n] - le1)*fact, 0.5001), edge1)
- zpos = fmin(fmax((posz[n] - le2)*fact, 0.5001), edge2)
-
- i1 = <int> (xpos + 0.5)
- j1 = <int> (ypos + 0.5)
- k1 = <int> (zpos + 0.5)
-
- # Compute the weights
- dx = (<float> i1) + 0.5 - xpos
- dy = (<float> j1) + 0.5 - ypos
- dz = (<float> k1) + 0.5 - zpos
- dx2 = 1.0 - dx
- dy2 = 1.0 - dy
- dz2 = 1.0 - dz
-
- # Interpolate from field onto the particle
- sample[n] = (field[i1-1,j1-1,k1-1] * dx * dy * dz +
- field[i1 ,j1-1,k1-1] * dx2 * dy * dz +
- field[i1-1,j1 ,k1-1] * dx * dy2 * dz +
- field[i1 ,j1 ,k1-1] * dx2 * dy2 * dz +
- field[i1-1,j1-1,k1 ] * dx * dy * dz2 +
- field[i1 ,j1-1,k1 ] * dx2 * dy * dz2 +
- field[i1-1,j1 ,k1 ] * dx * dy2 * dz2 +
- field[i1 ,j1 ,k1 ] * dx2 * dy2 * dz2)
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
@cython.cdivision(True)
def assign_particles_to_cells(np.ndarray[np.int32_t, ndim=1] levels, #for cells
np.ndarray[np.float32_t, ndim=2] left_edges, #many cells
https://bitbucket.org/yt_analysis/yt/changeset/05d098ee846f/
changeset: 05d098ee846f
branch: yt
user: jzuhone
date: 2012-11-20 20:28:11
summary: Tests for GridTree
affected #: 1 file
diff -r 82c32faef93404f2457578e97a1e6d9bd2675181 -r 05d098ee846f9cb7dad14515bbfc633d00a40dce yt/utilities/lib/tests/test_grid_tree.py
--- /dev/null
+++ b/yt/utilities/lib/tests/test_grid_tree.py
@@ -0,0 +1,75 @@
+import numpy as np
+
+from yt.testing import *
+from yt.frontends.stream.api import load_amr_grids
+
+def setup():
+
+ global pf
+
+ grid_data = [
+ dict(left_edge = [0.0, 0.0, 0.0], right_edge = [1.0, 1.0, 1.],
+ level = 0, dimensions = [16, 16, 16]),
+ dict(left_edge = [0.25, 0.25, 0.25], right_edge = [0.75, 0.75, 0.75],
+ level = 1, dimensions = [16, 16, 16]),
+ dict(left_edge = [0.25, 0.25, 0.375], right_edge = [0.5, 0.5, 0.625],
+ level = 2, dimensions = [16, 16, 16]),
+ dict(left_edge = [0.5, 0.5, 0.375], right_edge = [0.75, 0.75, 0.625],
+ level = 2, dimensions = [16, 16, 16]),
+ dict(left_edge = [0.3125, 0.3125, 0.4375], right_edge = [0.4375, 0.4375, 0.5625],
+ level = 3, dimensions = [16, 16, 16]),
+ dict(left_edge = [0.5625, 0.5625, 0.4375], right_edge = [0.6875, 0.6875, 0.5625],
+ level = 3, dimensions = [16, 16, 16])
+ ]
+
+ for g in grid_data: g["Density"] = np.random.random(g["dimensions"]) * 2**g["level"]
+ pf = load_amr_grids(grid_data, [16, 16, 16], 1.0)
+
+def test_grid_tree() :
+
+ grid_tree = pf.h.get_grid_tree()
+ indices, levels, nchild, children = grid_tree.return_tree_info()
+
+ grid_levels = [grid.Level for grid in pf.h.grids]
+ grid_indices = [grid.id-grid._id_offset for grid in pf.h.grids]
+ grid_nchild = [len(grid.Children) for grid in pf.h.grids]
+
+ print levels, grid_levels
+ assert_equal(levels, grid_levels)
+ assert_equal(indices, grid_indices)
+ assert_equal(nchild, grid_nchild)
+
+ for i, grid in enumerate(pf.h.grids) :
+ if grid_nchild[i] > 0:
+ grid_children = np.array([child.id-child._id_offset
+ for child in grid.Children])
+ assert_equal(grid_children, children[i])
+
+def test_find_points() :
+
+ num_points = 100
+
+ x = np.random.uniform(low=pf.domain_left_edge[0],
+ high=pf.domain_right_edge[0], size=num_points)
+ y = np.random.uniform(low=pf.domain_left_edge[1],
+ high=pf.domain_right_edge[1], size=num_points)
+ z = np.random.uniform(low=pf.domain_left_edge[2],
+ high=pf.domain_right_edge[2], size=num_points)
+
+ point_grids, point_grid_inds = pf.h.find_points(x,y,z)
+
+ grid_inds = np.zeros((num_points), dtype='int64')
+
+ for i, xx, yy, zz in zip(range(num_points), x, y, z) :
+
+ pt_level = -1
+
+ for grid in pf.h.grids:
+
+ if grid.is_in_grid(xx, yy, zz) :
+
+ if grid.Level > pt_level :
+ pt_level = grid.Level
+ grid_inds[i] = grid.id-grid._id_offset
+
+ assert_equal(point_grid_inds, grid_inds)
https://bitbucket.org/yt_analysis/yt/changeset/26e7f66ee689/
changeset: 26e7f66ee689
branch: yt
user: jzuhone
date: 2012-11-20 21:16:55
summary: Merging
affected #: 3 files
diff -r 05d098ee846f9cb7dad14515bbfc633d00a40dce -r 26e7f66ee6899c3cc8b57730c5d235f9508b89b3 yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -30,6 +30,7 @@
import h5py
import numpy as np
import weakref
+import glob #ST 9/12
from yt.funcs import *
from yt.data_objects.grid_patch import \
AMRGridPatch
@@ -112,6 +113,8 @@
grid['read_field'] = field
grid['read_type'] = 'vector'
+
+
class AthenaHierarchy(AMRHierarchy):
grid = AthenaGrid
@@ -215,34 +218,88 @@
(np.prod(grid['dimensions']), grid['ncells']))
raise TypeError
+ # Need to determine how many grids: self.num_grids
+ dname = self.hierarchy_filename
+ gridlistread = glob.glob('id*/%s-id*%s' % (dname[4:-9],dname[-9:] ))
+ gridlistread.insert(0,self.hierarchy_filename)
+ self.num_grids = len(gridlistread)
dxs=[]
self.grids = np.empty(self.num_grids, dtype='object')
levels = np.zeros(self.num_grids, dtype='int32')
- single_grid_width = grid['dds']*grid['dimensions']
- grids_per_dim = (self.parameter_file.domain_width/single_grid_width).astype('int32')
- glis = np.empty((self.num_grids,3), dtype='int64')
- for i in range(self.num_grids):
- procz = i/(grids_per_dim[0]*grids_per_dim[1])
- procy = (i - procz*(grids_per_dim[0]*grids_per_dim[1]))/grids_per_dim[0]
- procx = i - procz*(grids_per_dim[0]*grids_per_dim[1]) - procy*grids_per_dim[0]
- glis[i, 0] = procx*grid['dimensions'][0]
- glis[i, 1] = procy*grid['dimensions'][1]
- glis[i, 2] = procz*grid['dimensions'][2]
+ glis = np.empty((self.num_grids,3), dtype='float64')
+ gdds = np.empty((self.num_grids,3), dtype='float64')
gdims = np.ones_like(glis)
- gdims[:] = grid['dimensions']
+ j = 0
+ while j < (self.num_grids):
+ f = open(gridlistread[j],'rb')
+ f.close()
+ if j == 0:
+ f = open(dname,'rb')
+ if j != 0:
+ f = open('id%i/%s-id%i%s' % (j, dname[4:-9],j, dname[-9:]),'rb')
+ gridread = {}
+ gridread['read_field'] = None
+ gridread['read_type'] = None
+ table_read=False
+ line = f.readline()
+ while gridread['read_field'] is None:
+ parse_line(line, gridread)
+ if "SCALAR" in line.strip().split():
+ break
+ if "VECTOR" in line.strip().split():
+ break
+ if 'TABLE' in line.strip().split():
+ break
+ if len(line) == 0: break
+ line = f.readline()
+ f.close()
+ glis[j,0] = gridread['left_edge'][0]
+ glis[j,1] = gridread['left_edge'][1]
+ glis[j,2] = gridread['left_edge'][2]
+ # It seems some datasets have a mismatch between ncells and
+ # the actual grid dimensions.
+ if np.prod(gridread['dimensions']) != gridread['ncells']:
+ gridread['dimensions'] -= 1
+ gridread['dimensions'][gridread['dimensions']==0]=1
+ if np.prod(gridread['dimensions']) != gridread['ncells']:
+ mylog.error('product of dimensions %i not equal to number of cells %i' %
+ (np.prod(gridread['dimensions']), gridread['ncells']))
+ raise TypeError
+ gdims[j,0] = gridread['dimensions'][0]
+ gdims[j,1] = gridread['dimensions'][1]
+ gdims[j,2] = gridread['dimensions'][2]
+ gdds[j,:] = gridread['dds']
+
+ j=j+1
+
+ gres = glis + gdims*gdds
+ # Now we convert the glis, which were left edges (floats), to indices
+ # from the domain left edge. Then we do a bunch of fixing now that we
+ # know the extent of all the grids.
+ glis = np.round((glis - self.parameter_file.domain_left_edge)/gdds).astype('int')
+ new_dre = np.max(gres,axis=0)
+ self.parameter_file.domain_right_edge = np.round(new_dre, decimals=6)
+ self.parameter_file.domain_width = \
+ (self.parameter_file.domain_right_edge -
+ self.parameter_file.domain_left_edge)
+ self.parameter_file.domain_center = \
+ 0.5*(self.parameter_file.domain_left_edge +
+ self.parameter_file.domain_right_edge)
+ self.parameter_file.domain_dimensions = \
+ np.round(self.parameter_file.domain_width/gdds[0]).astype('int')
for i in range(levels.shape[0]):
- self.grids[i] = self.grid(i, self, levels[i],
+ self.grids[i] = self.grid(i,self,levels[i],
glis[i],
gdims[i])
-
dx = (self.parameter_file.domain_right_edge-
self.parameter_file.domain_left_edge)/self.parameter_file.domain_dimensions
dx = dx/self.parameter_file.refine_by**(levels[i])
- dxs.append(grid['dds'])
+ dxs.append(dx)
+
dx = np.array(dxs)
- self.grid_left_edge = self.parameter_file.domain_left_edge + dx*glis
+ self.grid_left_edge = np.round(self.parameter_file.domain_left_edge + dx*glis, decimals=6)
self.grid_dimensions = gdims.astype("int32")
- self.grid_right_edge = self.grid_left_edge + dx*self.grid_dimensions
+ self.grid_right_edge = np.round(self.grid_left_edge + dx*self.grid_dimensions, decimals=6)
self.grid_particle_count = np.zeros([self.num_grids, 1], dtype='int64')
def _populate_grid_objects(self):
@@ -256,9 +313,6 @@
g1.Parent.append(g)
self.max_level = self.grid_levels.max()
-# def _setup_derived_fields(self):
-# self.derived_field_list = []
-
def _get_grid_children(self, grid):
mask = np.zeros(self.num_grids, dtype='bool')
grids, grid_ind = self.get_box_grids(grid.LeftEdge, grid.RightEdge)
@@ -277,6 +331,11 @@
StaticOutput.__init__(self, filename, data_style)
self.filename = filename
self.storage_filename = filename[4:-4]
+
+ # Unfortunately we now have to mandate that the hierarchy gets
+ # instantiated so that we can make sure we have the correct left
+ # and right domain edges.
+ self.h
def _set_units(self):
"""
@@ -315,14 +374,11 @@
line = self._handle.readline()
self.domain_left_edge = grid['left_edge']
- if 'domain_right_edge' in self.specified_parameters:
- self.domain_right_edge = np.array(self.specified_parameters['domain_right_edge'])
- else:
- mylog.info("Please set 'domain_right_edge' in parameters dictionary argument " +
- "if it is not equal to -domain_left_edge.")
- self.domain_right_edge = -self.domain_left_edge
+ mylog.info("Temporarily setting domain_right_edge = -domain_left_edge."+
+ " This will be corrected automatically if it is not the case.")
+ self.domain_right_edge = -self.domain_left_edge
self.domain_width = self.domain_right_edge-self.domain_left_edge
- self.domain_dimensions = self.domain_width/grid['dds']
+ self.domain_dimensions = np.round(self.domain_width/grid['dds']).astype('int32')
refine_by = None
if refine_by is None: refine_by = 2
self.refine_by = refine_by
@@ -341,8 +397,9 @@
self.field_ordering = 'fortran'
self.boundary_conditions = [1]*6
- ND = self.dimensionality
- self.nvtk = int(np.product(self.domain_dimensions[:ND]/(grid['dimensions'][:ND]-1)))
+ dname = self.parameter_filename
+ gridlistread = glob.glob('id*/%s-id*%s' % (dname[4:-9],dname[-9:] ))
+ self.nvtk = len(gridlistread)+1
self.current_redshift = self.omega_lambda = self.omega_matter = \
self.hubble_constant = self.cosmological_simulation = 0.0
@@ -350,6 +407,7 @@
self.parameters["HydroMethod"] = 0 # Hardcode for now until field staggering is supported.
self._handle.close()
+
@classmethod
def _is_valid(self, *args, **kwargs):
try:
diff -r 05d098ee846f9cb7dad14515bbfc633d00a40dce -r 26e7f66ee6899c3cc8b57730c5d235f9508b89b3 yt/frontends/athena/fields.py
--- a/yt/frontends/athena/fields.py
+++ b/yt/frontends/athena/fields.py
@@ -72,13 +72,13 @@
units=r"")
add_athena_field("cell_centered_B_x", function=NullFunc, take_log=False,
- units=r"")
+ units=r"", display_name=r"$\rm{cell\ centered\ B_x}$")
add_athena_field("cell_centered_B_y", function=NullFunc, take_log=False,
- units=r"")
+ units=r"", display_name=r"$\rm{cell\ centered\ B_y}$")
add_athena_field("cell_centered_B_z", function=NullFunc, take_log=False,
- units=r"")
+ units=r"", display_name=r"$\rm{cell\ centered\ B_z}$")
for f,v in log_translation_dict.items():
add_field(f, TranslationFunc(v), take_log=True)
diff -r 05d098ee846f9cb7dad14515bbfc633d00a40dce -r 26e7f66ee6899c3cc8b57730c5d235f9508b89b3 yt/frontends/athena/io.py
--- a/yt/frontends/athena/io.py
+++ b/yt/frontends/athena/io.py
@@ -45,10 +45,15 @@
def _read_data_set(self,grid,field):
f = file(grid.filename, 'rb')
- dtype, offset = grid.hierarchy._field_map[field]
+ dtype, offsetr = grid.hierarchy._field_map[field]
grid_ncells = np.prod(grid.ActiveDimensions)
grid_dims = grid.ActiveDimensions
+ grid0_ncells = np.prod(grid.hierarchy.grid_dimensions[0,:])
read_table_offset = get_read_table_offset(f)
+ if grid_ncells != grid0_ncells:
+ offset = offsetr + ((grid_ncells-grid0_ncells) * (offsetr//grid0_ncells))
+ if grid_ncells == grid0_ncells:
+ offset = offsetr
f.seek(read_table_offset+offset)
if dtype == 'scalar':
data = np.fromfile(f, dtype='>f4',
@@ -74,10 +79,15 @@
sl.reverse()
f = file(grid.filename, 'rb')
- dtype, offset = grid.hierarchy._field_map[field]
+ dtype, offsetr = grid.hierarchy._field_map[field]
grid_ncells = np.prod(grid.ActiveDimensions)
-
+ grid_dims = grid.ActiveDimensions
+ grid0_ncells = np.prod(grid.hierarchy.grid_dimensions[0,:])
read_table_offset = get_read_table_offset(f)
+ if grid_ncells != grid0_ncells:
+ offset = offsetr + ((grid_ncells-grid0_ncells) * (offsetr//grid0_ncells))
+ if grid_ncells == grid0_ncells:
+ offset = offsetr
f.seek(read_table_offset+offset)
if dtype == 'scalar':
data = np.fromfile(f, dtype='>f4',
https://bitbucket.org/yt_analysis/yt/changeset/1ea61855c847/
changeset: 1ea61855c847
branch: yt
user: jzuhone
date: 2012-11-20 21:42:49
summary: Small fixes to get 1D/2D data working with Sam's new changes.
affected #: 1 file
diff -r 26e7f66ee6899c3cc8b57730c5d235f9508b89b3 -r 1ea61855c847d65670ceb8354c2803d891ecf400 yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -268,6 +268,8 @@
gdims[j,0] = gridread['dimensions'][0]
gdims[j,1] = gridread['dimensions'][1]
gdims[j,2] = gridread['dimensions'][2]
+ # Setting dds=1 for non-active dimensions in 1D/2D datasets
+ gridread['dds'][gridread['dimensions']==1] = 1.
gdds[j,:] = gridread['dds']
j=j+1
@@ -287,6 +289,10 @@
self.parameter_file.domain_right_edge)
self.parameter_file.domain_dimensions = \
np.round(self.parameter_file.domain_width/gdds[0]).astype('int')
+ if self.parameter_file.dimensionality <= 2 :
+ self.parameter_file.domain_dimensions[2] = np.int(1)
+ if self.parameter_file.dimensionality == 1 :
+ self.parameter_file.domain_dimensions[1] = np.int(1)
for i in range(levels.shape[0]):
self.grids[i] = self.grid(i,self,levels[i],
glis[i],
@@ -387,8 +393,8 @@
dimensionality = 2
if grid['dimensions'][1] == 1 :
dimensionality = 1
- if dimensionality <= 2 : self.domain_dimensions[2] = 1.
- if dimensionality == 1 : self.domain_dimensions[1] = 1.
+ if dimensionality <= 2 : self.domain_dimensions[2] = np.int32(1)
+ if dimensionality == 1 : self.domain_dimensions[1] = np.int32(1)
self.dimensionality = dimensionality
self.current_time = grid["time"]
self.unique_identifier = self._handle.__hash__()
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
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