Hey Ben, 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. -Nathan ---------- Forwarded message ---------- From: <commits-noreply@bitbucket.org> Date: Wed, Aug 19, 2015 at 12:09 PM Subject: [yt-svn] commit/yt: bwkeller: Merged in qobilidop/yt (pull request #1670) To: yt-svn@lists.spacepope.org 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. _______________________________________________ yt-svn mailing list yt-svn@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-svn-spacepope.org