Commits

Matthew Turk  committed 1b37645 Merge

Merging in the basics of particle deposition

  • Participants
  • Parent commits 77b261f, fd90ac9

Comments (0)

Files changed (11)

File yt/data_objects/data_containers.py

                         if f not in fields_to_generate:
                             fields_to_generate.append(f)
 
+    def deposit(self, positions, fields, op):
+        assert(self._current_chunk.chunk_type == "spatial")
+        fields = ensure_list(fields)
+        self.hierarchy._deposit_particle_fields(self, positions, fields, op)
+
     @contextmanager
     def _field_lock(self):
         self._locked = True

File yt/data_objects/derived_quantities.py

File contents unchanged.

File yt/frontends/art/data_structures.py

File contents unchanged.

File yt/frontends/art/io.py

File contents unchanged.

File yt/frontends/sph/smoothing_kernel.pyx

     for p in range(ngas):
         kernel_sum[p] = 0.0
         skip = 0
+        # Find the # of cells of the kernel
         for i in range(3):
             pos[i] = ppos[p, i]
+            # Get particle root grid integer index
             ind[i] = <int>((pos[i] - left_edge[i]) / dds[i])
+            # How many root grid cells does the smoothing length span + 1
             half_len = <int>(hsml[p]/dds[i]) + 1
+            # Left and right integer indices of the smoothing range
+            # If smoothing len is small could be inside the same bin
             ib0[i] = ind[i] - half_len
             ib1[i] = ind[i] + half_len
             #pos[i] = ppos[p, i] - left_edge[i]
             #ind[i] = <int>(pos[i] / dds[i])
             #ib0[i] = <int>((pos[i] - hsml[i]) / dds[i]) - 1
             #ib1[i] = <int>((pos[i] + hsml[i]) / dds[i]) + 1
+            # Skip if outside out root grid
             if ib0[i] >= dims[i] or ib1[i] < 0:
                 skip = 1
             ib0[i] = iclip(ib0[i], 0, dims[i] - 1)
             ib1[i] = iclip(ib1[i], 0, dims[i] - 1)
         if skip == 1: continue
+        # Having found the kernel shape, calculate the kernel weight
         for i from ib0[0] <= i <= ib1[0]:
             idist[0] = (ind[0] - i) * (ind[0] - i) * sdds[0]
             for j from ib0[1] <= j <= ib1[1]:
                 for k from ib0[2] <= k <= ib1[2]:
                     idist[2] = (ind[2] - k) * (ind[2] - k) * sdds[2]
                     dist = idist[0] + idist[1] + idist[2]
+                    # Calculate distance in multiples of the smoothing length
                     dist = sqrt(dist) / hsml[p]
+                    # Kernel is 3D but save the elements in a 1D array
                     gi = ((i * dims[1] + j) * dims[2]) + k
                     pdist[gi] = sph_kernel(dist)
+                    # Save sum to normalize later
                     kernel_sum[p] += pdist[gi]
+        # Having found the kernel, deposit accordingly into gdata
         for i from ib0[0] <= i <= ib1[0]:
             for j from ib0[1] <= j <= ib1[1]:
                 for k from ib0[2] <= k <= ib1[2]:

File yt/geometry/fake_octree.pyx

+"""
+Make a fake octree, deposit particle at every leaf
+
+Author: Christopher Moody <chris.e.moody@gmail.com>
+Affiliation: UC Santa Cruz
+Author: Matthew Turk <matthewturk@gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2013 Matthew Turk.  All Rights Reserved.
+
+  This file is part of yt.
+
+  yt is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from libc.stdlib cimport malloc, free, rand, RAND_MAX
+cimport numpy as np
+import numpy as np
+cimport cython
+
+from oct_container cimport Oct, RAMSESOctreeContainer
+
+# Create a balanced octree by a random walk that recursively
+# subdivides
+def create_fake_octree(RAMSESOctreeContainer oct_handler,
+                       long max_noct,
+                       long max_level,
+                       np.ndarray[np.int32_t, ndim=1] ndd,
+                       np.ndarray[np.float64_t, ndim=1] dle,
+                       np.ndarray[np.float64_t, ndim=1] dre,
+                       float fsubdivide):
+    cdef int[3] dd #hold the octant index
+    cdef int[3] ind #hold the octant index
+    cdef long i
+    cdef long cur_leaf = 0
+    cdef np.ndarray[np.uint8_t, ndim=2] mask
+    for i in range(3):
+        ind[i] = 0
+        dd[i] = ndd[i]
+    oct_handler.allocate_domains([max_noct])
+    parent = oct_handler.next_root(1, ind)
+    parent.domain = 1
+    cur_leaf = 8 #we've added one parent...
+    mask = np.ones((max_noct,8),dtype='uint8')
+    while oct_handler.domains[0].n_assigned < max_noct:
+        print "root: nocts ", oct_handler.domains[0].n_assigned
+        cur_leaf = subdivide(oct_handler, parent, ind, dd, cur_leaf, 0,
+                             max_noct, max_level, fsubdivide, mask)
+    return cur_leaf
+                             
+
+cdef long subdivide(RAMSESOctreeContainer oct_handler, 
+                    Oct *parent,
+                    int ind[3], int dd[3], 
+                    long cur_leaf, long cur_level, 
+                    long max_noct, long max_level, float fsubdivide,
+                    np.ndarray[np.uint8_t, ndim=2] mask):
+    print "child", parent.ind, ind[0], ind[1], ind[2], cur_leaf, cur_level
+    cdef int ddr[3]
+    cdef long i,j,k
+    cdef float rf #random float from 0-1
+    if cur_level >= max_level: 
+        return cur_leaf
+    if oct_handler.domains[0].n_assigned >= max_noct:
+        return cur_leaf
+    for i in range(3):
+        ind[i] = <int> ((rand() * 1.0 / RAND_MAX) * dd[i])
+        ddr[i] = 2
+    rf = rand() * 1.0 / RAND_MAX
+    if rf > fsubdivide:
+        if parent.children[ind[0]][ind[1]][ind[2]] == NULL:
+            cur_leaf += 7 
+        oct = oct_handler.next_child(1, ind, parent)
+        oct.domain = 1
+        cur_leaf = subdivide(oct_handler, oct, ind, ddr, cur_leaf, 
+                             cur_level + 1, max_noct, max_level, 
+                             fsubdivide, mask)
+    return cur_leaf

File yt/geometry/oct_container.pxd

     cdef np.float64_t DLE[3], DRE[3]
     cdef public int nocts
     cdef public int max_domain
-    cdef Oct* get(self, ppos)
+    cdef Oct* get(self, np.float64_t ppos[3], int *ii = ?)
     cdef void neighbors(self, Oct *, Oct **)
     cdef void oct_bounds(self, Oct *, np.float64_t *, np.float64_t *)
 

File yt/geometry/oct_container.pyx

     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
-    cdef Oct *get(self, ppos):
+    cdef Oct *get(self, np.float64_t ppos[3], int *ii = NULL):
         #Given a floating point position, retrieve the most
         #refined oct at that time
         cdef np.int64_t ind[3]
                     ind[i] = 1
                     cp[i] += dds[i]/2.0
             cur = cur.children[ind[0]][ind[1]][ind[2]]
+        if ii != NULL: return cur
+        for i in range(3):
+            if cp[i] > pp[i]:
+                ind[i] = 0
+            else:
+                ind[i] = 1
+        ii[0] = ((ind[2]*2)+ind[1])*2+ind[0]
         return cur
 
     @cython.boundscheck(False)
                 count[o.domain - 1] += mask[o.local_ind,i]
         return count
 
+    @cython.boundscheck(True)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    def count_leaves(self, np.ndarray[np.uint8_t, ndim=2, cast=True] mask):
+        # Modified to work when not all octs are assigned
+        cdef int i, j, k, ii
+        cdef np.int64_t oi
+        # pos here is CELL center, not OCT center.
+        cdef np.float64_t pos[3]
+        cdef int n = mask.shape[0]
+        cdef np.ndarray[np.int64_t, ndim=1] count
+        count = np.zeros(self.max_domain, 'int64')
+        # 
+        cur = self.cont
+        for oi in range(n):
+            if oi - cur.offset >= cur.n_assigned:
+                cur = cur.next
+                if cur == NULL:
+                    break
+            o = &cur.my_octs[oi - cur.offset]
+            # skip if unassigned
+            if o == NULL:
+                continue
+            if o.domain == -1: 
+                continue
+            for i in range(2):
+                for j in range(2):
+                    for k in range(2):
+                        if o.children[i][j][k] == NULL:
+                            ii = ((k*2)+j)*2+i
+                            count[o.domain - 1] += mask[o.local_ind,ii]
+        return count
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
-    def get_neighbor_boundaries(self, ppos):
+    def get_neighbor_boundaries(self, oppos):
+        cdef int i, ii
+        cdef np.float64_t ppos[3]
+        for i in range(3):
+            ppos[i] = oppos[i]
         cdef Oct *main = self.get(ppos)
         cdef Oct* neighbors[27]
         self.neighbors(main, neighbors)
         cdef np.ndarray[np.float64_t, ndim=2] bounds
         cdef np.float64_t corner[3], size[3]
         bounds = np.zeros((27,6), dtype="float64")
-        cdef int i, ii
         tnp = 0
         for i in range(27):
             self.oct_bounds(neighbors[i], corner, size)
                 m2[o.local_ind, i] = mask[o.local_ind, i]
         return m2
 
-    def check(self, int curdom):
+    def check(self, int curdom, int print_all = 0):
         cdef int dind, pi
         cdef Oct oct
         cdef OctAllocationContainer *cont = self.domains[curdom - 1]
         cdef int unassigned = 0
         for pi in range(cont.n_assigned):
             oct = cont.my_octs[pi]
+            if print_all==1:
+                print pi, oct.level, oct.domain,
+                print oct.pos[0],oct.pos[1],oct.pos[2]
             for i in range(2):
                 for j in range(2):
                     for k in range(2):
 
 cdef class ARTOctreeContainer(RAMSESOctreeContainer):
     #this class is specifically for the NMSU ART
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    def deposit_particle_cumsum(self,
+                                np.ndarray[np.float64_t, ndim=2] ppos, 
+                                np.ndarray[np.float64_t, ndim=1] pdata,
+                                np.ndarray[np.float64_t, ndim=1] mask,
+                                np.ndarray[np.float64_t, ndim=1] dest,
+                                fields, int domain):
+        cdef Oct *o
+        cdef OctAllocationContainer *dom = self.domains[domain - 1]
+        cdef np.float64_t pos[3]
+        cdef int ii
+        cdef int no = ppos.shape[0]
+        for n in range(no):
+            for j in range(3):
+                pos[j] = ppos[n,j]
+            o = self.get(pos, &ii) 
+            if mask[o.local_ind,ii]==0: continue
+            dest[o.ind+ii] += pdata[n]
+        return dest
+
     @cython.boundscheck(True)
     @cython.wraparound(False)
     @cython.cdivision(True)
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
-    def count_neighbor_particles(self, ppos):
+    def count_neighbor_particles(self, oppos):
         #How many particles are in my neighborhood
+        cdef int i, ni, dl, tnp
+        cdef np.float64_t ppos[3]
+        for i in range(3):
+            ppos[i] = oppos[i]
         cdef Oct *main = self.get(ppos)
         cdef Oct* neighbors[27]
         self.neighbors(main, neighbors)
-        cdef int i, ni, dl, tnp
         tnp = 0
         for i in range(27):
             if neighbors[i].sd != NULL:

File yt/geometry/oct_deposit.pyx

+"""
+Particle Deposition onto Octs
+
+Author: Christopher Moody <chris.e.moody@gmail.com>
+Affiliation: UC Santa Cruz
+Author: Matthew Turk <matthewturk@gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2013 Matthew Turk.  All Rights Reserved.
+
+  This file is part of yt.
+
+  yt is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from libc.stdlib cimport malloc, free
+cimport numpy as np
+import numpy as np
+cimport cython
+
+from oct_container cimport Oct, OctAllocationContainer, OctreeContainer
+
+# Mode functions
+ctypedef np.float64_t (*type_opt)(np.float64_t, np.float64_t)
+cdef np.float64_t opt_count(np.float64_t pdata,
+                            np.float64_t weight,
+                            np.int64_t index,
+                            np.ndarray[np.float64_t, ndim=2] data_out, 
+                            np.ndarray[np.float64_t, ndim=2] data_in):
+    data_out[index] += 1.0
+
+cdef np.float64_t opt_sum(np.float64_t pdata,
+                            np.float64_t weight,
+                            np.int64_t index,
+                            np.ndarray[np.float64_t, ndim=2] data_out, 
+                            np.ndarray[np.float64_t, ndim=2] data_in):
+    data_out[index] += pdata 
+
+cdef np.float64_t opt_diff(np.float64_t pdata,
+                            np.float64_t weight,
+                            np.int64_t index,
+                            np.ndarray[np.float64_t, ndim=2] data_out, 
+                            np.ndarray[np.float64_t, ndim=2] data_in):
+    data_out[index] += (data_in[index] - pdata) 
+
+cdef np.float64_t opt_wcount(np.float64_t pdata,
+                            np.float64_t weight,
+                            np.int64_t index,
+                            np.ndarray[np.float64_t, ndim=2] data_out, 
+                            np.ndarray[np.float64_t, ndim=2] data_in):
+    data_out[index] += weight
+
+cdef np.float64_t opt_wsum(np.float64_t pdata,
+                            np.float64_t weight,
+                            np.int64_t index,
+                            np.ndarray[np.float64_t, ndim=2] data_out, 
+                            np.ndarray[np.float64_t, ndim=2] data_in):
+    data_out[index] += pdata * weight
+
+cdef np.float64_t opt_wdiff(np.float64_t pdata,
+                            np.float64_t weight,
+                            np.int64_t index,
+                            np.ndarray[np.float64_t, ndim=2] data_out, 
+                            np.ndarray[np.float64_t, ndim=2] data_in):
+    data_out[index] += (data_in[index] - pdata) * weight
+
+# Selection functions
+ctypedef NOTSURE (*type_sel)(OctreeContainer, 
+                                np.ndarray[np.float64_t, ndim=1],
+                                np.float64_t)
+cdef NOTSURE select_nearest(OctreeContainer oct_handler,
+                            np.ndarray[np.float64_t, ndim=1] pos,
+                            np.float64_t radius):
+    #return only the nearest oct
+    pass
+
+
+cdef NOTSURE select_radius(OctreeContainer oct_handler,
+                            np.ndarray[np.float64_t, ndim=1] pos,
+                            np.float64_t radius):
+    #return a list of octs within the radius
+    pass
+    
+
+# Kernel functions
+ctypedef np.float64_t (*type_ker)(np.float64_t)
+cdef np.float64_t kernel_sph(np.float64_t x) nogil:
+    cdef np.float64_t kernel
+    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
+
+cdef np.float64_t kernel_null(np.float64_t x) nogil: return 0.0
+
+cdef deposit(OctreeContainer oct_handler, 
+        np.ndarray[np.float64_t, ndim=2] ppos, #positions,columns are x,y,z
+        np.ndarray[np.float64_t, ndim=2] pd, # particle fields
+        np.ndarray[np.float64_t, ndim=1] pr, # particle radius
+        np.ndarray[np.float64_t, ndim=2] data_in, #used to calc diff, same shape as data_out
+        np.ndarray[np.float64_t, ndim=2] data_out, #write deposited here
+        mode='count', selection='nearest', kernel='null'):
+    cdef type_opt fopt
+    cdef type_sel fsel
+    cdef type_ker fker
+    cdef long pi #particle index
+    cdef long nocts #number of octs in selection
+    cdef Oct oct 
+    cdef np.float64_t w
+    # Can we do this with dicts?
+    # Setup the function pointers
+    if mode == 'count':
+        fopt = opt_count
+    elif mode == 'sum':
+        fopt = opt_sum
+    elif mode == 'diff':
+        fopt = opt_diff
+    if mode == 'wcount':
+        fopt = opt_count
+    elif mode == 'wsum':
+        fopt = opt_sum
+    elif mode == 'wdiff':
+        fopt = opt_diff
+    if selection == 'nearest':
+        fsel = select_nearest
+    elif selection == 'radius':
+        fsel = select_radius
+    if kernel == 'null':
+        fker = kernel_null
+    if kernel == 'sph':
+        fker = kernel_sph
+    for pi in range(particles):
+        octs = fsel(oct_handler, ppos[pi], pr[pi])
+        for oct in octs:
+            for cell in oct.cells:
+                w = fker(pr[pi],cell) 
+                weights.append(w)
+        norm = weights.sum()
+        for w, oct in zip(weights, octs):
+            for cell in oct.cells:
+                fopt(pd[pi], w/norm, oct.index, data_in, data_out)
+
+

File yt/geometry/setup.py

                 depends=["yt/utilities/lib/fp_utils.pxd",
                          "yt/geometry/oct_container.pxd",
                          "yt/geometry/selection_routines.pxd"])
+    config.add_extension("fake_octree", 
+                ["yt/geometry/fake_octree.pyx"],
+                include_dirs=["yt/utilities/lib/"],
+                libraries=["m"],
+                depends=["yt/utilities/lib/fp_utils.pxd",
+                         "yt/geometry/oct_container.pxd",
+                         "yt/geometry/selection_routines.pxd"])
     config.make_config_py() # installs __config__.py
     #config.make_svn_version_py()
     return config

File yt/geometry/tests/fake_octree.py

+from yt.geometry.fake_octree import create_fake_octree
+from yt.geometry.oct_container import RAMSESOctreeContainer, ParticleOctreeContainer
+import numpy as np
+
+nocts= 3
+max_level = 12
+dn = 2
+dd = np.ones(3,dtype='i4')*dn
+dle = np.ones(3,dtype='f8')*0.0
+dre = np.ones(3,dtype='f8')
+fsub = 0.25
+domain = 1
+
+oct_handler = RAMSESOctreeContainer(dd,dle,dre)
+leaves = create_fake_octree(oct_handler, nocts, max_level, dd, dle, dre, fsub)
+mask = np.ones((nocts,8),dtype='bool')
+cell_count = nocts*8
+oct_counts = oct_handler.count_levels(max_level, 1, mask)
+level_counts = np.concatenate(([0,],np.cumsum(oct_counts)))
+fc = oct_handler.fcoords(domain,mask,cell_count, level_counts.copy())
+leavesb = oct_handler.count_leaves(mask)
+assert leaves == leavesb
+
+#Now take the fcoords, call them particles and recreate the same octree
+print "particle-based recreate"
+oct_handler2 = ParticleOctreeContainer(dd,dle,dre)
+oct_handler2.allocate_domains([nocts])
+oct_handler2.n_ref = 1 #specifically make a maximum of 1 particle per oct
+oct_handler2.add(fc, 1)
+print "added particles"
+cell_count2 = nocts*8
+oct_counts2 = oct_handler.count_levels(max_level, 1, mask)
+level_counts2 = np.concatenate(([0,],np.cumsum(oct_counts)))
+fc2 = oct_handler.fcoords(domain,mask,cell_count, level_counts.copy())
+leaves2 = oct_handler2.count_leaves(mask)
+assert leaves == leaves2
+
+print "success"