Commits

Matthew Turk committed 9a8cdca

This fixes a missing factor of dds[i] in calculating smoothing kernel radius.

  • Participants
  • Parent commits bc401ef

Comments (0)

Files changed (1)

File yt/frontends/sph/smoothing_kernel.pyx

 import numpy as np
 from libc.stdlib cimport free, malloc
 from libc.math cimport sqrt
+from fp_utils cimport iclip
 cimport cython
 
 #@cython.boundscheck(False)
 #@cython.wraparound(False)
 @cython.cdivision(True)
-def grid(input_fields, output_grids,
-         np.ndarray[np.float64_t, ndim=1] left_edge,
-         np.ndarray[np.float64_t, ndim=1] right_edge,
-         np.ndarray[np.float64_t, ndim=2] ppos,
-         np.ndarray[np.float64_t, ndim=1] hsml):
+def smooth_particles(
+        input_fields, output_grids,
+        np.ndarray[np.float64_t, ndim=1] left_edge,
+        np.ndarray[np.float64_t, ndim=1] right_edge,
+        np.ndarray[np.float64_t, ndim=2] ppos,
+        np.ndarray[np.float64_t, ndim=1] hsml):
 
     cdef np.int64_t ngas
-    cdef np.float64_t dds[3], pos[3], idist[3], kern
+    cdef np.float64_t dds[3], sdds[3], pos[3], idist[3], kern
     cdef int p, i, j, k, d, ind[3], ib0[3], ib1[3], dims[3]
     cdef int nf, half_len, skip
     cdef np.float64_t dist
     for i in range(3):
         dims[i] = output_grids[0].shape[i]
         dds[i] = (right_edge[i] - left_edge[i])/dims[i]
+        sdds[i] = dds[i] * dds[i]
 
     cdef np.float64_t *kernel_sum = \
         <np.float64_t *>malloc(ngas * sizeof(np.float64_t))
         for i in range(3):
             pos[i] = ppos[p, i]
             ind[i] = <int>((pos[i] - left_edge[i]) / dds[i])
-            half_len = <int>(hsml[p]/dds[i])
+            half_len = <int>(hsml[p]/dds[i]) + 1
             ib0[i] = ind[i] - half_len
             ib1[i] = ind[i] + half_len
             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
-
         for i from ib0[0] <= i <= ib1[0]:
-            idist[0] = (ind[0] - i) * (ind[0] - i) * dds[0]
+            idist[0] = (ind[0] - i) * (ind[0] - i) * sdds[0]
             for j from ib0[1] <= j <= ib1[1]:
-                idist[1] = (ind[1] - j) * (ind[1] - j) * dds[1] 
-                idist[1] += idist[0]
+                idist[1] = (ind[1] - j) * (ind[1] - j) * sdds[1] 
                 for k from ib0[2] <= k <= ib1[2]:
-                    idist[2] = (ind[2] - k) * (ind[2] - k) * dds[2]
-                    idist[2] += idist[1]
-                    dist = sqrt(idist[2]) / hsml[p]
+                    idist[2] = (ind[2] - k) * (ind[2] - k) * sdds[2]
+                    dist = idist[0] + idist[1] + idist[2]
+                    dist = sqrt(dist) / hsml[p]
                     kernel_sum[p] += sph_kernel(dist)
-
         for i from ib0[0] <= i <= ib1[0]:
-            idist[0] = (ind[0] - i) * (ind[0] - i) * dds[0]
+            idist[0] = (ind[0] - i) * (ind[0] - i) * sdds[0]
             for j from ib0[1] <= j <= ib1[1]:
-                idist[1] = (ind[1] - j) * (ind[1] - j) * dds[1] 
-                idist[1] += idist[0]
+                idist[1] = (ind[1] - j) * (ind[1] - j) * sdds[1] 
                 for k from ib0[2] <= k <= ib1[2]:
-                    idist[2] = (ind[2] - k) * (ind[2] - k) * dds[2]
-                    idist[2] += idist[1]
-                    dist = sqrt(idist[2]) / hsml[p]
+                    idist[2] = (ind[2] - k) * (ind[2] - k) * sdds[2]
+                    dist = idist[0] + idist[1] + idist[2]
+                    dist = sqrt(dist) / hsml[p]
                     kern = sph_kernel(dist)
                     gi = ((i * dims[1] + j) * dims[2]) + k
                     for d in range(nf):
                         gdata[d][gi] += pdata[d][p] * kern / kernel_sum[p]
-
     free(kernel_sum)
     free(gdata)
     free(pdata)