From 0bf0d6d9552ccfdba0db499746f0631a7c3d7b4e Mon Sep 17 00:00:00 2001 From: Prabhu Ramachandran Date: Tue, 26 Apr 2011 22:52:00 +0530 Subject: [PATCH] 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 | 50 +++++++++++++++++++++++++ source/pysph/base/kernels.pyx | 49 +++++++++++++++++++++++- source/pysph/solver/cl_common.h | 2 +- source/pysph/sph/funcs/density_funcs.cl | 14 +------ 4 files changed, 100 insertions(+), 15 deletions(-) diff --git a/source/pysph/base/kernels.h b/source/pysph/base/kernels.h index 41b0b4b..ef5f97d 100644 --- a/source/pysph/base/kernels.h +++ b/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 + diff --git a/source/pysph/base/kernels.pyx b/source/pysph/base/kernels.pyx index debb0f7..f6290b3 100755 --- a/source/pysph/base/kernels.pyx +++ b/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 diff --git a/source/pysph/solver/cl_common.h b/source/pysph/solver/cl_common.h index a106eb5..2823604 100644 --- a/source/pysph/solver/cl_common.h +++ b/source/pysph/solver/cl_common.h @@ -58,4 +58,4 @@ unsigned int get_gid(unsigned int work_dim) } -#endif CL_COMMON +#endif // CL_COMMON diff --git a/source/pysph/sph/funcs/density_funcs.cl b/source/pysph/sph/funcs/density_funcs.cl index aa51bcb..95d705a 100644 --- a/source/pysph/sph/funcs/density_funcs.cl +++ b/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