[yt-svn] commit/yt: 18 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Sep 15 13:21:28 PDT 2015
18 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/58ca9d4b7c4e/
Changeset: 58ca9d4b7c4e
Branch: yt
User: MatthewTurk
Date: 2015-08-24 01:30:32+00:00
Summary: Returning a closure for smoothing length.
Affected #: 1 file
diff -r 0d479d8c21a4557c8fc17329c0171ff6dec4e8e8 -r 58ca9d4b7c4e87f14025c5d14e27e79aea64c4ce yt/frontends/tipsy/fields.py
--- a/yt/frontends/tipsy/fields.py
+++ b/yt/frontends/tipsy/fields.py
@@ -60,15 +60,19 @@
def setup_gas_particle_fields(self, ptype):
- def _smoothing_length(field, data):
- # For now, we hardcode num_neighbors. We should make this configurable
- # in the future.
- num_neighbors = 64
- fn, = add_nearest_neighbor_field(ptype, "particle_position", self, num_neighbors)
- return data[ptype, 'nearest_neighbor_distance_%d' % num_neighbors]
+ num_neighbors = 65
+ fn, = add_nearest_neighbor_field(ptype, "particle_position", self, num_neighbors)
+ def _func():
+ def _smoothing_length(field, data):
+ # For now, we hardcode num_neighbors. We should make this configurable
+ # in the future.
+ rv = data[ptype, 'nearest_neighbor_distance_%d' % num_neighbors]
+ #np.maximum(rv, 0.5*data[ptype, "Epsilon"], rv)
+ return rv
+ return _smoothing_length
self.add_field(
(ptype, "smoothing_length"),
- function=_smoothing_length,
+ function=_func(),
particle_type=True,
units="code_length")
https://bitbucket.org/yt_analysis/yt/commits/9cbff01341c1/
Changeset: 9cbff01341c1
Branch: yt
User: MatthewTurk
Date: 2015-08-25 21:31:31+00:00
Summary: Adding neighbors-of-neighbors to smoothing.
Affected #: 3 files
diff -r 58ca9d4b7c4e87f14025c5d14e27e79aea64c4ce -r 9cbff01341c1b0efa3ab089659d488938cb0718d yt/frontends/tipsy/fields.py
--- a/yt/frontends/tipsy/fields.py
+++ b/yt/frontends/tipsy/fields.py
@@ -38,7 +38,8 @@
'FeMassFrac':("FeMassFrac", ("dimensionless", ["Fe_fraction"], None)),
'c':("c", ("code_velocity", [""], None)),
'acc':("acc", ("code_velocity / code_time", [""], None)),
- 'accg':("accg", ("code_velocity / code_time", [""], None))}
+ 'accg':("accg", ("code_velocity / code_time", [""], None)),
+ 'smoothlength':('smoothlength', ("code_length", ["smoothing_length"], None))}
def __init__(self, ds, field_list, slice_info = None):
for field in field_list:
diff -r 58ca9d4b7c4e87f14025c5d14e27e79aea64c4ce -r 9cbff01341c1b0efa3ab089659d488938cb0718d yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -58,7 +58,8 @@
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 = ?)
+ np.int64_t nneighbors, np.int64_t domain_id,
+ Oct **oct = ?, int extra_layer = ?)
cdef void neighbor_process_particle(self, np.float64_t cpos[3],
np.float64_t *ppos,
np.float64_t **fields,
diff -r 58ca9d4b7c4e87f14025c5d14e27e79aea64c4ce -r 9cbff01341c1b0efa3ab089659d488938cb0718d yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -394,11 +394,12 @@
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):
+ Oct **oct = NULL, int extra_layer = 0):
cdef OctInfo oi
cdef Oct *ooct
- cdef Oct **neighbors
- cdef int j
+ cdef Oct **neighbors, **first_layer
+ cdef int j, total_neighbors = 0, initial_layer = 0
+ cdef int layer_ind = 0
cdef np.int64_t moff = octree.get_domain_offset(domain_id)
ooct = octree.get(pos, &oi)
if oct != NULL and ooct == oct[0]:
@@ -407,24 +408,53 @@
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, self.periodicity)
- # 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
+ # This is our "seed" set of neighbors. If we are asked to, we will
+ # create a master list of neighbors that is much bigger and includes
+ # everything.
+ layer_ind = 0
+ first_layer = NULL
+ while 1:
+ neighbors = octree.neighbors(&oi, &nneighbors, ooct, self.periodicity)
+ # Now we have all our neighbors. And, we should be set for what
+ # else we need to do.
+ if total_neighbors + nneighbors > nsize[0]:
+ nind[0] = <np.int64_t *> realloc(
+ nind[0], sizeof(np.int64_t)*(nneighbors + total_neighbors))
+ nsize[0] = nneighbors + total_neighbors
+ for j in range(nneighbors):
+ # Particle octree neighbor indices
+ nind[0][j + total_neighbors] = neighbors[j].domain_ind - moff
+ total_neighbors += nneighbors
+ if extra_layer == 0:
+ # Not adding on any additional layers here.
+ free(neighbors)
+ neighbors = NULL
+ break
+ if initial_layer == 0:
+ initial_layer = nneighbors
+ first_layer = neighbors
+ else:
+ # Allocated internally; we free this in the loops if we aren't
+ # tracking it
+ free(neighbors)
+ neighbors = NULL
+ ooct = first_layer[layer_ind]
+ layer_ind += 1
+ if layer_ind == initial_layer:
+ neighbors
+ break
+
- for j in range(nneighbors):
+ for j in range(total_neighbors):
# 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
+ if first_layer != NULL:
+ free(first_layer)
+ return total_neighbors
@cython.cdivision(True)
@cython.boundscheck(False)
@@ -572,7 +602,7 @@
cdef np.float64_t opos[3]
self.pos_setup(cpos, opos)
nneighbors = self.neighbor_search(opos, octree,
- nind, nsize, nneighbors, domain_id, &oct)
+ nind, nsize, nneighbors, domain_id, &oct, 1)
self.neighbor_find(nneighbors, nind[0], doffs, pcounts, pinds, ppos, opos)
self.process(offset, i, j, k, dim, opos, fields, index_fields)
https://bitbucket.org/yt_analysis/yt/commits/b36a9df1bc87/
Changeset: b36a9df1bc87
Branch: yt
User: MatthewTurk
Date: 2015-08-30 02:27:04+00:00
Summary: Fixing extra-layer neighbor calculation, and make smoothing use it.
Affected #: 1 file
diff -r 9cbff01341c1b0efa3ab089659d488938cb0718d -r b36a9df1bc8713dbf080d1a080ccfee8ca04f3aa yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -447,10 +447,10 @@
for j in range(total_neighbors):
# Particle octree neighbor indices
+ if nind[0][j] == -1: continue
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.
if first_layer != NULL:
free(first_layer)
@@ -564,7 +564,7 @@
for k in range(dim[2]):
self.pos_setup(cpos, opos)
nneighbors = self.neighbor_search(opos, octree,
- nind, nsize, nneighbors, domain_id, &oct)
+ nind, nsize, nneighbors, domain_id, &oct, 1)
self.neighbor_find(nneighbors, nind[0], doffs, pcounts,
pinds, ppos, opos)
# Now we have all our neighbors in our neighbor list.
https://bitbucket.org/yt_analysis/yt/commits/363e2f859540/
Changeset: 363e2f859540
Branch: yt
User: MatthewTurk
Date: 2015-08-30 02:36:26+00:00
Summary: Fix the smoothing length application for SPH
Affected #: 1 file
diff -r b36a9df1bc8713dbf080d1a080ccfee8ca04f3aa -r 363e2f859540ae7ad515d65a60ee6479ddac972d yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -642,6 +642,7 @@
# We also have a list of neighboring particles with particle numbers.
cdef int n, fi
cdef np.float64_t weight, r2, val, hsml, dens, mass, coeff, max_r
+ cdef np.float64_t max_hsml, ihsml, ihsml3
coeff = 0.0
cdef np.int64_t pn
# We get back our mass
@@ -658,10 +659,13 @@
if hsml < 0:
hsml = max_r
if hsml == 0: continue
+ ihsml = 1.0/hsml
+ ihsml3 = ihsml*ihsml*ihsml
# Usually this density has been computed
dens = fields[2][pn]
if dens == 0.0: continue
- weight = mass * self.sph_kernel(sqrt(r2) / hsml) / dens
+ weight = mass * self.sph_kernel(sqrt(r2) * ihsml) / dens
+ weight *= ihsml3
# Mass of the particle times the value
for fi in range(self.nfields - 3):
val = fields[fi + 3][pn]
https://bitbucket.org/yt_analysis/yt/commits/fc4a2ed4f537/
Changeset: fc4a2ed4f537
Branch: yt
User: MatthewTurk
Date: 2015-09-02 22:38:56+00:00
Summary: Apply constant to cubic spline, add in hooks for normalization.
Affected #: 3 files
diff -r 363e2f859540ae7ad515d65a60ee6479ddac972d -r fc4a2ed4f537ceacb2ae9236497f244e3af5f091 yt/fields/particle_fields.py
--- a/yt/fields/particle_fields.py
+++ b/yt/fields/particle_fields.py
@@ -790,13 +790,19 @@
kwargs = {}
if nneighbors:
kwargs['nneighbors'] = nneighbors
+ # This is for applying cutoffs, similar to in the SPLASH paper.
+ smooth_cutoff = data["index","cell_volume"]**(1./3)
+ smooth_cutoff.convert_to_units("code_length")
# volume_weighted smooth operations return lists of length 1.
rv = data.smooth(pos, [mass, hsml, dens, quan],
method="volume_weighted",
create_octree=True,
+ index_fields=[smooth_cutoff],
kernel_name=kernel_name)[0]
rv[np.isnan(rv)] = 0.0
# Now some quick unit conversions.
+ # This should be used when seeking a non-normalized value:
+ #rv /= hsml.uq**3 / hsml.uq.in_cgs().uq**3
rv = data.apply_units(rv, field_units)
return rv
registry.add_field(field_name, function = _vol_weight,
diff -r 363e2f859540ae7ad515d65a60ee6479ddac972d -r fc4a2ed4f537ceacb2ae9236497f244e3af5f091 yt/geometry/particle_deposit.pxd
--- a/yt/geometry/particle_deposit.pxd
+++ b/yt/geometry/particle_deposit.pxd
@@ -40,13 +40,14 @@
cdef inline np.float64_t sph_kernel_cubic(np.float64_t x) nogil:
cdef np.float64_t kernel
+ cdef np.float64_t C = 2.5464790894703255
if x <= 0.5:
kernel = 1.-6.*x*x*(1.-x)
elif x>0.5 and x<=1.0:
kernel = 2.*(1.-x)*(1.-x)*(1.-x)
else:
kernel = 0.
- return kernel
+ return kernel * C
########################################################
# Alternative SPH kernels for use with the Grid method #
diff -r 363e2f859540ae7ad515d65a60ee6479ddac972d -r fc4a2ed4f537ceacb2ae9236497f244e3af5f091 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -642,12 +642,13 @@
# We also have a list of neighboring particles with particle numbers.
cdef int n, fi
cdef np.float64_t weight, r2, val, hsml, dens, mass, coeff, max_r
- cdef np.float64_t max_hsml, ihsml, ihsml3
+ cdef np.float64_t max_hsml, ihsml, ihsml3, kern
coeff = 0.0
cdef np.int64_t pn
# We get back our mass
# rho_i = sum(j = 1 .. n) m_j * W_ij
max_r = sqrt(self.neighbors[self.curn-1].r2)
+ max_hsml = index_fields[0][gind(i,j,k,dim) + offset]
for n in range(self.curn):
# No normalization for the moment.
# fields[0] is the smoothing length.
@@ -656,21 +657,25 @@
# Smoothing kernel weight function
mass = fields[0][pn]
hsml = fields[1][pn]
+ dens = fields[2][pn]
if hsml < 0:
hsml = max_r
if hsml == 0: continue
ihsml = 1.0/hsml
+ hsml = fmax(max_hsml/2.0, hsml)
ihsml3 = ihsml*ihsml*ihsml
# Usually this density has been computed
- dens = fields[2][pn]
if dens == 0.0: continue
- weight = mass * self.sph_kernel(sqrt(r2) * ihsml) / dens
- weight *= ihsml3
+ weight = (mass / dens) * ihsml3
+ kern = self.sph_kernel(sqrt(r2) * ihsml)
+ weight *= kern
# Mass of the particle times the value
for fi in range(self.nfields - 3):
val = fields[fi + 3][pn]
self.fp[fi][gind(i,j,k,dim) + offset] += val * weight
self.fp[self.nfields - 3][gind(i,j,k,dim) + offset] += weight
+ # When looking for non-normalized values, uncomment:
+ #self.fp[self.nfields - 3][gind(i,j,k,dim) + offset] = 1.0
return
volume_weighted_smooth = VolumeWeightedSmooth
https://bitbucket.org/yt_analysis/yt/commits/6a0faa0d0fc8/
Changeset: 6a0faa0d0fc8
Branch: yt
User: ngoldbaum
Date: 2015-09-04 16:12:45+00:00
Summary: Uncommenting lines that remove normalization from SPH smoothing
Affected #: 2 files
diff -r fc4a2ed4f537ceacb2ae9236497f244e3af5f091 -r 6a0faa0d0fc85b1d97b118e4b2affd20e105242f yt/fields/particle_fields.py
--- a/yt/fields/particle_fields.py
+++ b/yt/fields/particle_fields.py
@@ -802,7 +802,7 @@
rv[np.isnan(rv)] = 0.0
# Now some quick unit conversions.
# This should be used when seeking a non-normalized value:
- #rv /= hsml.uq**3 / hsml.uq.in_cgs().uq**3
+ rv /= hsml.uq**3 / hsml.uq.in_cgs().uq**3
rv = data.apply_units(rv, field_units)
return rv
registry.add_field(field_name, function = _vol_weight,
diff -r fc4a2ed4f537ceacb2ae9236497f244e3af5f091 -r 6a0faa0d0fc85b1d97b118e4b2affd20e105242f yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -675,7 +675,7 @@
self.fp[fi][gind(i,j,k,dim) + offset] += val * weight
self.fp[self.nfields - 3][gind(i,j,k,dim) + offset] += weight
# When looking for non-normalized values, uncomment:
- #self.fp[self.nfields - 3][gind(i,j,k,dim) + offset] = 1.0
+ self.fp[self.nfields - 3][gind(i,j,k,dim) + offset] = 1.0
return
volume_weighted_smooth = VolumeWeightedSmooth
https://bitbucket.org/yt_analysis/yt/commits/c0ad27544548/
Changeset: c0ad27544548
Branch: yt
User: ngoldbaum
Date: 2015-09-04 17:02:00+00:00
Summary: Add (still broken) early termination logic to neighbor finding
Affected #: 2 files
diff -r 6a0faa0d0fc85b1d97b118e4b2affd20e105242f -r c0ad27544548666e24bb2e951b6838a62cbf121e yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -79,7 +79,8 @@
np.int64_t *pcounts,
np.int64_t *pinds,
np.float64_t *ppos,
- np.float64_t cpos[3])
+ np.float64_t cpos[3],
+ np.float64_t *dds)
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)
diff -r 6a0faa0d0fc85b1d97b118e4b2affd20e105242f -r c0ad27544548666e24bb2e951b6838a62cbf121e yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -1,3 +1,4 @@
+# cython: profile=True
"""
Particle smoothing in cells
@@ -522,14 +523,26 @@
np.int64_t *pcounts,
np.int64_t *pinds,
np.float64_t *ppos,
- np.float64_t cpos[3]
+ np.float64_t cpos[3],
+ np.float64_t *dds,
):
# We are now given the number of neighbors, the indices into the
# domains for them, and the number of particles for each.
cdef int ni, i, j
cdef np.int64_t offset, pn, pc
- cdef np.float64_t pos[3]
+ cdef np.float64_t pos[3], lpos[3], rpos[3], r2_trunc, lr2, rr2
self.neighbor_reset()
+ # terminate early if left and right edge r2 for this oct
+ # is bigger than final r2
+ if dds != NULL and self.curn == self.maxn:
+ for j in range(3):
+ lpos[j] = cpos[j] - dds[j]
+ rpos[j] = cpos[j] + dds[j]
+ r2_trunc = self.neighbors[self.curn - 1].r2
+ lr2 = r2dist(ppos, lpos, self.DW, self.periodicity, r2_trunc)
+ rr2 = r2dist(ppos, rpos, self.DW, self.periodicity, r2_trunc)
+ if lr2 == -1 and rr2 == -1:
+ return
for ni in range(nneighbors):
if nind[ni] == -1: continue
offset = doffs[nind[ni]]
@@ -566,7 +579,7 @@
nneighbors = self.neighbor_search(opos, octree,
nind, nsize, nneighbors, domain_id, &oct, 1)
self.neighbor_find(nneighbors, nind[0], doffs, pcounts,
- pinds, ppos, opos)
+ pinds, ppos, opos, dds)
# Now we have all our neighbors in our neighbor list.
if self.curn <-1*self.maxn:
ntot = nntot = 0
@@ -603,7 +616,8 @@
self.pos_setup(cpos, opos)
nneighbors = self.neighbor_search(opos, octree,
nind, nsize, nneighbors, domain_id, &oct, 1)
- self.neighbor_find(nneighbors, nind[0], doffs, pcounts, pinds, ppos, opos)
+ self.neighbor_find(nneighbors, nind[0], doffs, pcounts, pinds, ppos,
+ opos, NULL)
self.process(offset, i, j, k, dim, opos, fields, index_fields)
cdef class VolumeWeightedSmooth(ParticleSmoothOperation):
https://bitbucket.org/yt_analysis/yt/commits/2aec0c5340dd/
Changeset: 2aec0c5340dd
Branch: yt
User: ngoldbaum
Date: 2015-09-04 16:04:05+00:00
Summary: Create cache of oct edges to send to neighbor_find. Still doesn't work.
Affected #: 2 files
diff -r c0ad27544548666e24bb2e951b6838a62cbf121e -r 2aec0c5340ddca8513a345313b252bad1388d96a yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -55,7 +55,8 @@
np.int64_t *pinds, np.int64_t *pcounts,
np.int64_t offset, np.float64_t **index_fields,
OctreeContainer octree, np.int64_t domain_id,
- int *nsize)
+ int *nsize, np.float64_t *oct_left_edges,
+ np.float64_t *oct_dds)
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,
@@ -80,7 +81,8 @@
np.int64_t *pinds,
np.float64_t *ppos,
np.float64_t cpos[3],
- np.float64_t *dds)
+ np.float64_t* oct_left_edges,
+ np.float64_t* oct_dds)
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)
diff -r c0ad27544548666e24bb2e951b6838a62cbf121e -r 2aec0c5340ddca8513a345313b252bad1388d96a yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -155,6 +155,8 @@
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
+ cdef np.ndarray[np.float64_t, ndim=2] oct_left_edges, oct_dds
+ cdef OctInfo oinfo
if geometry == "cartesian":
self.pos_setup = cart_coord_setup
cart_positions = positions
@@ -178,6 +180,8 @@
numpart = positions.shape[0]
# pcount is the number of particles per oct.
pcount = np.zeros_like(pdom_ind)
+ oct_left_edges = np.zeros((positions.shape[0], 3), dtype='float64')
+ oct_dds = np.zeros_like(oct_left_edges)
# doff is the offset to a given oct in the sorted particles.
doff = np.zeros_like(pdom_ind) - 1
doff_m = np.zeros((mdom_ind.shape[0], 2), dtype="int64")
@@ -215,7 +219,12 @@
offset = oct.domain_ind - moff_p
pcount[offset] += 1
pdoms[i] = offset # We store the *actual* offset.
- oct = mesh_octree.get(pos)
+ oct = mesh_octree.get(pos, &oinfo)
+ # store oct positions and dds to avoid searching for neighbors
+ # in octs that we know are too far away
+ for j in range(3):
+ oct_left_edges[i, j] = oinfo.left_edge[j]
+ oct_dds[i, j] = oinfo.dds[j]
offset = oct.domain_ind - moff_m
# Now we have oct assignments. Let's sort them.
# Note that what we will be providing to our processing functions will
@@ -255,10 +264,11 @@
visited[oct.domain_ind - moff_m] = 1
if offset < 0: continue
nproc += 1
- self.neighbor_process(dims, moi.left_edge, moi.dds,
- cart_pos, field_pointers, doffs, &nind,
- pinds, pcounts, offset, index_field_pointers,
- particle_octree, domain_id, &nsize)
+ self.neighbor_process(
+ dims, moi.left_edge, moi.dds, cart_pos, field_pointers, doffs,
+ &nind, pinds, pcounts, offset, index_field_pointers,
+ particle_octree, domain_id, &nsize, &oct_left_edges[0, 0],
+ &oct_dds[0, 0])
#print "VISITED", visited.sum(), visited.size,
#print 100.0*float(visited.sum())/visited.size
if nind != NULL:
@@ -524,7 +534,8 @@
np.int64_t *pinds,
np.float64_t *ppos,
np.float64_t cpos[3],
- np.float64_t *dds,
+ np.float64_t *oct_left_edges,
+ np.float64_t *oct_dds,
):
# We are now given the number of neighbors, the indices into the
# domains for them, and the number of particles for each.
@@ -532,19 +543,20 @@
cdef np.int64_t offset, pn, pc
cdef np.float64_t pos[3], lpos[3], rpos[3], r2_trunc, lr2, rr2
self.neighbor_reset()
- # terminate early if left and right edge r2 for this oct
- # is bigger than final r2
- if dds != NULL and self.curn == self.maxn:
- for j in range(3):
- lpos[j] = cpos[j] - dds[j]
- rpos[j] = cpos[j] + dds[j]
- r2_trunc = self.neighbors[self.curn - 1].r2
- lr2 = r2dist(ppos, lpos, self.DW, self.periodicity, r2_trunc)
- rr2 = r2dist(ppos, rpos, self.DW, self.periodicity, r2_trunc)
- if lr2 == -1 and rr2 == -1:
- return
for ni in range(nneighbors):
if nind[ni] == -1: continue
+ # terminate early if left and right edge r2 for this oct
+ # is bigger than final r2
+ if oct_left_edges != NULL and self.curn == self.maxn:
+ for j in range(3):
+ lpos[j] = oct_left_edges[3*nind[ni] + j]
+ rpos[j] = lpos[j] + oct_dds[3*nind[ni] + j]
+ r2_trunc = self.neighbors[self.curn - 1].r2
+ lr2 = r2dist(ppos, lpos, self.DW, self.periodicity, r2_trunc)
+ rr2 = r2dist(ppos, rpos, self.DW, self.periodicity, r2_trunc)
+ if lr2 == -1 and rr2 == -1:
+ continue
+
offset = doffs[nind[ni]]
pc = pcounts[nind[ni]]
for i in range(pc):
@@ -561,7 +573,8 @@
np.int64_t offset,
np.float64_t **index_fields,
OctreeContainer octree, np.int64_t domain_id,
- int *nsize):
+ int *nsize, np.float64_t *oct_left_edges,
+ np.float64_t *oct_dds):
# 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.
@@ -579,7 +592,7 @@
nneighbors = self.neighbor_search(opos, octree,
nind, nsize, nneighbors, domain_id, &oct, 1)
self.neighbor_find(nneighbors, nind[0], doffs, pcounts,
- pinds, ppos, opos, dds)
+ pinds, ppos, opos, oct_left_edges, oct_dds)
# Now we have all our neighbors in our neighbor list.
if self.curn <-1*self.maxn:
ntot = nntot = 0
@@ -617,7 +630,7 @@
nneighbors = self.neighbor_search(opos, octree,
nind, nsize, nneighbors, domain_id, &oct, 1)
self.neighbor_find(nneighbors, nind[0], doffs, pcounts, pinds, ppos,
- opos, NULL)
+ opos, NULL, NULL)
self.process(offset, i, j, k, dim, opos, fields, index_fields)
cdef class VolumeWeightedSmooth(ParticleSmoothOperation):
https://bitbucket.org/yt_analysis/yt/commits/910e220f8d81/
Changeset: 910e220f8d81
Branch: yt
User: ngoldbaum
Date: 2015-09-04 16:23:12+00:00
Summary: Check all 8 corners
Affected #: 2 files
diff -r 2aec0c5340ddca8513a345313b252bad1388d96a -r 910e220f8d81e534b8dbcbdd8e0a45136e243540 yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -32,6 +32,20 @@
np.int64_t pn # Particle number
np.float64_t r2 # radius**2
+cdef np.int64_t corner_permutations[8][3]
+
+corner_permutations = \
+ [[0, 0, 0],
+ [0, 0, 1],
+ [0, 1, 0],
+ [0, 1, 1],
+ [0, 1, 1],
+ [1, 0, 0],
+ [1, 0, 1],
+ [1, 1, 0],
+ [1, 1, 1]]
+
+
cdef class ParticleSmoothOperation:
# We assume each will allocate and define their own temporary storage
cdef kernel_func sph_kernel
diff -r 2aec0c5340ddca8513a345313b252bad1388d96a -r 910e220f8d81e534b8dbcbdd8e0a45136e243540 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -540,21 +540,25 @@
# We are now given the number of neighbors, the indices into the
# domains for them, and the number of particles for each.
cdef int ni, i, j
- cdef np.int64_t offset, pn, pc
- cdef np.float64_t pos[3], lpos[3], rpos[3], r2_trunc, lr2, rr2
+ cdef np.int64_t offset, pn, pc, bad_corner[8]
+ cdef np.float64_t pos[3], corner_pos[3], r2_trunc, lr2, rr2
self.neighbor_reset()
for ni in range(nneighbors):
if nind[ni] == -1: continue
# terminate early if left and right edge r2 for this oct
# is bigger than final r2
if oct_left_edges != NULL and self.curn == self.maxn:
- for j in range(3):
- lpos[j] = oct_left_edges[3*nind[ni] + j]
- rpos[j] = lpos[j] + oct_dds[3*nind[ni] + j]
r2_trunc = self.neighbors[self.curn - 1].r2
- lr2 = r2dist(ppos, lpos, self.DW, self.periodicity, r2_trunc)
- rr2 = r2dist(ppos, rpos, self.DW, self.periodicity, r2_trunc)
- if lr2 == -1 and rr2 == -1:
+ for j in range(8):
+ for k in range(3):
+ corner_pos[k] = oct_left_edges[3*nind[ni] + k]
+ corner_pos[k] += \
+ corner_permutations[j][k] * oct_dds[3*nind[ni] + k]
+ r2 = r2dist(ppos, corner_pos, self.DW, self.periodicity,
+ r2_trunc)
+ if r2 != -1:
+ break
+ if r2 == -1:
continue
offset = doffs[nind[ni]]
https://bitbucket.org/yt_analysis/yt/commits/fc58b430583b/
Changeset: fc58b430583b
Branch: yt
User: MatthewTurk
Date: 2015-09-08 15:10:45+00:00
Summary: Changing early-termination to use all eight corners
Affected #: 1 file
diff -r 910e220f8d81e534b8dbcbdd8e0a45136e243540 -r fc58b430583bd04b4e91243b5afb02efb4894306 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -539,7 +539,7 @@
):
# We are now given the number of neighbors, the indices into the
# domains for them, and the number of particles for each.
- cdef int ni, i, j
+ cdef int ni, i, j, out_of_range
cdef np.int64_t offset, pn, pc, bad_corner[8]
cdef np.float64_t pos[3], corner_pos[3], r2_trunc, lr2, rr2
self.neighbor_reset()
@@ -549,6 +549,7 @@
# is bigger than final r2
if oct_left_edges != NULL and self.curn == self.maxn:
r2_trunc = self.neighbors[self.curn - 1].r2
+ out_of_range = 0
for j in range(8):
for k in range(3):
corner_pos[k] = oct_left_edges[3*nind[ni] + k]
@@ -557,10 +558,9 @@
r2 = r2dist(ppos, corner_pos, self.DW, self.periodicity,
r2_trunc)
if r2 != -1:
- break
- if r2 == -1:
- continue
-
+ out_of_range += 1
+ if out_of_range == 8:
+ continue
offset = doffs[nind[ni]]
pc = pcounts[nind[ni]]
for i in range(pc):
https://bitbucket.org/yt_analysis/yt/commits/e1ef1cb64209/
Changeset: e1ef1cb64209
Branch: yt
User: ngoldbaum
Date: 2015-09-08 16:04:22+00:00
Summary: Remove profiling header
Affected #: 1 file
diff -r fc58b430583bd04b4e91243b5afb02efb4894306 -r e1ef1cb64209d56db7d643abdebf95706addeb07 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -1,4 +1,3 @@
-# cython: profile=True
"""
Particle smoothing in cells
https://bitbucket.org/yt_analysis/yt/commits/7784b13b739c/
Changeset: 7784b13b739c
Branch: yt
User: ngoldbaum
Date: 2015-09-08 16:05:04+00:00
Summary: Ensure all variables have types to avoid python overhead
Affected #: 1 file
diff -r e1ef1cb64209d56db7d643abdebf95706addeb07 -r 7784b13b739cd819a05eb2e899f8005902394aab yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -538,9 +538,9 @@
):
# We are now given the number of neighbors, the indices into the
# domains for them, and the number of particles for each.
- cdef int ni, i, j, out_of_range
+ cdef int ni, i, j, k, out_of_range
cdef np.int64_t offset, pn, pc, bad_corner[8]
- cdef np.float64_t pos[3], corner_pos[3], r2_trunc, lr2, rr2
+ cdef np.float64_t pos[3], corner_pos[3], r2_trunc, r2
self.neighbor_reset()
for ni in range(nneighbors):
if nind[ni] == -1: continue
https://bitbucket.org/yt_analysis/yt/commits/1bde4640c1e6/
Changeset: 1bde4640c1e6
Branch: yt
User: ngoldbaum
Date: 2015-09-08 16:05:34+00:00
Summary: Add early termination to avoid unnnecessarily checking all 8 corners
Affected #: 1 file
diff -r 7784b13b739cd819a05eb2e899f8005902394aab -r 1bde4640c1e60ace26dd3e38a3db0a8ceec0d9d0 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -558,6 +558,8 @@
r2_trunc)
if r2 != -1:
out_of_range += 1
+ else:
+ break
if out_of_range == 8:
continue
offset = doffs[nind[ni]]
https://bitbucket.org/yt_analysis/yt/commits/0e2417cc7db4/
Changeset: 0e2417cc7db4
Branch: yt
User: ngoldbaum
Date: 2015-09-08 16:10:06+00:00
Summary: Add better comment
Affected #: 1 file
diff -r 1bde4640c1e60ace26dd3e38a3db0a8ceec0d9d0 -r 0e2417cc7db40108aa8a26393016574e43d32351 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -544,8 +544,8 @@
self.neighbor_reset()
for ni in range(nneighbors):
if nind[ni] == -1: continue
- # terminate early if left and right edge r2 for this oct
- # is bigger than final r2
+ # terminate early if all 8 corners of oct are farther away than
+ # most distant currently known neighbor
if oct_left_edges != NULL and self.curn == self.maxn:
r2_trunc = self.neighbors[self.curn - 1].r2
out_of_range = 0
https://bitbucket.org/yt_analysis/yt/commits/465514c5db2b/
Changeset: 465514c5db2b
Branch: yt
User: MatthewTurk
Date: 2015-09-08 18:42:03+00:00
Summary: Remove corner checks, add quick dist check, fix mesh/particle disparity.
Affected #: 2 files
diff -r 0e2417cc7db40108aa8a26393016574e43d32351 -r 465514c5db2bb36d2dcd99a8734cc59b96abccf5 yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -32,20 +32,6 @@
np.int64_t pn # Particle number
np.float64_t r2 # radius**2
-cdef np.int64_t corner_permutations[8][3]
-
-corner_permutations = \
- [[0, 0, 0],
- [0, 0, 1],
- [0, 1, 0],
- [0, 1, 1],
- [0, 1, 1],
- [1, 0, 0],
- [1, 0, 1],
- [1, 1, 0],
- [1, 1, 1]]
-
-
cdef class ParticleSmoothOperation:
# We assume each will allocate and define their own temporary storage
cdef kernel_func sph_kernel
diff -r 0e2417cc7db40108aa8a26393016574e43d32351 -r 465514c5db2bb36d2dcd99a8734cc59b96abccf5 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -179,7 +179,7 @@
numpart = positions.shape[0]
# pcount is the number of particles per oct.
pcount = np.zeros_like(pdom_ind)
- oct_left_edges = np.zeros((positions.shape[0], 3), dtype='float64')
+ oct_left_edges = np.zeros((pdom_ind.shape[0], 3), dtype='float64')
oct_dds = np.zeros_like(oct_left_edges)
# doff is the offset to a given oct in the sorted particles.
doff = np.zeros_like(pdom_ind) - 1
@@ -206,10 +206,11 @@
for i in range(3):
self.DW[i] = (mesh_octree.DRE[i] - mesh_octree.DLE[i])
self.periodicity[i] = periodicity[i]
+ cdef np.float64_t factor = (1 << (particle_octree.oref))
for i in range(positions.shape[0]):
for j in range(3):
pos[j] = positions[i, j]
- oct = particle_octree.get(pos)
+ oct = particle_octree.get(pos, &oinfo)
if oct == NULL or (domain_id > 0 and oct.domain != domain_id):
continue
# Note that this has to be our local index, not our in-file index.
@@ -218,13 +219,11 @@
offset = oct.domain_ind - moff_p
pcount[offset] += 1
pdoms[i] = offset # We store the *actual* offset.
- oct = mesh_octree.get(pos, &oinfo)
# store oct positions and dds to avoid searching for neighbors
# in octs that we know are too far away
for j in range(3):
- oct_left_edges[i, j] = oinfo.left_edge[j]
- oct_dds[i, j] = oinfo.dds[j]
- offset = oct.domain_ind - moff_m
+ oct_left_edges[offset, j] = oinfo.left_edge[j]
+ oct_dds[offset, j] = oinfo.dds[j] * factor
# Now we have oct assignments. Let's sort them.
# Note that what we will be providing to our processing functions will
# actually be indirectly-sorted fields. This preserves memory at the
@@ -538,9 +537,9 @@
):
# We are now given the number of neighbors, the indices into the
# domains for them, and the number of particles for each.
- cdef int ni, i, j, k, out_of_range
- cdef np.int64_t offset, pn, pc, bad_corner[8]
- cdef np.float64_t pos[3], corner_pos[3], r2_trunc, r2
+ cdef int ni, i, j, k
+ cdef np.int64_t offset, pn, pc
+ cdef np.float64_t pos[3], cp, r2_trunc, r2, ex[2], DR[2], dist
self.neighbor_reset()
for ni in range(nneighbors):
if nind[ni] == -1: continue
@@ -548,20 +547,35 @@
# most distant currently known neighbor
if oct_left_edges != NULL and self.curn == self.maxn:
r2_trunc = self.neighbors[self.curn - 1].r2
- out_of_range = 0
- for j in range(8):
- for k in range(3):
- corner_pos[k] = oct_left_edges[3*nind[ni] + k]
- corner_pos[k] += \
- corner_permutations[j][k] * oct_dds[3*nind[ni] + k]
- r2 = r2dist(ppos, corner_pos, self.DW, self.periodicity,
- r2_trunc)
- if r2 != -1:
- out_of_range += 1
- else:
- break
- if out_of_range == 8:
- continue
+ # iterate over each dimension in the outer loop so we can
+ # consolidate temporary storage
+ # What this next bit does is figure out which component is the
+ # closest, of each possible permutation.
+ # k here is the dimension
+ r2 = 0.0
+ for k in range(3):
+ # We start at left edge, then do halfway, then right edge.
+ ex[0] = oct_left_edges[3*nind[ni] + k]
+ ex[1] = ex[0] + oct_dds[3*nind[ni] + k]
+ # There are three possibilities; we are between, left-of,
+ # or right-of the extrema. Thanks to
+ # http://stackoverflow.com/questions/5254838/calculating-distance-between-a-point-and-a-rectangular-box-nearest-point
+ # for some help. This has been modified to account for
+ # periodicity.
+ dist = 0.0
+ DR[0] = (ex[0] - cpos[k])
+ DR[1] = (cpos[k] - ex[1])
+ for j in range(2):
+ if not self.periodicity[k]:
+ pass
+ elif (DR[j] > self.DW[k]/2.0):
+ DR[j] -= self.DW[k]
+ elif (DR[j] < -self.DW[k]/2.0):
+ DR[j] += self.DW[k]
+ dist = fmax(dist, DR[j])
+ r2 += dist*dist
+ if r2 > r2_trunc:
+ continue
offset = doffs[nind[ni]]
pc = pcounts[nind[ni]]
for i in range(pc):
https://bitbucket.org/yt_analysis/yt/commits/a653b6ad133a/
Changeset: a653b6ad133a
Branch: yt
User: MatthewTurk
Date: 2015-09-08 22:27:47+00:00
Summary: We do not need an extra layer in almost all cases.
Affected #: 1 file
diff -r 465514c5db2bb36d2dcd99a8734cc59b96abccf5 -r a653b6ad133ade1b1db331223a42858d5186e965 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -609,7 +609,7 @@
for k in range(dim[2]):
self.pos_setup(cpos, opos)
nneighbors = self.neighbor_search(opos, octree,
- nind, nsize, nneighbors, domain_id, &oct, 1)
+ nind, nsize, nneighbors, domain_id, &oct, 0)
self.neighbor_find(nneighbors, nind[0], doffs, pcounts,
pinds, ppos, opos, oct_left_edges, oct_dds)
# Now we have all our neighbors in our neighbor list.
@@ -647,7 +647,7 @@
cdef np.float64_t opos[3]
self.pos_setup(cpos, opos)
nneighbors = self.neighbor_search(opos, octree,
- nind, nsize, nneighbors, domain_id, &oct, 1)
+ nind, nsize, nneighbors, domain_id, &oct, 0)
self.neighbor_find(nneighbors, nind[0], doffs, pcounts, pinds, ppos,
opos, NULL, NULL)
self.process(offset, i, j, k, dim, opos, fields, index_fields)
https://bitbucket.org/yt_analysis/yt/commits/07bdfafdedda/
Changeset: 07bdfafdedda
Branch: yt
User: MatthewTurk
Date: 2015-09-09 14:23:57+00:00
Summary: Add comments, remove weight field.
Affected #: 1 file
diff -r a653b6ad133ade1b1db331223a42858d5186e965 -r 07bdfafdedda39f8324344b0df83105645854e33 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -653,6 +653,12 @@
self.process(offset, i, j, k, dim, opos, fields, index_fields)
cdef class VolumeWeightedSmooth(ParticleSmoothOperation):
+ # This smoothing function evaluates the field, *without* normalization, at
+ # every point in the *mesh*. Applying a normalization results in
+ # non-conservation of mass when smoothing density; to avoid this, we do not
+ # apply this normalization factor. The SPLASH paper
+ # (http://arxiv.org/abs/0709.0832v1) discusses this in detail; what we are
+ # applying here is equation 6, with variable smoothing lengths (eq 13).
cdef np.float64_t **fp
cdef public object vals
def initialize(self):
@@ -664,18 +670,17 @@
raise RuntimeError
cdef np.ndarray tarr
self.fp = <np.float64_t **> malloc(
- sizeof(np.float64_t *) * (self.nfields - 2))
+ sizeof(np.float64_t *) * (self.nfields - 3))
self.vals = []
- for i in range(self.nfields - 2):
+ # We usually only allocate one field; if we are doing multiple field,
+ # single-pass smoothing, then we might have more.
+ for i in range(self.nfields - 3):
tarr = np.zeros(self.nvals, dtype="float64", order="F")
self.vals.append(tarr)
self.fp[i] = <np.float64_t *> tarr.data
def finalize(self):
free(self.fp)
- vv = self.vals.pop(-1)
- for v in self.vals:
- v /= vv
return self.vals
@cython.cdivision(True)
@@ -719,9 +724,6 @@
for fi in range(self.nfields - 3):
val = fields[fi + 3][pn]
self.fp[fi][gind(i,j,k,dim) + offset] += val * weight
- self.fp[self.nfields - 3][gind(i,j,k,dim) + offset] += weight
- # When looking for non-normalized values, uncomment:
- self.fp[self.nfields - 3][gind(i,j,k,dim) + offset] = 1.0
return
volume_weighted_smooth = VolumeWeightedSmooth
https://bitbucket.org/yt_analysis/yt/commits/50b4d07251c6/
Changeset: 50b4d07251c6
Branch: yt
User: xarthisius
Date: 2015-09-15 20:21:17+00:00
Summary: Merged in MatthewTurk/yt (pull request #1740)
Remove normalization from smoothing kernel
Affected #: 5 files
diff -r 4979406431c662b038db914287ae227092199b78 -r 50b4d07251c61fabb1a47837041769cdb911bcb7 yt/fields/particle_fields.py
--- a/yt/fields/particle_fields.py
+++ b/yt/fields/particle_fields.py
@@ -788,13 +788,19 @@
kwargs = {}
if nneighbors:
kwargs['nneighbors'] = nneighbors
+ # This is for applying cutoffs, similar to in the SPLASH paper.
+ smooth_cutoff = data["index","cell_volume"]**(1./3)
+ smooth_cutoff.convert_to_units("code_length")
# volume_weighted smooth operations return lists of length 1.
rv = data.smooth(pos, [mass, hsml, dens, quan],
method="volume_weighted",
create_octree=True,
+ index_fields=[smooth_cutoff],
kernel_name=kernel_name)[0]
rv[np.isnan(rv)] = 0.0
# Now some quick unit conversions.
+ # This should be used when seeking a non-normalized value:
+ rv /= hsml.uq**3 / hsml.uq.in_cgs().uq**3
rv = data.apply_units(rv, field_units)
return rv
registry.add_field(field_name, function = _vol_weight,
diff -r 4979406431c662b038db914287ae227092199b78 -r 50b4d07251c61fabb1a47837041769cdb911bcb7 yt/frontends/tipsy/fields.py
--- a/yt/frontends/tipsy/fields.py
+++ b/yt/frontends/tipsy/fields.py
@@ -38,7 +38,8 @@
'FeMassFrac':("FeMassFrac", ("dimensionless", ["Fe_fraction"], None)),
'c':("c", ("code_velocity", [""], None)),
'acc':("acc", ("code_velocity / code_time", [""], None)),
- 'accg':("accg", ("code_velocity / code_time", [""], None))}
+ 'accg':("accg", ("code_velocity / code_time", [""], None)),
+ 'smoothlength':('smoothlength', ("code_length", ["smoothing_length"], None))}
def __init__(self, ds, field_list, slice_info = None):
for field in field_list:
@@ -60,15 +61,19 @@
def setup_gas_particle_fields(self, ptype):
- def _smoothing_length(field, data):
- # For now, we hardcode num_neighbors. We should make this configurable
- # in the future.
- num_neighbors = 64
- fn, = add_nearest_neighbor_field(ptype, "particle_position", self, num_neighbors)
- return data[ptype, 'nearest_neighbor_distance_%d' % num_neighbors]
+ num_neighbors = 65
+ fn, = add_nearest_neighbor_field(ptype, "particle_position", self, num_neighbors)
+ def _func():
+ def _smoothing_length(field, data):
+ # For now, we hardcode num_neighbors. We should make this configurable
+ # in the future.
+ rv = data[ptype, 'nearest_neighbor_distance_%d' % num_neighbors]
+ #np.maximum(rv, 0.5*data[ptype, "Epsilon"], rv)
+ return rv
+ return _smoothing_length
self.add_field(
(ptype, "smoothing_length"),
- function=_smoothing_length,
+ function=_func(),
particle_type=True,
units="code_length")
diff -r 4979406431c662b038db914287ae227092199b78 -r 50b4d07251c61fabb1a47837041769cdb911bcb7 yt/geometry/particle_deposit.pxd
--- a/yt/geometry/particle_deposit.pxd
+++ b/yt/geometry/particle_deposit.pxd
@@ -40,13 +40,14 @@
cdef inline np.float64_t sph_kernel_cubic(np.float64_t x) nogil:
cdef np.float64_t kernel
+ cdef np.float64_t C = 2.5464790894703255
if x <= 0.5:
kernel = 1.-6.*x*x*(1.-x)
elif x>0.5 and x<=1.0:
kernel = 2.*(1.-x)*(1.-x)*(1.-x)
else:
kernel = 0.
- return kernel
+ return kernel * C
########################################################
# Alternative SPH kernels for use with the Grid method #
diff -r 4979406431c662b038db914287ae227092199b78 -r 50b4d07251c61fabb1a47837041769cdb911bcb7 yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -55,10 +55,12 @@
np.int64_t *pinds, np.int64_t *pcounts,
np.int64_t offset, np.float64_t **index_fields,
OctreeContainer octree, np.int64_t domain_id,
- int *nsize)
+ int *nsize, np.float64_t *oct_left_edges,
+ np.float64_t *oct_dds)
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 = ?)
+ np.int64_t nneighbors, np.int64_t domain_id,
+ Oct **oct = ?, int extra_layer = ?)
cdef void neighbor_process_particle(self, np.float64_t cpos[3],
np.float64_t *ppos,
np.float64_t **fields,
@@ -78,7 +80,9 @@
np.int64_t *pcounts,
np.int64_t *pinds,
np.float64_t *ppos,
- np.float64_t cpos[3])
+ np.float64_t cpos[3],
+ np.float64_t* oct_left_edges,
+ np.float64_t* oct_dds)
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)
diff -r 4979406431c662b038db914287ae227092199b78 -r 50b4d07251c61fabb1a47837041769cdb911bcb7 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -154,6 +154,8 @@
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
+ cdef np.ndarray[np.float64_t, ndim=2] oct_left_edges, oct_dds
+ cdef OctInfo oinfo
if geometry == "cartesian":
self.pos_setup = cart_coord_setup
cart_positions = positions
@@ -177,6 +179,8 @@
numpart = positions.shape[0]
# pcount is the number of particles per oct.
pcount = np.zeros_like(pdom_ind)
+ oct_left_edges = np.zeros((pdom_ind.shape[0], 3), dtype='float64')
+ oct_dds = np.zeros_like(oct_left_edges)
# doff is the offset to a given oct in the sorted particles.
doff = np.zeros_like(pdom_ind) - 1
doff_m = np.zeros((mdom_ind.shape[0], 2), dtype="int64")
@@ -202,10 +206,11 @@
for i in range(3):
self.DW[i] = (mesh_octree.DRE[i] - mesh_octree.DLE[i])
self.periodicity[i] = periodicity[i]
+ cdef np.float64_t factor = (1 << (particle_octree.oref))
for i in range(positions.shape[0]):
for j in range(3):
pos[j] = positions[i, j]
- oct = particle_octree.get(pos)
+ oct = particle_octree.get(pos, &oinfo)
if oct == NULL or (domain_id > 0 and oct.domain != domain_id):
continue
# Note that this has to be our local index, not our in-file index.
@@ -214,8 +219,11 @@
offset = oct.domain_ind - moff_p
pcount[offset] += 1
pdoms[i] = offset # We store the *actual* offset.
- oct = mesh_octree.get(pos)
- offset = oct.domain_ind - moff_m
+ # store oct positions and dds to avoid searching for neighbors
+ # in octs that we know are too far away
+ for j in range(3):
+ oct_left_edges[offset, j] = oinfo.left_edge[j]
+ oct_dds[offset, j] = oinfo.dds[j] * factor
# Now we have oct assignments. Let's sort them.
# Note that what we will be providing to our processing functions will
# actually be indirectly-sorted fields. This preserves memory at the
@@ -254,10 +262,11 @@
visited[oct.domain_ind - moff_m] = 1
if offset < 0: continue
nproc += 1
- self.neighbor_process(dims, moi.left_edge, moi.dds,
- cart_pos, field_pointers, doffs, &nind,
- pinds, pcounts, offset, index_field_pointers,
- particle_octree, domain_id, &nsize)
+ self.neighbor_process(
+ dims, moi.left_edge, moi.dds, cart_pos, field_pointers, doffs,
+ &nind, pinds, pcounts, offset, index_field_pointers,
+ particle_octree, domain_id, &nsize, &oct_left_edges[0, 0],
+ &oct_dds[0, 0])
#print "VISITED", visited.sum(), visited.size,
#print 100.0*float(visited.sum())/visited.size
if nind != NULL:
@@ -394,11 +403,12 @@
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):
+ Oct **oct = NULL, int extra_layer = 0):
cdef OctInfo oi
cdef Oct *ooct
- cdef Oct **neighbors
- cdef int j
+ cdef Oct **neighbors, **first_layer
+ cdef int j, total_neighbors = 0, initial_layer = 0
+ cdef int layer_ind = 0
cdef np.int64_t moff = octree.get_domain_offset(domain_id)
ooct = octree.get(pos, &oi)
if oct != NULL and ooct == oct[0]:
@@ -407,24 +417,53 @@
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, self.periodicity)
- # 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
+ # This is our "seed" set of neighbors. If we are asked to, we will
+ # create a master list of neighbors that is much bigger and includes
+ # everything.
+ layer_ind = 0
+ first_layer = NULL
+ while 1:
+ neighbors = octree.neighbors(&oi, &nneighbors, ooct, self.periodicity)
+ # Now we have all our neighbors. And, we should be set for what
+ # else we need to do.
+ if total_neighbors + nneighbors > nsize[0]:
+ nind[0] = <np.int64_t *> realloc(
+ nind[0], sizeof(np.int64_t)*(nneighbors + total_neighbors))
+ nsize[0] = nneighbors + total_neighbors
+ for j in range(nneighbors):
+ # Particle octree neighbor indices
+ nind[0][j + total_neighbors] = neighbors[j].domain_ind - moff
+ total_neighbors += nneighbors
+ if extra_layer == 0:
+ # Not adding on any additional layers here.
+ free(neighbors)
+ neighbors = NULL
+ break
+ if initial_layer == 0:
+ initial_layer = nneighbors
+ first_layer = neighbors
+ else:
+ # Allocated internally; we free this in the loops if we aren't
+ # tracking it
+ free(neighbors)
+ neighbors = NULL
+ ooct = first_layer[layer_ind]
+ layer_ind += 1
+ if layer_ind == initial_layer:
+ neighbors
+ break
+
- for j in range(nneighbors):
+ for j in range(total_neighbors):
# Particle octree neighbor indices
- nind[0][j] = neighbors[j].domain_ind - moff
+ if nind[0][j] == -1: continue
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
+ if first_layer != NULL:
+ free(first_layer)
+ return total_neighbors
@cython.cdivision(True)
@cython.boundscheck(False)
@@ -492,16 +531,51 @@
np.int64_t *pcounts,
np.int64_t *pinds,
np.float64_t *ppos,
- np.float64_t cpos[3]
+ np.float64_t cpos[3],
+ np.float64_t *oct_left_edges,
+ np.float64_t *oct_dds,
):
# We are now given the number of neighbors, the indices into the
# domains for them, and the number of particles for each.
- cdef int ni, i, j
+ cdef int ni, i, j, k
cdef np.int64_t offset, pn, pc
- cdef np.float64_t pos[3]
+ cdef np.float64_t pos[3], cp, r2_trunc, r2, ex[2], DR[2], dist
self.neighbor_reset()
for ni in range(nneighbors):
if nind[ni] == -1: continue
+ # terminate early if all 8 corners of oct are farther away than
+ # most distant currently known neighbor
+ if oct_left_edges != NULL and self.curn == self.maxn:
+ r2_trunc = self.neighbors[self.curn - 1].r2
+ # iterate over each dimension in the outer loop so we can
+ # consolidate temporary storage
+ # What this next bit does is figure out which component is the
+ # closest, of each possible permutation.
+ # k here is the dimension
+ r2 = 0.0
+ for k in range(3):
+ # We start at left edge, then do halfway, then right edge.
+ ex[0] = oct_left_edges[3*nind[ni] + k]
+ ex[1] = ex[0] + oct_dds[3*nind[ni] + k]
+ # There are three possibilities; we are between, left-of,
+ # or right-of the extrema. Thanks to
+ # http://stackoverflow.com/questions/5254838/calculating-distance-between-a-point-and-a-rectangular-box-nearest-point
+ # for some help. This has been modified to account for
+ # periodicity.
+ dist = 0.0
+ DR[0] = (ex[0] - cpos[k])
+ DR[1] = (cpos[k] - ex[1])
+ for j in range(2):
+ if not self.periodicity[k]:
+ pass
+ elif (DR[j] > self.DW[k]/2.0):
+ DR[j] -= self.DW[k]
+ elif (DR[j] < -self.DW[k]/2.0):
+ DR[j] += self.DW[k]
+ dist = fmax(dist, DR[j])
+ r2 += dist*dist
+ if r2 > r2_trunc:
+ continue
offset = doffs[nind[ni]]
pc = pcounts[nind[ni]]
for i in range(pc):
@@ -518,7 +592,8 @@
np.int64_t offset,
np.float64_t **index_fields,
OctreeContainer octree, np.int64_t domain_id,
- int *nsize):
+ int *nsize, np.float64_t *oct_left_edges,
+ np.float64_t *oct_dds):
# 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.
@@ -534,9 +609,9 @@
for k in range(dim[2]):
self.pos_setup(cpos, opos)
nneighbors = self.neighbor_search(opos, octree,
- nind, nsize, nneighbors, domain_id, &oct)
+ nind, nsize, nneighbors, domain_id, &oct, 0)
self.neighbor_find(nneighbors, nind[0], doffs, pcounts,
- pinds, ppos, opos)
+ pinds, ppos, opos, oct_left_edges, oct_dds)
# Now we have all our neighbors in our neighbor list.
if self.curn <-1*self.maxn:
ntot = nntot = 0
@@ -572,11 +647,18 @@
cdef np.float64_t opos[3]
self.pos_setup(cpos, opos)
nneighbors = self.neighbor_search(opos, octree,
- nind, nsize, nneighbors, domain_id, &oct)
- self.neighbor_find(nneighbors, nind[0], doffs, pcounts, pinds, ppos, opos)
+ nind, nsize, nneighbors, domain_id, &oct, 0)
+ self.neighbor_find(nneighbors, nind[0], doffs, pcounts, pinds, ppos,
+ opos, NULL, NULL)
self.process(offset, i, j, k, dim, opos, fields, index_fields)
cdef class VolumeWeightedSmooth(ParticleSmoothOperation):
+ # This smoothing function evaluates the field, *without* normalization, at
+ # every point in the *mesh*. Applying a normalization results in
+ # non-conservation of mass when smoothing density; to avoid this, we do not
+ # apply this normalization factor. The SPLASH paper
+ # (http://arxiv.org/abs/0709.0832v1) discusses this in detail; what we are
+ # applying here is equation 6, with variable smoothing lengths (eq 13).
cdef np.float64_t **fp
cdef public object vals
def initialize(self):
@@ -588,18 +670,17 @@
raise RuntimeError
cdef np.ndarray tarr
self.fp = <np.float64_t **> malloc(
- sizeof(np.float64_t *) * (self.nfields - 2))
+ sizeof(np.float64_t *) * (self.nfields - 3))
self.vals = []
- for i in range(self.nfields - 2):
+ # We usually only allocate one field; if we are doing multiple field,
+ # single-pass smoothing, then we might have more.
+ for i in range(self.nfields - 3):
tarr = np.zeros(self.nvals, dtype="float64", order="F")
self.vals.append(tarr)
self.fp[i] = <np.float64_t *> tarr.data
def finalize(self):
free(self.fp)
- vv = self.vals.pop(-1)
- for v in self.vals:
- v /= vv
return self.vals
@cython.cdivision(True)
@@ -612,11 +693,13 @@
# We also have a list of neighboring particles with particle numbers.
cdef int n, fi
cdef np.float64_t weight, r2, val, hsml, dens, mass, coeff, max_r
+ cdef np.float64_t max_hsml, ihsml, ihsml3, kern
coeff = 0.0
cdef np.int64_t pn
# We get back our mass
# rho_i = sum(j = 1 .. n) m_j * W_ij
max_r = sqrt(self.neighbors[self.curn-1].r2)
+ max_hsml = index_fields[0][gind(i,j,k,dim) + offset]
for n in range(self.curn):
# No normalization for the moment.
# fields[0] is the smoothing length.
@@ -625,18 +708,22 @@
# Smoothing kernel weight function
mass = fields[0][pn]
hsml = fields[1][pn]
+ dens = fields[2][pn]
if hsml < 0:
hsml = max_r
if hsml == 0: continue
+ ihsml = 1.0/hsml
+ hsml = fmax(max_hsml/2.0, hsml)
+ ihsml3 = ihsml*ihsml*ihsml
# Usually this density has been computed
- dens = fields[2][pn]
if dens == 0.0: continue
- weight = mass * self.sph_kernel(sqrt(r2) / hsml) / dens
+ weight = (mass / dens) * ihsml3
+ kern = self.sph_kernel(sqrt(r2) * ihsml)
+ weight *= kern
# Mass of the particle times the value
for fi in range(self.nfields - 3):
val = fields[fi + 3][pn]
self.fp[fi][gind(i,j,k,dim) + offset] += val * weight
- self.fp[self.nfields - 3][gind(i,j,k,dim) + offset] += weight
return
volume_weighted_smooth = VolumeWeightedSmooth
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