[yt-svn] commit/yt: 5 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Sep 16 20:36:39 PDT 2014
5 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/f8fc80f9afb3/
Changeset: f8fc80f9afb3
Branch: yt
User: MatthewTurk
Date: 2014-09-09 10:46:41+00:00
Summary: Fix bug with mismatched smoothing trees.
This fixes a somewhat pernicious bug. It usually showed up with the
over_refine_factor for a particle dataset was greater than one, although in
principle it could also affect when it was equal to one. What can happen is
that when iterating over mesh items in the mesh octree, sometimes we would walk
outside the particle octree oct that we had identified for the neighbor
searching. This fixes the issue by adding in a walk over the particle search
octree during the course of the iteration over mesh elements. This could in
principle lead to considerably more runtime, but I have added a simple early
termination process to try to ease that.
Affected #: 4 files
diff -r c8042b5d05d064f8e5c38f03dd58a680a52036f3 -r f8fc80f9afb3cfd097034555889a64cd609fdc85 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 c8042b5d05d064f8e5c38f03dd58a680a52036f3 -r f8fc80f9afb3cfd097034555889a64cd609fdc85 yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -418,7 +418,7 @@
cdef np.ndarray[np.uint8_t, ndim=1] coords
cdef OctVisitorData data
self.setup_data(&data, domain_id)
- coords = np.zeros((num_cells*8), dtype="uint8")
+ coords = np.zeros((num_cells*data.nz), dtype="uint8")
data.array = <void *> coords.data
self.visit_all_octs(selector, oct_visitors.mask_octs, &data)
return coords.astype("bool")
diff -r c8042b5d05d064f8e5c38f03dd58a680a52036f3 -r f8fc80f9afb3cfd097034555889a64cd609fdc85 yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -53,10 +53,14 @@
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)
+ 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 c8042b5d05d064f8e5c38f03dd58a680a52036f3 -r f8fc80f9afb3cfd097034555889a64cd609fdc85 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):
@@ -247,33 +244,49 @@
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)
#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 +382,18 @@
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):
# 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, nsize
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,6 +401,8 @@
cpos[2] = left_edge[2] + 0.5*dds[2]
for k in range(dim[2]):
self.pos_setup(cpos, opos)
+ nneighbors = self.neighbor_search(opos, octree,
+ &nind, &nsize, nneighbors, domain_id, &oct)
self.neighbor_find(nneighbors, nind, doffs, pcounts,
pinds, ppos, opos)
# Now we have all our neighbors in our neighbor list.
https://bitbucket.org/yt_analysis/yt/commits/d23a479f5ebe/
Changeset: d23a479f5ebe
Branch: yt
User: MatthewTurk
Date: 2014-09-09 17:04:49+00:00
Summary: These are offset by one.
Affected #: 1 file
diff -r f8fc80f9afb3cfd097034555889a64cd609fdc85 -r d23a479f5ebe66124b50d57a87afed0f19df407a yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -926,7 +926,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
https://bitbucket.org/yt_analysis/yt/commits/dbcd5e67e513/
Changeset: dbcd5e67e513
Branch: yt
User: MatthewTurk
Date: 2014-09-09 17:09:01+00:00
Summary: Minor cleanups, reducing allocations.
Affected #: 3 files
diff -r d23a479f5ebe66124b50d57a87afed0f19df407a -r dbcd5e67e513bd44e83e7bbb37a212fb91edfbe7 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] = ind[i]
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
diff -r d23a479f5ebe66124b50d57a87afed0f19df407a -r dbcd5e67e513bd44e83e7bbb37a212fb91edfbe7 yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -54,10 +54,11 @@
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 *doffs, np.int64_t *nind,
+ 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,
- OctreeContainer octree, np.int64_t domain_id)
+ 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 = ?)
diff -r d23a479f5ebe66124b50d57a87afed0f19df407a -r dbcd5e67e513bd44e83e7bbb37a212fb91edfbe7 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -243,12 +243,11 @@
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.
nproc += 1
self.neighbor_process(dims, moi.left_edge, moi.dds,
- cart_pos, field_pointers, doffs, nind,
+ cart_pos, field_pointers, doffs, &nind,
pinds, pcounts, offset, index_field_pointers,
- particle_octree, domain_id)
+ particle_octree, domain_id, &nsize)
#print "VISITED", visited.sum(), visited.size,
#print 100.0*float(visited.sum())/visited.size
if nind != NULL:
@@ -383,15 +382,16 @@
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 *doffs, np.int64_t *nind,
+ 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,
- OctreeContainer octree, np.int64_t domain_id):
+ 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, nneighbors, nsize
+ 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]
@@ -402,16 +402,16 @@
for k in range(dim[2]):
self.pos_setup(cpos, opos)
nneighbors = self.neighbor_search(opos, octree,
- &nind, &nsize, nneighbors, domain_id, &oct)
- self.neighbor_find(nneighbors, nind, doffs, pcounts,
+ 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)
https://bitbucket.org/yt_analysis/yt/commits/9ac61672fc1d/
Changeset: 9ac61672fc1d
Branch: yt
User: MatthewTurk
Date: 2014-09-16 14:38:40+00:00
Summary: Fixing segfault for ARTIO and RAMSES.
Affected #: 1 file
diff -r dbcd5e67e513bd44e83e7bbb37a212fb91edfbe7 -r 9ac61672fc1d0614e8b3e1e98764d3b062f16912 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] = ind[i]
+ 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
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