<div dir="ltr">Hey Ben,<div><br></div><div>What happened here? Did you mean to merge this PR? I think it was still marked as work in progress, e.g. [WIP], and Billi had more commits to push to it.</div><div><br></div><div>-Nathan</div><div><br><div class="gmail_quote">---------- Forwarded message ----------<br>From: <b class="gmail_sendername"></b> <span dir="ltr"><<a href="mailto:commits-noreply@bitbucket.org">commits-noreply@bitbucket.org</a>></span><br>Date: Wed, Aug 19, 2015 at 12:09 PM<br>Subject: [yt-svn] commit/yt: bwkeller: Merged in qobilidop/yt (pull request #1670)<br>To: <a href="mailto:yt-svn@lists.spacepope.org">yt-svn@lists.spacepope.org</a><br><br><br>1 new commit in yt:<br>
<br>
<a href="https://bitbucket.org/yt_analysis/yt/commits/cfc98fddd9d7/" rel="noreferrer" target="_blank">https://bitbucket.org/yt_analysis/yt/commits/cfc98fddd9d7/</a><br>
Changeset: cfc98fddd9d7<br>
Branch: yt<br>
User: bwkeller<br>
Date: 2015-08-19 17:09:19+00:00<br>
Summary: Merged in qobilidop/yt (pull request #1670)<br>
<br>
Alternative Smoothing Kernels<br>
Affected #: 11 files<br>
<br>
diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/data_objects/construction_data_containers.py<br>
--- a/yt/data_objects/construction_data_containers.py<br>
+++ b/yt/data_objects/construction_data_containers.py<br>
@@ -682,11 +682,13 @@<br>
def RightEdge(self):<br>
return self.right_edge<br>
<br>
- def deposit(self, positions, fields = None, method = None):<br>
+ def deposit(self, positions, fields = None, method = None,<br>
+ kernel_name = 'cubic'):<br>
cls = getattr(particle_deposit, "deposit_%s" % method, None)<br>
if cls is None:<br>
raise YTParticleDepositionNotImplemented(method)<br>
- op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs<br>
+ # We allocate number of zones, not number of octs<br>
+ op = cls(self.ActiveDimensions.prod(), kernel_name)<br>
op.initialize()<br>
op.process_grid(self, positions, fields)<br>
vals = op.finalize()<br>
@@ -1787,5 +1789,3 @@<br>
else:<br>
mylog.error("Problem uploading.")<br>
return upload_id<br>
-<br>
-<br>
<br>
diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/data_objects/grid_patch.py<br>
--- a/yt/data_objects/grid_patch.py<br>
+++ b/yt/data_objects/grid_patch.py<br>
@@ -330,12 +330,14 @@<br>
def particle_operation(self, *args, **kwargs):<br>
raise NotImplementedError<br>
<br>
- def deposit(self, positions, fields = None, method = None):<br>
+ def deposit(self, positions, fields = None, method = None,<br>
+ kernel_name = 'cubic'):<br>
# Here we perform our particle deposition.<br>
cls = getattr(particle_deposit, "deposit_%s" % method, None)<br>
if cls is None:<br>
raise YTParticleDepositionNotImplemented(method)<br>
- op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs<br>
+ # We allocate number of zones, not number of octs<br>
+ op = cls(self.ActiveDimensions.prod(), kernel_name)<br>
op.initialize()<br>
op.process_grid(self, positions, fields)<br>
vals = op.finalize()<br>
<br>
diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/data_objects/octree_subset.py<br>
--- a/yt/data_objects/octree_subset.py<br>
+++ b/yt/data_objects/octree_subset.py<br>
@@ -141,7 +141,8 @@<br>
self._domain_ind = di<br>
return self._domain_ind<br>
<br>
- def deposit(self, positions, fields = None, method = None):<br>
+ def deposit(self, positions, fields = None, method = None,<br>
+ kernel_name='cubic'):<br>
r"""Operate on the mesh, in a particle-against-mesh fashion, with<br>
exclusively local input.<br>
<br>
@@ -176,7 +177,8 @@<br>
raise YTParticleDepositionNotImplemented(method)<br>
nz = <a href="http://self.nz" rel="noreferrer" target="_blank">self.nz</a><br>
nvals = (nz, nz, nz, (self.domain_ind >= 0).sum())<br>
- op = cls(nvals) # We allocate number of zones, not number of octs<br>
+ # We allocate number of zones, not number of octs<br>
+ op = cls(nvals, kernel_name)<br>
op.initialize()<br>
mylog.debug("Depositing %s (%s^3) particles into %s Octs",<br>
positions.shape[0], positions.shape[0]**0.3333333, nvals[-1])<br>
@@ -192,7 +194,8 @@<br>
return np.asfortranarray(vals)<br>
<br>
def smooth(self, positions, fields = None, index_fields = None,<br>
- method = None, create_octree = False, nneighbors = 64):<br>
+ method = None, create_octree = False, nneighbors = 64,<br>
+ kernel_name = 'cubic'):<br>
r"""Operate on the mesh, in a particle-against-mesh fashion, with<br>
non-local input.<br>
<br>
@@ -258,7 +261,7 @@<br>
nz = <a href="http://self.nz" rel="noreferrer" target="_blank">self.nz</a><br>
mdom_ind = self.domain_ind<br>
nvals = (nz, nz, nz, (mdom_ind >= 0).sum())<br>
- op = cls(nvals, len(fields), nneighbors)<br>
+ op = cls(nvals, len(fields), nneighbors, kernel_name)<br>
op.initialize()<br>
mylog.debug("Smoothing %s particles into %s Octs",<br>
positions.shape[0], nvals[-1])<br>
@@ -280,7 +283,7 @@<br>
return vals<br>
<br>
def particle_operation(self, positions, fields = None,<br>
- method = None, nneighbors = 64):<br>
+ method = None, nneighbors = 64, kernel_name = 'cubic'):<br>
r"""Operate on particles, in a particle-against-particle fashion.<br>
<br>
This uses the octree indexing system to call a "smoothing" operation<br>
@@ -335,7 +338,7 @@<br>
nz = <a href="http://self.nz" rel="noreferrer" target="_blank">self.nz</a><br>
mdom_ind = self.domain_ind<br>
nvals = (nz, nz, nz, (mdom_ind >= 0).sum())<br>
- op = cls(nvals, len(fields), nneighbors)<br>
+ op = cls(nvals, len(fields), nneighbors, kernel_name)<br>
op.initialize()<br>
mylog.debug("Smoothing %s particles into %s Octs",<br>
positions.shape[0], nvals[-1])<br>
<br>
diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/data_objects/unstructured_mesh.py<br>
--- a/yt/data_objects/unstructured_mesh.py<br>
+++ b/yt/data_objects/unstructured_mesh.py<br>
@@ -105,13 +105,15 @@<br>
def select_tcoords(self, dobj):<br>
raise NotImplementedError<br>
<br>
- def deposit(self, positions, fields = None, method = None):<br>
+ def deposit(self, positions, fields = None, method = None,<br>
+ kernel_name = 'cubic'):<br>
raise NotImplementedError<br>
# Here we perform our particle deposition.<br>
cls = getattr(particle_deposit, "deposit_%s" % method, None)<br>
if cls is None:<br>
raise YTParticleDepositionNotImplemented(method)<br>
- op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs<br>
+ # We allocate number of zones, not number of octs<br>
+ op = cls(self.ActiveDimensions.prod(), kernel_name)<br>
op.initialize()<br>
op.process_grid(self, positions, fields)<br>
vals = op.finalize()<br>
@@ -200,4 +202,3 @@<br>
else:<br>
self._last_count = mask.sum()<br>
return mask<br>
-<br>
<br>
diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/fields/particle_fields.py<br>
--- a/yt/fields/particle_fields.py<br>
+++ b/yt/fields/particle_fields.py<br>
@@ -766,8 +766,12 @@<br>
<br>
def add_volume_weighted_smoothed_field(ptype, coord_name, mass_name,<br>
smoothing_length_name, density_name, smoothed_field, registry,<br>
- nneighbors = None):<br>
- field_name = ("deposit", "%s_smoothed_%s" % (ptype, smoothed_field))<br>
+ nneighbors = None, kernel_name = 'cubic'):<br>
+ if kernel_name == 'cubic':<br>
+ field_name = ("deposit", "%s_smoothed_%s" % (ptype, smoothed_field))<br>
+ else:<br>
+ field_name = ("deposit", "%s_%s_smoothed_%s" % (ptype, kernel_name,<br>
+ smoothed_field))<br>
field_units = registry[ptype, smoothed_field].units<br>
def _vol_weight(field, data):<br>
pos = data[ptype, coord_name]<br>
@@ -789,7 +793,8 @@<br>
# volume_weighted smooth operations return lists of length 1.<br>
rv = data.smooth(pos, [mass, hsml, dens, quan],<br>
method="volume_weighted",<br>
- create_octree=True)[0]<br>
+ create_octree=True,<br>
+ kernel_name=kernel_name)[0]<br>
rv[np.isnan(rv)] = 0.0<br>
# Now some quick unit conversions.<br>
rv = data.apply_units(rv, field_units)<br>
@@ -816,8 +821,12 @@<br>
units = "code_length")<br>
return [field_name]<br>
<br>
-def add_density_kernel(ptype, coord_name, mass_name, registry, nneighbors = 64):<br>
- field_name = (ptype, "smoothed_density")<br>
+def add_density_kernel(ptype, coord_name, mass_name, registry, nneighbors = 64,<br>
+ kernel_name = 'cubic'):<br>
+ if kernel_name == 'cubic':<br>
+ field_name = (ptype, "smoothed_density")<br>
+ else:<br>
+ field_name = (ptype, "%s_smoothed_density" % (kernel_name))<br>
field_units = registry[ptype, mass_name].units<br>
def _nth_neighbor(field, data):<br>
pos = data[ptype, coord_name]<br>
@@ -827,7 +836,8 @@<br>
densities = mass * 0.0<br>
data.particle_operation(pos, [mass, densities],<br>
method="density",<br>
- nneighbors = nneighbors)<br>
+ nneighbors = nneighbors,<br>
+ kernel_name = kernel_name)<br>
ones = pos.prod(axis=1) # Get us in code_length**3<br>
ones[:] = 1.0<br>
densities /= ones<br>
<br>
diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/frontends/artio/data_structures.py<br>
--- a/yt/frontends/artio/data_structures.py<br>
+++ b/yt/frontends/artio/data_structures.py<br>
@@ -133,7 +133,8 @@<br>
tr = dict((field, v) for field, v in zip(fields, tr))<br>
return tr<br>
<br>
- def deposit(self, positions, fields = None, method = None):<br>
+ def deposit(self, positions, fields = None, method = None,<br>
+ kernel_name = 'cubic'):<br>
# Here we perform our particle deposition.<br>
if fields is None: fields = []<br>
cls = getattr(particle_deposit, "deposit_%s" % method, None)<br>
@@ -141,7 +142,8 @@<br>
raise YTParticleDepositionNotImplemented(method)<br>
nz = <a href="http://self.nz" rel="noreferrer" target="_blank">self.nz</a><br>
nvals = (nz, nz, nz, self.ires.size)<br>
- op = cls(nvals) # We allocate number of zones, not number of octs<br>
+ # We allocate number of zones, not number of octs<br>
+ op = cls(nvals, kernel_name)<br>
op.initialize()<br>
mylog.debug("Depositing %s (%s^3) particles into %s Root Mesh",<br>
positions.shape[0], positions.shape[0]**0.3333333, nvals[-1])<br>
<br>
diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/frontends/artio/setup.py<br>
--- a/yt/frontends/artio/setup.py<br>
+++ b/yt/frontends/artio/setup.py<br>
@@ -19,7 +19,8 @@<br>
depends=artio_sources +<br>
["yt/utilities/lib/fp_utils.pxd",<br>
"yt/geometry/oct_container.pxd",<br>
- "yt/geometry/selection_routines.pxd"])<br>
+ "yt/geometry/selection_routines.pxd",<br>
+ "yt/geometry/particle_deposit.pxd"])<br>
config.make_config_py() # installs __config__.py<br>
#config.make_svn_version_py()<br>
return config<br>
<br>
diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/geometry/particle_deposit.pxd<br>
--- a/yt/geometry/particle_deposit.pxd<br>
+++ b/yt/geometry/particle_deposit.pxd<br>
@@ -38,7 +38,7 @@<br>
# Standard SPH kernel for use with the Grid method #<br>
####################################################<br>
<br>
-cdef inline np.float64_t sph_kernel(np.float64_t x) nogil:<br>
+cdef inline np.float64_t sph_kernel_cubic(np.float64_t x) nogil:<br>
cdef np.float64_t kernel<br>
if x <= 0.5:<br>
kernel = 1.-6.*x*x*(1.-x)<br>
@@ -48,8 +48,55 @@<br>
kernel = 0.<br>
return kernel<br>
<br>
+########################################################<br>
+# Alternative SPH kernels for use with the Grid method #<br>
+########################################################<br>
+<br>
+# quartic spline<br>
+cdef inline np.float64_t sph_kernel_quartic(np.float64_t x):<br>
+ cdef np.float64_t kernel<br>
+ cdef np.float64_t C = 5.**6/512/np.pi<br>
+ if x < 1:<br>
+ kernel = (1.-x)**4<br>
+ if x < 3./5:<br>
+ kernel -= 5*(3./5-x)**4<br>
+ if x < 1./5:<br>
+ kernel += 10*(1./5-x)**4<br>
+ else:<br>
+ kernel = 0.<br>
+ return kernel * C<br>
+<br>
+# quintic spline<br>
+cdef inline np.float64_t sph_kernel_quintic(np.float64_t x):<br>
+ cdef np.float64_t kernel<br>
+ cdef np.float64_t C = 3.**7/40/np.pi<br>
+ if x < 1:<br>
+ kernel = (1.-x)**5<br>
+ if x < 2./3:<br>
+ kernel -= 6*(2./3-x)**5<br>
+ if x < 1./3:<br>
+ kernel += 15*(1./3-x)**5<br>
+ else:<br>
+ kernel = 0.<br>
+ return kernel * C<br>
+<br>
+# I don't know the way to use a dict in a cdef class.<br>
+# So in order to mimic a registry functionality,<br>
+# I manually created a function to lookup the kernel functions.<br>
+ctypedef np.float64_t (*kernel_func) (np.float64_t)<br>
+cdef inline kernel_func get_kernel_func(str kernel_name):<br>
+ if kernel_name == 'cubic':<br>
+ return sph_kernel_cubic<br>
+ elif kernel_name == 'quartic':<br>
+ return sph_kernel_quartic<br>
+ elif kernel_name == 'quintic':<br>
+ return sph_kernel_quintic<br>
+ else:<br>
+ raise NotImplementedError<br>
+<br>
cdef class ParticleDepositOperation:<br>
# We assume each will allocate and define their own temporary storage<br>
+ cdef kernel_func sph_kernel<br>
cdef public object nvals<br>
cdef public int update_values<br>
cdef void process(self, int dim[3], np.float64_t left_edge[3],<br>
<br>
diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/geometry/particle_deposit.pyx<br>
--- a/yt/geometry/particle_deposit.pyx<br>
+++ b/yt/geometry/particle_deposit.pyx<br>
@@ -25,9 +25,10 @@<br>
OctreeContainer, OctInfo<br>
<br>
cdef class ParticleDepositOperation:<br>
- def __init__(self, nvals):<br>
+ def __init__(self, nvals, kernel_name):<br>
self.nvals = nvals<br>
self.update_values = 0 # This is the default<br>
+ self.sph_kernel = get_kernel_func(kernel_name)<br>
<br>
def initialize(self, *args):<br>
raise NotImplementedError<br>
@@ -227,7 +228,7 @@<br>
dist = idist[0] + idist[1] + idist[2]<br>
# Calculate distance in multiples of the smoothing length<br>
dist = sqrt(dist) / fields[0]<br>
- self.temp[gind(i,j,k,dim) + offset] = sph_kernel(dist)<br>
+ self.temp[gind(i,j,k,dim) + offset] = self.sph_kernel(dist)<br>
kernel_sum += self.temp[gind(i,j,k,dim) + offset]<br>
# Having found the kernel, deposit accordingly into gdata<br>
for i from ib0[0] <= i <= ib1[0]:<br>
@@ -493,4 +494,3 @@<br>
return self.onnfield<br>
<br>
deposit_nearest = NNParticleField<br>
-<br>
<br>
diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/geometry/particle_smooth.pxd<br>
--- a/yt/geometry/particle_smooth.pxd<br>
+++ b/yt/geometry/particle_smooth.pxd<br>
@@ -22,7 +22,7 @@<br>
<br>
from fp_utils cimport *<br>
from oct_container cimport Oct, OctAllocationContainer, OctreeContainer<br>
-from .particle_deposit cimport sph_kernel, gind<br>
+from .particle_deposit cimport kernel_func, get_kernel_func, gind<br>
<br>
cdef extern from "platform_dep.h":<br>
void *alloca(int)<br>
@@ -34,6 +34,7 @@<br>
<br>
cdef class ParticleSmoothOperation:<br>
# We assume each will allocate and define their own temporary storage<br>
+ cdef kernel_func sph_kernel<br>
cdef public object nvals<br>
cdef np.float64_t DW[3]<br>
cdef int nfields<br>
<br>
diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/geometry/particle_smooth.pyx<br>
--- a/yt/geometry/particle_smooth.pyx<br>
+++ b/yt/geometry/particle_smooth.pyx<br>
@@ -74,7 +74,7 @@<br>
opos[2] = ipos[2]<br>
<br>
cdef class ParticleSmoothOperation:<br>
- def __init__(self, nvals, nfields, max_neighbors):<br>
+ def __init__(self, nvals, nfields, max_neighbors, kernel_name):<br>
# This is the set of cells, in grids, blocks or octs, we are handling.<br>
cdef int i<br>
self.nvals = nvals<br>
@@ -83,6 +83,7 @@<br>
self.neighbors = <NeighborList *> malloc(<br>
sizeof(NeighborList) * self.maxn)<br>
self.neighbor_reset()<br>
+ self.sph_kernel = get_kernel_func(kernel_name)<br>
<br>
def initialize(self, *args):<br>
raise NotImplementedError<br>
@@ -630,7 +631,7 @@<br>
# Usually this density has been computed<br>
dens = fields[2][pn]<br>
if dens == 0.0: continue<br>
- weight = mass * sph_kernel(sqrt(r2) / hsml) / dens<br>
+ weight = mass * self.sph_kernel(sqrt(r2) / hsml) / dens<br>
# Mass of the particle times the value<br>
for fi in range(self.nfields - 3):<br>
val = fields[fi + 3][pn]<br>
@@ -756,7 +757,7 @@<br>
for pn in range(self.curn):<br>
mass = fields[0][self.neighbors[pn].pn]<br>
r2 = self.neighbors[pn].r2<br>
- lw = sph_kernel(sqrt(r2) / hsml)<br>
+ lw = self.sph_kernel(sqrt(r2) / hsml)<br>
dens += mass * lw<br>
weight = (4.0/3.0) * 3.1415926 * hsml**3<br>
fields[1][offset] = dens/weight<br>
<br>
Repository URL: <a href="https://bitbucket.org/yt_analysis/yt/" rel="noreferrer" target="_blank">https://bitbucket.org/yt_analysis/yt/</a><br>
<br>
--<br>
<br>
This is a commit notification from <a href="http://bitbucket.org" rel="noreferrer" target="_blank">bitbucket.org</a>. You are receiving<br>
this because you have the service enabled, addressing the recipient of<br>
this email.<br>
_______________________________________________<br>
yt-svn mailing list<br>
<a href="mailto:yt-svn@lists.spacepope.org">yt-svn@lists.spacepope.org</a><br>
<a href="http://lists.spacepope.org/listinfo.cgi/yt-svn-spacepope.org" rel="noreferrer" target="_blank">http://lists.spacepope.org/listinfo.cgi/yt-svn-spacepope.org</a><br>
</div><br></div></div>