ENH: Returning a correct and unique kernel type from the Cython SPH

kernels.  Using the same values for the OpenCL code and creating generic
functions for kernel function and gradient evaluation in OpenCL using
the kernel type.  This simplifies the OpenCL kernel codes considerably.
Also fixed a small error in the cl_common.h file.

source/pysph/base/kernels.h

@@ -1,7 +1,23 @@
 // C declarations for the SPH smoothing kernels
+#ifndef CL_KERNEL_H
+#define CL_KERNEL_H
 
 #define PI 2.0f*acos(0.0f)
 
+enum KernelType {
+    CUBIC_SPLINE = 1,
+    GAUSSIAN = 2,
+    QUINTIC_SPLINE = 3,
+    WENDLAND_QUINTIC_SPLINE = 4,
+    HARMONIC = 5,
+    M6_SPLINE = 6,
+    W8 = 7,
+    W10 = 8,
+    REPULSIVE = 9,
+    POLY6 = 10
+};
+
+
 REAL cubic_spline_fac(unsigned int dim, REAL h)
 {
   REAL fac = 0.0;
@@ -72,3 +88,37 @@ void cubic_spline_gradient(REAL4 pa, REAL4 pb, REAL4* grad,
   grad[0].s2 = val * rab.s2;
 
 }
+
+REAL kernel_function(REAL4 pa, REAL4 pb, unsigned int dim, int kernel_type)
+{
+    REAL w;
+
+    switch (kernel_type) {
+        case CUBIC_SPLINE:
+            w = cubic_spline_function(pa, pb, dim);
+            break;
+        // FIXME: implement the rest!
+        default:
+            w = cubic_spline_function(pa, pb, dim);
+            break;
+    } // switch (kernel_type).
+
+    return w;
+}
+
+void kernel_gradient(REAL4 pa, REAL4 pb, REAL4* grad, 
+                     unsigned int dim, int kernel_type)
+{
+    switch (kernel_type) {
+        case CUBIC_SPLINE:
+            cubic_spline_gradient(pa, pb, grad, dim);
+            break;
+        // FIXME: implement the rest!
+        default:
+            cubic_spline_gradient(pa, pb, grad, dim);
+            break;
+    } // switch (kernel_type).
+}
+
+#endif // CL_KERNEL_H
+

source/pysph/base/kernels.pyx

@@ -25,6 +25,25 @@ cdef:
     double SQRT_1_PI = 1.0/sqrt(PI)
     double infty = numpy.inf
 
+
+##############################################################################
+# These are simple ints to define the different types of kernels.  This
+# is what is returned by the `KernelBase.get_type` method. The only
+# reason we do this is to unify the numbers used between Cython and
+# OpenCL.
+cdef:
+    int CUBIC_SPLINE = 1
+    int GAUSSIAN = 2
+    int QUINTIC_SPLINE = 3
+    int WENDLAND_QUINTIC_SPLINE = 4
+    int HARMONIC = 5
+    int M6_SPLINE = 6
+    int W8 = 7
+    int W10 = 8
+    int REPULSIVE = 9
+    int POLY6 = 10
+
+
 cdef inline double h_dim(double h, int dim):
     if dim == 1:
         return 1/h
@@ -300,6 +319,9 @@ cdef class Poly6Kernel(KernelBase):
     cpdef double radius(self):
         return 1.0
 
+    def get_type(self):
+        return POLY6 
+
 ##############################################################################
 # DummyKernel` class.
 ##############################################################################
@@ -424,7 +446,8 @@ cdef class CubicSplineKernel(KernelBase):
         return 2
 
     def get_type(self):
-        return 1
+        return CUBIC_SPLINE
+
 ##############################################################################
 
 
@@ -526,6 +549,8 @@ cdef class QuinticSplineKernel(KernelBase):
     cpdef double radius(self):
         return 3
 
+    def get_type(self):
+        return QUINTIC_SPLINE 
 
 ##############################################################################
 #`WendlandQuinticSplineKernel`
@@ -605,6 +630,9 @@ cdef class WendlandQuinticSplineKernel(KernelBase):
     cpdef double radius(self):
         return 2
 
+    def get_type(self):
+        return WENDLAND_QUINTIC_SPLINE 
+
 ##############################################################################
 #`HarmonicKernel`
 ##############################################################################
@@ -697,6 +725,10 @@ cdef class HarmonicKernel(KernelBase):
 
     cpdef double radius(self):
         return 2
+
+    def get_type(self):
+        return HARMONIC
+
 #############################################################################
 
 ##############################################################################
@@ -769,6 +801,10 @@ cdef class M6SplineKernel(KernelBase):
 
     cpdef double radius(self):
         return 2
+
+    def get_type(self):
+        return M6_SPLINE
+
 #############################################################################
 
 ##############################################################################
@@ -826,6 +862,9 @@ cdef class GaussianKernel(KernelBase):
 
     cpdef double radius(self):
         return 3.0
+
+    def get_type(self):
+        return GAUSSIAN 
 ##############################################################################
 
 
@@ -903,6 +942,9 @@ cdef class W8Kernel(KernelBase):
     cpdef double radius(self):
         return 2
 
+    def get_type(self):
+        return W8 
+
 ##############################################################################
 
 
@@ -978,6 +1020,9 @@ cdef class W10Kernel(KernelBase):
     cpdef double radius(self):
         return 2
 
+    def get_type(self):
+        return W10 
+
 
 ################################################################################
 # `RepulsiveBoundaryKernel` class.
@@ -1056,5 +1101,7 @@ cdef class RepulsiveBoundaryKernel(KernelBase):
             raise ValueError
         else: return (1.0/h)
 
+    def get_type(self):
+        return REPULSIVE 
 
 

source/pysph/solver/cl_common.h

@@ -58,4 +58,4 @@ unsigned int get_gid(unsigned int work_dim)
 
 }
 
-#endif CL_COMMON
+#endif // CL_COMMON

source/pysph/sph/funcs/density_funcs.cl

@@ -24,19 +24,7 @@ __kernel void SPHRho(int const kernel_type, int const dim, int const nbrs,
     {
       REAL4 pb = (REAL4)( s_x[i], s_y[i], s_z[i], s_h[i] );
       REAL mb = s_m[i];
-
-      switch (kernel_type)
-	{
-	case 1:
-	  w = cubic_spline_function(pa, pb, dim);
-	  break;
-
-	default:
-	  w = cubic_spline_function(pa, pb, dim);
-	  break;
-
-	} // switch
-
+      w = kernel_function(pa, pb, dim, kernel_type); 
       wmb += w*mb;
 
     } // for i