[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