[yt-svn] commit/yt: ngoldbaum: Merged in MatthewTurk/yt (pull request #1193)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Sep 16 20:36:52 PDT 2014
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/b3c9e119a5b6/
Changeset: b3c9e119a5b6
Branch: yt
User: ngoldbaum
Date: 2014-09-17 03:36:33+00:00
Summary: Merged in MatthewTurk/yt (pull request #1193)
Fix bug with mismatched smoothing trees.
Affected #: 4 files
diff -r c7ff1e9f725774c2a2cdc757755ada9ec752386a -r b3c9e119a5b61c46aad9467107d27dc0fb6b5b50 yt/data_objects/octree_subset.py
--- a/yt/data_objects/octree_subset.py
+++ b/yt/data_objects/octree_subset.py
@@ -177,7 +177,7 @@
self.ds.domain_left_edge,
self.ds.domain_right_edge,
over_refine = self._oref)
- particle_octree.n_ref = nneighbors / 2
+ particle_octree.n_ref = nneighbors
particle_octree.add(morton)
particle_octree.finalize()
pdom_ind = particle_octree.domain_ind(self.selector)
diff -r c7ff1e9f725774c2a2cdc757755ada9ec752386a -r b3c9e119a5b61c46aad9467107d27dc0fb6b5b50 yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -281,7 +281,7 @@
dds[i] = (self.DRE[i] - self.DLE[i])/self.nn[i]
ind[i] = <np.int64_t> (floor((ppos[i] - self.DLE[i])/dds[i]))
cp[i] = (ind[i] + 0.5) * dds[i] + self.DLE[i]
- ipos[i] = 0
+ ipos[i] = 0 # We add this to ind later, so it should be zero.
ind32[i] = ind[i]
self.get_root(ind32, &next)
# We want to stop recursing when there's nowhere else to go
@@ -307,12 +307,6 @@
cdef np.float64_t factor = 1.0 / (1 << (self.oref-1))
if self.oref == 0: factor = 2.0
for i in range(3):
- # This will happen *after* we quit out, so we need to back out the
- # last change to cp
- if ind[i] == 1:
- cp[i] -= dds[i]/2.0 # Now centered
- else:
- cp[i] += dds[i]/2.0
# We don't normally need to change dds[i] as it has been halved
# from the oct width, thus making it already the cell width.
# But, since not everything has the cell width equal to have the
@@ -926,7 +920,7 @@
# to zero; i.e., the number of dimensions along which we are aligned.
# For now, we assume we will not be doing this along all three zeros,
# because that would be pretty tricky.
- if i == j == k == 0: return olist
+ if i == j == k == 1: return olist
cdef np.int64_t n[3], ind[3], off[3][2], ii, ij, ik, ci
ind[0] = 1 - i
ind[1] = 1 - j
diff -r c7ff1e9f725774c2a2cdc757755ada9ec752386a -r b3c9e119a5b61c46aad9467107d27dc0fb6b5b50 yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -53,10 +53,15 @@
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,
- np.int64_t *nind, np.int64_t *doffs,
+ np.float64_t **fields,
+ np.int64_t *doffs, np.int64_t **nind,
np.int64_t *pinds, np.int64_t *pcounts,
- np.int64_t offset, np.float64_t **index_fields)
+ np.int64_t offset, np.float64_t **index_fields,
+ OctreeContainer octree, np.int64_t domain_id,
+ int *nsize)
+ cdef int neighbor_search(self, np.float64_t pos[3], OctreeContainer octree,
+ np.int64_t **nind, int *nsize,
+ np.int64_t nneighbors, np.int64_t domain_id, Oct **oct = ?)
cdef void neighbor_eval(self, np.int64_t pn, np.float64_t ppos[3],
np.float64_t cpos[3])
cdef void neighbor_reset(self)
diff -r c7ff1e9f725774c2a2cdc757755ada9ec752386a -r b3c9e119a5b61c46aad9467107d27dc0fb6b5b50 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -135,8 +135,8 @@
cdef int nsize = 0
cdef np.int64_t *nind = NULL
cdef OctInfo moi, poi
- cdef Oct *oct, **neighbors = NULL
- cdef np.int64_t nneighbors, numpart, offset, local_ind
+ cdef Oct *oct
+ cdef np.int64_t numpart, offset, local_ind
cdef np.int64_t moff_p, moff_m
cdef np.int64_t *doffs, *pinds, *pcounts, poff
cdef np.ndarray[np.int64_t, ndim=1] pind, doff, pdoms, pcount
@@ -232,11 +232,8 @@
doffs = <np.int64_t*> doff.data
pinds = <np.int64_t*> pind.data
pcounts = <np.int64_t*> pcount.data
- nsize = 27
- nind = <np.int64_t *> malloc(sizeof(np.int64_t)*nsize)
cdef np.ndarray[np.uint8_t, ndim=1] visited
visited = np.zeros(mdom_ind.shape[0], dtype="uint8")
- cdef int maxnei = 0
cdef int nproc = 0
for i in range(oct_positions.shape[0]):
for j in range(3):
@@ -246,34 +243,49 @@
if visited[oct.domain_ind - moff_m] == 1: continue
visited[oct.domain_ind - moff_m] = 1
if offset < 0: continue
- # These will be PARTICLE octree neighbors.
- oct = particle_octree.get(pos, &poi)
- neighbors = particle_octree.neighbors(&poi, &nneighbors, oct)
- if nneighbors > maxnei:
- maxnei = nneighbors
- # Now we have all our neighbors. And, we should be set for what
- # else we need to do.
- if nneighbors > nsize:
- nind = <np.int64_t *> realloc(
- nind, sizeof(np.int64_t)*nneighbors)
- nsize = nneighbors
- for j in range(nneighbors):
- # Particle octree neighbor indices
- nind[j] = neighbors[j].domain_ind - moff_p
- for n in range(j):
- if nind[j] == nind[n]:
- nind[j] = -1
- break
- # This is allocated by the neighbors function, so we deallocate it.
- free(neighbors)
nproc += 1
self.neighbor_process(dims, moi.left_edge, moi.dds,
- cart_pos, field_pointers, nneighbors, nind, doffs,
- pinds, pcounts, offset, index_field_pointers)
+ cart_pos, field_pointers, doffs, &nind,
+ pinds, pcounts, offset, index_field_pointers,
+ particle_octree, domain_id, &nsize)
#print "VISITED", visited.sum(), visited.size,
#print 100.0*float(visited.sum())/visited.size
if nind != NULL:
free(nind)
+
+ cdef int neighbor_search(self, np.float64_t pos[3], OctreeContainer octree,
+ np.int64_t **nind, int *nsize,
+ np.int64_t nneighbors, np.int64_t domain_id,
+ Oct **oct = NULL):
+ cdef OctInfo oi
+ cdef Oct *ooct, **neighbors
+ cdef int j
+ cdef np.int64_t moff = octree.get_domain_offset(domain_id)
+ ooct = octree.get(pos, &oi)
+ if oct != NULL and ooct == oct[0]:
+ return nneighbors
+ oct[0] = ooct
+ if nind[0] == NULL:
+ nsize[0] = 27
+ nind[0] = <np.int64_t *> malloc(sizeof(np.int64_t)*nsize[0])
+ neighbors = octree.neighbors(&oi, &nneighbors, ooct)
+ # Now we have all our neighbors. And, we should be set for what
+ # else we need to do.
+ if nneighbors > nsize[0]:
+ nind[0] = <np.int64_t *> realloc(
+ nind[0], sizeof(np.int64_t)*nneighbors)
+ nsize[0] = nneighbors
+
+ for j in range(nneighbors):
+ # Particle octree neighbor indices
+ nind[0][j] = neighbors[j].domain_ind - moff
+ for n in range(j):
+ if nind[0][j] == nind[0][n]:
+ nind[0][j] = -1
+ break
+ # This is allocated by the neighbors function, so we deallocate it.
+ free(neighbors)
+ return nneighbors
@cython.cdivision(True)
@cython.boundscheck(False)
@@ -369,16 +381,19 @@
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,
- np.int64_t *nind, np.int64_t *doffs,
+ np.float64_t **fields,
+ np.int64_t *doffs, np.int64_t **nind,
np.int64_t *pinds, np.int64_t *pcounts,
np.int64_t offset,
- np.float64_t **index_fields):
+ np.float64_t **index_fields,
+ OctreeContainer octree, np.int64_t domain_id,
+ int *nsize):
# Note that we assume that fields[0] == smoothing length in the native
# 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 int i, j, k, ntot, nntot, m, nneighbors
cdef np.float64_t cpos[3], opos[3]
+ cdef Oct* oct = NULL
cpos[0] = left_edge[0] + 0.5*dds[0]
for i in range(dim[0]):
cpos[1] = left_edge[1] + 0.5*dds[1]
@@ -386,15 +401,17 @@
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,
+ nneighbors = self.neighbor_search(opos, octree,
+ nind, nsize, nneighbors, domain_id, &oct)
+ self.neighbor_find(nneighbors, nind[0], doffs, pcounts,
pinds, ppos, opos)
# Now we have all our neighbors in our neighbor list.
if self.curn <-1*self.maxn:
ntot = nntot = 0
for m in range(nneighbors):
- if nind[m] < 0: continue
+ if nind[0][m] < 0: continue
nntot += 1
- ntot += pcounts[nind[m]]
+ ntot += pcounts[nind[0][m]]
print "SOMETHING WRONG", self.curn, nneighbors, ntot, nntot
self.process(offset, i, j, k, dim, opos, fields,
index_fields)
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