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