Adding OpenCL support for energy functions.

Modifying tests for pressure functions to include viscosity.

source/pysph/base/kernels.h

@@ -86,6 +86,7 @@ void cubic_spline_gradient(REAL4 pa, REAL4 pb, REAL4* grad,
   grad[0].s0 = val * rab.s0;
   grad[0].s1 = val * rab.s1;
   grad[0].s2 = val * rab.s2;
+  grad[0].s3 = 0.0F;
 
 }
 

source/pysph/sph/funcs/energy_funcs.clt

@@ -0,0 +1,144 @@
+
+
+$EnergyEquationNoVisc
+
+/* Thermal energy equation (eq 10) described in "A shock-captiting SPH
+   scheme based on adaptive kernel estimation" by Leonardo
+   Di. G. Sigalotti et al.
+
+   The equation without artificial viscosity and artificial heat is
+
+   \frac{DU_a}{Dt} = \frac{1}{2}\sum_{b=1}^{N}m_b\left[
+   \left(\frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2}\right)\,(v_a -
+   v_b)\right]\,\nabla_a \cdot W_{ab}
+
+*/
+
+#include "cl_common.h"
+#include "cl_common.cl"
+#include "kernels.h"
+
+__kernel void EnergyEquationNoVisc(%(kernel_args)s)
+{
+    %(workgroup_code)s
+
+    // The term `dest_id` will be suitably defined at this point.
+
+    REAL4 ra = (REAL4)( d_x[dest_id], d_y[dest_id], 
+                        d_z[dest_id], d_h[dest_id] );
+
+    REAL rhoa = d_rho[dest_id];
+    REAL Pa = d_p[dest_id];
+
+    REAL4 va = (REAL4)(d_u[dest_id], d_v[dest_id], d_w[dest_id], 0.0F);
+
+    REAL4 grad;
+    REAL temp;
+
+    %(neighbor_loop_code)s 
+
+       {
+
+        // SPH innermost loop code goes here.  The index `src_id` will
+        // be available and looped over, this index.
+
+        REAL4 rb = (REAL4)(s_x[src_id], s_y[src_id], s_z[src_id], s_h[src_id]);
+        REAL Pb = s_p[src_id];
+        REAL mb = s_m[src_id];
+	REAL rhob = s_rho[src_id];
+
+	REAL4 vb = (REAL4)(s_u[src_id], s_v[src_id], s_w[src_id], 0.0F);
+
+	temp = 0.5F*mb*( Pa/(rhoa*rhoa) + Pb/(rhob*rhob) ) ;
+
+	kernel_gradient(ra, rb, &grad, dim, kernel_type);
+
+	tmpx[dest_id] += temp * ( dot(va-vb, grad) );
+
+       } // neighbor loop
+
+} // __kernel
+
+$EnergyEquationNoVisc
+
+
+$EnergyEquationWithVisc
+
+/* Thermal energy equation (eq 10) described in "A shock-captiting SPH
+   scheme based on adaptive kernel estimation" by Leonardo
+   Di. G. Sigalotti et al.
+
+   The equation without artificial heat is
+
+   \frac{DU_a}{Dt} = \frac{1}{2}\sum_{b=1}^{N}m_b\left[
+   \left(\frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} +
+   \Pi_{ab}\right)\,(v_a - v_b)\right]\,\nabla_a \cdot W_{ab}
+
+*/
+
+#include "cl_common.h"
+#include "cl_common.cl"
+#include "kernels.h"
+
+__kernel void EnergyEquationWithVisc(%(kernel_args)s)
+{
+    %(workgroup_code)s
+
+    // The term `dest_id` will be suitably defined at this point.
+
+    REAL4 ra = (REAL4)( d_x[dest_id], d_y[dest_id], 
+                        d_z[dest_id], d_h[dest_id] );
+
+    REAL4 va = (REAL4)(d_u[dest_id], d_v[dest_id], d_w[dest_id], 0.0F);
+
+    REAL rhoa = d_rho[dest_id];
+    REAL Pa = d_p[dest_id];
+    REAL4 grad;
+    REAL temp, Pi_ab;
+
+    %(neighbor_loop_code)s 
+
+       {
+
+        // SPH innermost loop code goes here.  The index `src_id` will
+        // be available and looped over, this index.
+
+        REAL4 rb = (REAL4)(s_x[src_id], s_y[src_id], s_z[src_id], s_h[src_id]);
+        REAL Pb = s_p[src_id];
+	REAL rhob = s_rho[src_id];
+
+	REAL4 vb = (REAL4)(s_u[src_id], s_v[src_id], s_w[src_id], 0.0F);
+
+	REAL dot_product = dot( (va-vb), (ra-rb) );
+
+	temp =  Pa/(rhoa*rhoa) + Pb/(rhob*rhob)  ;
+
+	kernel_gradient(ra, rb, &grad, dim, kernel_type);
+
+	Pi_ab = 0.0F;	
+	if ( dot_product < 0.0F )
+	  {
+	    REAL cab = 0.5F * ( d_cs[dest_id] + s_cs[src_id] );
+	    REAL rhoab = 0.5F * (rhoa + rhob);
+
+	    REAL hab = 0.5F * ( d_h[dest_id] + s_h[src_id] );
+	    REAL mu = dot_product*hab;
+	    REAL norm = length(ra-rb);
+
+	    mu /= ( norm*norm + eta*eta*hab*hab );
+
+	    Pi_ab = -alpha*cab*mu + beta*mu*mu;
+	    Pi_ab /= rhoab;
+
+	  } //if
+
+	temp += Pi_ab;
+	temp *= 0.5F*s_m[src_id];
+
+	tmpx[dest_id] += temp * dot( va-vb, grad );
+
+       } //neighbor_loop_code
+
+} //__kernel
+
+$EnergyEquationWithVisc

source/pysph/sph/funcs/energy_funcs.pyx

@@ -2,6 +2,8 @@
 from pysph.base.point cimport cPoint, cPoint_dot, cPoint_new, cPoint_sub,\
      cPoint_norm
 
+from pysph.solver.cl_utils import get_real
+
 cdef extern from "math.h":
     double sqrt(double)
     double fabs(double)
@@ -31,6 +33,9 @@ cdef class EnergyEquationNoVisc(SPHFunctionParticle):
 
         self.src_reads.extend( ['u','v','w','p'] )
         self.dst_reads.extend( ['u','v','w','p'] )
+
+    def _set_extra_cl_args(self):
+        pass
 
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
                        KernelBase kernel, double *nr):
@@ -89,7 +94,14 @@ cdef class EnergyEquationNoVisc(SPHFunctionParticle):
         tmp = 0.5*mb*(pa/(rhoa*rhoa) + pb/(rhob*rhob))
 
         nr[0] += tmp*dot
-##############################################################################
+
+    def cl_eval(self, object queue, object context, object kernel):
+
+        self.set_cl_kernel_args()        
+
+        self.cl_program.EnergyEquationNoVisc(
+            queue, self.global_sizes, self.local_sizes, *self.cl_args).wait()
+
 
 ##############################################################################
 cdef class EnergyEquationAVisc(SPHFunctionParticle):
@@ -232,12 +244,13 @@ cdef class EnergyEquation(SPHFunctionParticle):
         self.alpha = alpha
         self.beta = beta
         self.gamma = gamma
+        self.eta = eta
 
         self.id = 'energyequation'
         self.tag = "energy"
 
         self.cl_kernel_src_file = "energy_funcs.cl"
-        self.cl_kernel_function_name = "EnergyEquation"
+        self.cl_kernel_function_name = "EnergyEquationWithVisc"
 
     def set_src_dst_reads(self):
         self.src_reads = []
@@ -247,7 +260,20 @@ cdef class EnergyEquation(SPHFunctionParticle):
         self.dst_reads.extend( ['x','y','z','h','rho','tag'] )
 
         self.src_reads.extend( ['u','v','w','p','cs'] )
-        self.dst_reads.extend( ['u','v','w','p','cs'] )        
+        self.dst_reads.extend( ['u','v','w','p','cs'] )
+
+    def _set_extra_cl_args(self):
+        self.cl_args.append( get_real(self.alpha, self.dest.cl_precision) )
+        self.cl_args_name.append( 'REAL const alpha' )
+
+        self.cl_args.append( get_real(self.beta, self.dest.cl_precision) )
+        self.cl_args_name.append( 'REAL const beta' )
+
+        self.cl_args.append( get_real(self.gamma, self.dest.cl_precision) )
+        self.cl_args_name.append( 'REAL const gamma' )
+
+        self.cl_args.append( get_real(self.eta, self.dest.cl_precision) )
+        self.cl_args_name.append( 'REAL const eta' )        
 
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
                        KernelBase kernel, double *nr):    
@@ -332,8 +358,14 @@ cdef class EnergyEquation(SPHFunctionParticle):
         tmp = cPoint_dot(grad, vab) * tmp
 
         nr[0] += 0.5*mb*tmp
-
 
+    def cl_eval(self, object queue, object context, object kernel):
+
+        self.set_cl_kernel_args()        
+
+        self.cl_program.EnergyEquationWithVisc(
+            queue, self.global_sizes, self.local_sizes, *self.cl_args).wait()
+
 ################################################################################
 # `ArtificialHeat` class.
 ################################################################################

source/pysph/sph/funcs/pressure_funcs.clt

@@ -69,7 +69,6 @@ __kernel void MomentumEquation(%(kernel_args)s)
     REAL4 grad;
     REAL temp, Pi_ab;
 
-    Pi_ab = 0.0F;
     %(neighbor_loop_code)s 
 
        {
@@ -80,6 +79,7 @@ __kernel void MomentumEquation(%(kernel_args)s)
         REAL4 rb = (REAL4)(s_x[src_id], s_y[src_id], s_z[src_id], s_h[src_id]);
         REAL Pb = s_p[src_id];
 	REAL rhob = s_rho[src_id];
+	REAL mb = s_m[src_id];
 
 	REAL4 vb = (REAL4)(s_u[src_id], s_v[src_id], s_w[src_id], 0.0F);
 
@@ -88,24 +88,26 @@ __kernel void MomentumEquation(%(kernel_args)s)
 	temp =  Pa/(rhoa*rhoa) + Pb/(rhob*rhob)  ;
 
 	kernel_gradient(ra, rb, &grad, dim, kernel_type);
-
+
+	Pi_ab = 0.0F;	
 	if ( dot_product < 0.0F )
 	  {
 	    REAL cab = 0.5F * ( d_cs[dest_id] + s_cs[src_id] );
 	    REAL rhoab = 0.5F * (rhoa + rhob);
 
 	    REAL hab = 0.5F * ( d_h[dest_id] + s_h[src_id] );
 	    REAL mu = dot_product*hab;
-
-	    mu /= ( length(ra-rb) + eta*eta*hab*hab );
+	    REAL norm2 = length(ra-rb) * length(ra-rb);
+
+	    mu /= ( norm2 + eta*eta*hab*hab );
 
 	    Pi_ab = -alpha*cab*mu + beta*mu*mu;
 	    Pi_ab /= rhoab;
 
 	  } //if
 
 	temp += Pi_ab;
-	temp *= -s_m[src_id];
+	temp = -mb*temp;
 
 	tmpx[dest_id] += temp*grad.x;
 	tmpy[dest_id] += temp*grad.y;