[yt-svn] commit/yt: ngoldbaum: Merged in MatthewTurk/yt/yt-3.0 (pull request #990)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Jul 3 10:40:29 PDT 2014
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/0050516c82a9/
Changeset: 0050516c82a9
Branch: yt-3.0
User: ngoldbaum
Date: 2014-07-03 19:40:18
Summary: Merged in MatthewTurk/yt/yt-3.0 (pull request #990)
Non-Cartesian particles and NN-deposition
Affected #: 7 files
diff -r b699b4b430a41bacab986e82c472fa2efcac02a4 -r 0050516c82a9622e381ef5ecc64ac444d64f9bf1 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -626,7 +626,7 @@
self.ActiveDimensions = np.array(dims, dtype='int32')
if self.ActiveDimensions.size == 1:
self.ActiveDimensions = np.array([dims, dims, dims], dtype="int32")
- self.dds = (self.right_edge - self.left_edge)/self.ActiveDimensions
+ self.dds = self.base_dds = (self.right_edge - self.left_edge)/self.ActiveDimensions
self.level = 99
self._setup_data_source()
diff -r b699b4b430a41bacab986e82c472fa2efcac02a4 -r 0050516c82a9622e381ef5ecc64ac444d64f9bf1 yt/data_objects/octree_subset.py
--- a/yt/data_objects/octree_subset.py
+++ b/yt/data_objects/octree_subset.py
@@ -199,7 +199,7 @@
op.process_octree(self.oct_handler, mdom_ind, positions,
self.fcoords, fields,
self.domain_id, self._domain_offset, self.pf.periodicity,
- index_fields, particle_octree, pdom_ind)
+ index_fields, particle_octree, pdom_ind, self.pf.geometry)
vals = op.finalize()
if vals is None: return
if isinstance(vals, list):
diff -r b699b4b430a41bacab986e82c472fa2efcac02a4 -r 0050516c82a9622e381ef5ecc64ac444d64f9bf1 yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -970,7 +970,7 @@
def load_particles(data, length_unit = None, bbox=None,
sim_time=0.0, mass_unit = None, time_unit = None,
velocity_unit=None, periodicity=(True, True, True),
- n_ref = 64, over_refine_factor = 1):
+ n_ref = 64, over_refine_factor = 1, geometry = "cartesian"):
r"""Load a set of particles into yt as a
:class:`~yt.frontends.stream.data_structures.StreamParticleHandler`.
@@ -1083,7 +1083,7 @@
handler.simulation_time = sim_time
handler.cosmology_simulation = 0
- spf = StreamParticlesDataset(handler)
+ spf = StreamParticlesDataset(handler, geometry = geometry)
spf.n_ref = n_ref
spf.over_refine_factor = over_refine_factor
diff -r b699b4b430a41bacab986e82c472fa2efcac02a4 -r 0050516c82a9622e381ef5ecc64ac444d64f9bf1 yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -95,7 +95,7 @@
if self.update_values == 1:
for j in range(nf):
field_pointers[j][i] = field_vals[j]
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
def process_grid(self, gobj,
@@ -424,3 +424,56 @@
return
deposit_mesh_id = MeshIdentifier
+
+cdef class NNParticleField(ParticleDepositOperation):
+ cdef np.float64_t *nnfield
+ cdef np.float64_t *distfield
+ cdef public object onnfield
+ cdef public object odistfield
+ def initialize(self):
+ self.onnfield = np.zeros(self.nvals, dtype="float64", order='F')
+ cdef np.ndarray arr = self.onnfield
+ self.nnfield = <np.float64_t*> arr.data
+
+ self.odistfield = np.zeros(self.nvals, dtype="float64", order='F')
+ self.odistfield[:] = np.inf
+ arr = self.odistfield
+ self.distfield = <np.float64_t*> arr.data
+
+ @cython.cdivision(True)
+ cdef void process(self, int dim[3],
+ np.float64_t left_edge[3],
+ np.float64_t dds[3],
+ np.int64_t offset,
+ np.float64_t ppos[3],
+ np.float64_t *fields,
+ np.int64_t domain_ind
+ ):
+ # This one is a bit slow. Every grid cell is going to be iterated
+ # over, and we're going to deposit particles in it.
+ cdef int ii[3], i, j, k
+ cdef np.int64_t ggind
+ cdef np.float64_t r2, gpos[3]
+ gpos[0] = left_edge[0] + 0.5 * dds[0]
+ for i in range(dim[0]):
+ gpos[1] = left_edge[1] + 0.5 * dds[1]
+ for j in range(dim[1]):
+ gpos[2] = left_edge[2] + 0.5 * dds[2]
+ for k in range(dim[2]):
+ ggind = gind(i, j, k, dim) + offset
+ r2 = ((ppos[0] - gpos[0])*(ppos[0] - gpos[0]) +
+ (ppos[1] - gpos[1])*(ppos[1] - gpos[1]) +
+ (ppos[2] - gpos[2])*(ppos[2] - gpos[2]))
+ if r2 < self.distfield[ggind]:
+ self.distfield[ggind] = r2
+ self.nnfield[ggind] = fields[0]
+ gpos[2] += dds[2]
+ gpos[1] += dds[1]
+ gpos[0] += dds[0]
+ return
+
+ def finalize(self):
+ return self.onnfield
+
+deposit_nearest = NNParticleField
+
diff -r b699b4b430a41bacab986e82c472fa2efcac02a4 -r 0050516c82a9622e381ef5ecc64ac444d64f9bf1 yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -36,27 +36,6 @@
np.int64_t pn # Particle number
np.float64_t r2 # radius**2
- at cython.cdivision(True)
- at cython.boundscheck(False)
- at cython.wraparound(False)
-cdef inline np.float64_t r2dist(np.float64_t ppos[3],
- np.float64_t cpos[3],
- np.float64_t DW[3],
- bint periodicity[3]):
- cdef int i
- cdef np.float64_t r2, DR
- r2 = 0.0
- for i in range(3):
- DR = (ppos[i] - cpos[i])
- if not periodicity[i]:
- pass
- elif (DR > DW[i]/2.0):
- DR -= DW[i]
- elif (DR < -DW[i]/2.0):
- DR += DW[i]
- r2 += DR * DR
- return r2
-
cdef class ParticleSmoothOperation:
# We assume each will allocate and define their own temporary storage
cdef public object nvals
@@ -71,6 +50,7 @@
cdef np.float64_t *ppos
# Note that we are preallocating here, so this is *not* threadsafe.
cdef NeighborList *neighbors
+ cdef void (*pos_setup)(np.float64_t ipos[3], np.float64_t opos[3])
cdef void neighbor_process(self, int dim[3], np.float64_t left_edge[3],
np.float64_t dds[3], np.float64_t *ppos,
np.float64_t **fields, np.int64_t nneighbors,
diff -r b699b4b430a41bacab986e82c472fa2efcac02a4 -r 0050516c82a9622e381ef5ecc64ac444d64f9bf1 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -18,7 +18,7 @@
import numpy as np
from libc.stdlib cimport malloc, free, realloc
cimport cython
-from libc.math cimport sqrt, fabs
+from libc.math cimport sqrt, fabs, sin, cos
from fp_utils cimport *
from oct_container cimport Oct, OctAllocationContainer, \
@@ -37,6 +37,40 @@
else:
return 1
+ at cython.cdivision(True)
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+cdef np.float64_t r2dist(np.float64_t ppos[3],
+ np.float64_t cpos[3],
+ np.float64_t DW[3],
+ bint periodicity[3],
+ np.float64_t max_dist2):
+ cdef int i
+ cdef np.float64_t r2, DR
+ r2 = 0.0
+ for i in range(3):
+ DR = (ppos[i] - cpos[i])
+ if not periodicity[i]:
+ pass
+ elif (DR > DW[i]/2.0):
+ DR -= DW[i]
+ elif (DR < -DW[i]/2.0):
+ DR += DW[i]
+ r2 += DR * DR
+ if max_dist2 >= 0.0 and r2 > max_dist2:
+ return -1.0
+ return r2
+
+cdef void spherical_coord_setup(np.float64_t ipos[3], np.float64_t opos[3]):
+ opos[0] = ipos[0] * sin(ipos[1]) * cos(ipos[2])
+ opos[1] = ipos[0] * sin(ipos[1]) * sin(ipos[2])
+ opos[2] = ipos[0] * cos(ipos[1])
+
+cdef void cart_coord_setup(np.float64_t ipos[3], np.float64_t opos[3]):
+ opos[0] = ipos[0]
+ opos[1] = ipos[1]
+ opos[2] = ipos[2]
+
cdef class ParticleSmoothOperation:
def __init__(self, nvals, nfields, max_neighbors):
# This is the set of cells, in grids, blocks or octs, we are handling.
@@ -66,7 +100,8 @@
periodicity = (True, True, True),
index_fields = None,
OctreeContainer particle_octree = None,
- np.ndarray[np.int64_t, ndim=1] pdom_ind = None):
+ np.ndarray[np.int64_t, ndim=1] pdom_ind = None,
+ geometry = "cartesian"):
# This will be a several-step operation.
#
# We first take all of our particles and assign them to Octs. If they
@@ -108,6 +143,25 @@
cdef np.ndarray[np.int64_t, ndim=2] doff_m
cdef np.ndarray[np.float64_t, ndim=1] tarr
cdef np.ndarray[np.float64_t, ndim=4] iarr
+ cdef np.ndarray[np.float64_t, ndim=2] cart_positions
+ if geometry == "cartesian":
+ self.pos_setup = cart_coord_setup
+ cart_positions = positions
+ elif geometry == "spherical":
+ self.pos_setup = spherical_coord_setup
+ cart_positions = np.empty((positions.shape[0], 3), dtype="float64")
+
+ cart_positions[:,0] = positions[:,0] * \
+ np.sin(positions[:,1]) * \
+ np.cos(positions[:,2])
+ cart_positions[:,1] = positions[:,0] * \
+ np.sin(positions[:,1]) * \
+ np.sin(positions[:,2])
+ cart_positions[:,2] = positions[:,0] * \
+ np.cos(positions[:,1])
+ periodicity = (False, False, False)
+ else:
+ raise NotImplementedError
dims[0] = dims[1] = dims[2] = (1 << mesh_octree.oref)
cdef int nz = dims[0] * dims[1] * dims[2]
numpart = positions.shape[0]
@@ -174,6 +228,7 @@
# Now doff is full of offsets to the first entry in the pind that
# refers to that oct's particles.
ppos = <np.float64_t *> positions.data
+ cart_pos = <np.float64_t *> cart_positions.data
doffs = <np.int64_t*> doff.data
pinds = <np.int64_t*> pind.data
pcounts = <np.int64_t*> pcount.data
@@ -213,7 +268,7 @@
free(neighbors)
nproc += 1
self.neighbor_process(dims, moi.left_edge, moi.dds,
- ppos, field_pointers, nneighbors, nind, doffs,
+ cart_pos, field_pointers, nneighbors, nind, doffs,
pinds, pcounts, offset, index_field_pointers)
#print "VISITED", visited.sum(), visited.size,
#print 100.0*float(visited.sum())/visited.size
@@ -253,7 +308,7 @@
if self.curn < self.maxn:
cur = &self.neighbors[self.curn]
cur.pn = pn
- cur.r2 = r2dist(ppos, cpos, self.DW, self.periodicity)
+ cur.r2 = r2dist(ppos, cpos, self.DW, self.periodicity, -1)
self.curn += 1
if self.curn == self.maxn:
# This time we sort it, so that future insertions will be able
@@ -262,7 +317,10 @@
Neighbor_compare)
return
# This will go (curn - 1) through 0.
- r2_c = r2dist(ppos, cpos, self.DW, self.periodicity)
+ r2_o = self.neighbors[self.curn - 1].r2
+ r2_c = r2dist(ppos, cpos, self.DW, self.periodicity, r2_o)
+ # Early terminate
+ if r2_c < 0: return
pn_c = pn
for i in range((self.curn - 1), -1, -1):
# First we evaluate against i. If our candidate radius is greater
@@ -320,15 +378,16 @@
# units supplied. We can now iterate over every cell in the block and
# every particle to find the nearest. We will use a priority heap.
cdef int i, j, k, ntot, nntot, m
- cdef np.float64_t cpos[3]
+ cdef np.float64_t cpos[3], opos[3]
cpos[0] = left_edge[0] + 0.5*dds[0]
for i in range(dim[0]):
cpos[1] = left_edge[1] + 0.5*dds[1]
for j in range(dim[1]):
cpos[2] = left_edge[2] + 0.5*dds[2]
for k in range(dim[2]):
+ self.pos_setup(cpos, opos)
self.neighbor_find(nneighbors, nind, doffs, pcounts,
- pinds, ppos, cpos)
+ pinds, ppos, opos)
# Now we have all our neighbors in our neighbor list.
if self.curn <-1*self.maxn:
ntot = nntot = 0
@@ -337,7 +396,7 @@
nntot += 1
ntot += pcounts[nind[m]]
print "SOMETHING WRONG", self.curn, nneighbors, ntot, nntot
- self.process(offset, i, j, k, dim, cpos, fields,
+ self.process(offset, i, j, k, dim, opos, fields,
index_fields)
cpos[2] += dds[2]
cpos[1] += dds[1]
@@ -407,3 +466,76 @@
return
volume_weighted_smooth = VolumeWeightedSmooth
+
+cdef class NearestNeighborSmooth(ParticleSmoothOperation):
+ cdef np.float64_t *fp
+ cdef public object vals
+ def initialize(self):
+ cdef np.ndarray tarr
+ assert(self.nfields == 1)
+ tarr = np.zeros(self.nvals, dtype="float64", order="F")
+ self.vals = tarr
+ self.fp = <np.float64_t *> tarr.data
+
+ def finalize(self):
+ return self.vals
+
+ @cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ cdef void process(self, np.int64_t offset, int i, int j, int k,
+ int dim[3], np.float64_t cpos[3], np.float64_t **fields,
+ np.float64_t **index_fields):
+ # We have our i, j, k for our cell, as well as the cell position.
+ # We also have a list of neighboring particles with particle numbers.
+ cdef np.int64_t pn
+ # We get back our mass
+ # rho_i = sum(j = 1 .. n) m_j * W_ij
+ pn = self.neighbors[0].pn
+ self.fp[gind(i,j,k,dim) + offset] = fields[0][pn]
+ #self.fp[gind(i,j,k,dim) + offset] = self.neighbors[0].r2
+ return
+
+nearest_smooth = NearestNeighborSmooth
+
+cdef class IDWInterpolationSmooth(ParticleSmoothOperation):
+ cdef np.float64_t *fp
+ cdef public int p2
+ cdef public object vals
+ def initialize(self):
+ cdef np.ndarray tarr
+ assert(self.nfields == 1)
+ tarr = np.zeros(self.nvals, dtype="float64", order="F")
+ self.vals = tarr
+ self.fp = <np.float64_t *> tarr.data
+ self.p2 = 2 # Power, for IDW, in units of 2. So we only do even p's.
+
+ def finalize(self):
+ return self.vals
+
+ @cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ cdef void process(self, np.int64_t offset, int i, int j, int k,
+ int dim[3], np.float64_t cpos[3], np.float64_t **fields,
+ np.float64_t **index_fields):
+ # We have our i, j, k for our cell, as well as the cell position.
+ # We also have a list of neighboring particles with particle numbers.
+ cdef np.int64_t pn, ni, di
+ cdef np.float64_t total_weight = 0.0, total_value = 0.0, r2, val, w
+ # We're going to do a very simple IDW average
+ if self.neighbors[0].r2 == 0.0:
+ pn = self.neighbors[0].pn
+ self.fp[gind(i,j,k,dim) + offset] = fields[0][pn]
+ for ni in range(self.curn):
+ r2 = self.neighbors[ni].r2
+ val = fields[0][self.neighbors[ni].pn]
+ w = r2
+ for di in range(self.p2 - 1):
+ w *= r2
+ total_value += w * val
+ total_weight += w
+ self.fp[gind(i,j,k,dim) + offset] = total_value / total_weight
+ return
+
+idw_smooth = IDWInterpolationSmooth
diff -r b699b4b430a41bacab986e82c472fa2efcac02a4 -r 0050516c82a9622e381ef5ecc64ac444d64f9bf1 yt/utilities/lib/ContourFinding.pyx
--- a/yt/utilities/lib/ContourFinding.pyx
+++ b/yt/utilities/lib/ContourFinding.pyx
@@ -24,7 +24,6 @@
OctreeContainer, OctInfo
from yt.geometry.oct_visitors cimport \
Oct
-from yt.geometry.particle_smooth cimport r2dist
from .amr_kdtools cimport _find_node, Node
from .grid_traversal cimport VolumeContainer, PartitionedGrid, \
vc_index, vc_pos_index
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