diff --git a/source/pysph/sph/funcs/adke_funcs.pyx b/source/pysph/sph/funcs/adke_funcs.pyx
index 62d9506..400646e 100644
--- a/source/pysph/sph/funcs/adke_funcs.pyx
+++ b/source/pysph/sph/funcs/adke_funcs.pyx
@@ -85,6 +85,14 @@ cdef class ADKESmoothingUpdate(ADKEPilotRho):
         self.cl_kernel_src_file = "adke_funcs.cl"
         self.cl_kernel_function_name = "ADKESmoothingUpdate"
 
+    def set_src_dst_reads(self):
+
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','m','rho'] )
+        self.dst_reads.extend( ['x','y','z','tag'] )
+
     cpdef eval(self, KernelBase kernel, DoubleArray output1,
                DoubleArray output2, DoubleArray output3):
         """ Evaluate the store the results in the output arrays """
@@ -147,12 +155,16 @@ cdef class SPHVelocityDivergence(SPHFunctionParticle):
         self.id = "vdivergence"
         self.tag = "vdivergence"
 
-        self.src_reads.extend( ['u','v','w'] )
-        self.dst_reads.extend( ['u','v','w'] )
-
         self.cl_kernel_src_file = "adke_funcs.cl"
         self.cl_kernel_function_name = "SPHVelocityDivergence"
 
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','u','v','w','m'] )
+        self.dst_reads.extend( ['x','y','z','h','u','v','w','rho','tag'] )
+
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid,
                        KernelBase kernel, double *nr):
         """ Compute the contribution from source_pid on dest_pid.
@@ -166,7 +178,6 @@ cdef class SPHVelocityDivergence(SPHFunctionParticle):
 
 
         """
-
         cdef double ha = self.d_h.data[dest_pid]
         cdef double hb = self.s_h.data[source_pid]
 
@@ -232,6 +243,13 @@ cdef class ADKEConductionCoeffUpdate(SPHVelocityDivergence):
 
         self.cl_kernel_function_name = "ADKEConductionCoeffUpdate"
 
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','u','v','w','m'] )
+        self.dst_reads.extend( ['x','y','z','h','u','v','w','rho','cs','tag'] )
+        
     cpdef eval(self, KernelBase kernel, DoubleArray output1,
                DoubleArray output2, DoubleArray output3):
         """ Evaluate the store the results in the output arrays Note
diff --git a/source/pysph/sph/funcs/basic_funcs.pyx b/source/pysph/sph/funcs/basic_funcs.pyx
index 5aa5036..79cc29b 100644
--- a/source/pysph/sph/funcs/basic_funcs.pyx
+++ b/source/pysph/sph/funcs/basic_funcs.pyx
@@ -33,6 +33,13 @@ cdef class SPH(CSPHFunctionParticle):
         self.cl_kernel_src_file = "basic_funcs.clt"
         self.cl_kernel_function_name = "SPH"
 
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','m','rho',self.prop_name] )
+        self.dst_reads.extend( ['x','y','z','h','tag'] )
+
     cpdef setup_arrays(self):
         """ Setup the arrays required to read data from source and dest. """
 
@@ -129,6 +136,13 @@ cdef class SPHSimpleGradient(SPHFunctionParticle):
         self.cl_kernel_src_file = "basic_funcs.cl"
         self.cl_kernel_function_name = "SPHSimpleGradient"
 
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','m','rho',self.prop_name] )
+        self.dst_reads.extend( ['x','y','z','h','tag'] )
+
     cpdef setup_arrays(self):
         """ Setup the arrays required to read data from source and dest. """
 
@@ -138,9 +152,6 @@ cdef class SPHSimpleGradient(SPHFunctionParticle):
         self.d_prop = self.dest.get_carray(self.prop_name)
         self.s_prop = self.source.get_carray(self.prop_name)
 
-        self.src_reads.append(self.prop_name)
-        self.dst_reads.append(self.prop_name)
-
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid,
                        KernelBase kernel, double *nr):
         """ 
@@ -227,6 +238,13 @@ cdef class SPHGradient(SPHFunctionParticle):
         self.cl_kernel_src_file = "basic_funcs.cl"
         self.cl_kernel_function_name = "SPHGradient"
 
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','m','rho',self.prop_name] )
+        self.dst_reads.extend( ['x','y','z','h','tag'] )
+
     cpdef setup_arrays(self):
         """ Setup the arrays required to read data from source and dest. """
 
@@ -236,9 +254,6 @@ cdef class SPHGradient(SPHFunctionParticle):
         self.d_prop = self.dest.get_carray(self.prop_name)
         self.s_prop = self.source.get_carray(self.prop_name)
 
-        self.src_reads.append(self.prop_name)
-        self.dst_reads.append(self.prop_name)
-
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
                        KernelBase kernel, double *nr):
         """ Perform an SPH interpolation of the property `prop_name` 
@@ -323,6 +338,13 @@ cdef class SPHLaplacian(SPHFunctionParticle):
         self.cl_kernel_src_file = "basic_funcs.cl"
         self.cl_kernel_function_name = "SPHLaplacian"
 
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','m','rho',self.prop_name] )
+        self.dst_reads.extend( ['x','y','z','h','tag',self.prop_name] )        
+
     cpdef setup_arrays(self):
         """ Setup the arrays required to read data from source and dest. """
 
@@ -331,9 +353,6 @@ cdef class SPHLaplacian(SPHFunctionParticle):
         self.s_prop = self.source.get_carray(self.prop_name)
         self.d_prop = self.dest.get_carray(self.prop_name)
 
-        self.src_reads.append(self.prop_name)
-        self.dst_reads.append(self.prop_name)
-
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
                        KernelBase kernel, double *nr):
         """ 
@@ -411,6 +430,9 @@ cdef class CountNeighbors(SPHFunctionParticle):
         SPHFunctionParticle.__init__(self, source, dest, setup_arrays = True)
         self.id = 'nbrs'
 
+    def set_src_dst_reads(self):
+        pass
+
     cdef void eval_single(self, size_t dest_pid,
                           KernelBase kernel, double *result):
         result[0] += self.nbr_locator.get_nearest_particles(dest_pid).length
diff --git a/source/pysph/sph/funcs/boundary_funcs.pyx b/source/pysph/sph/funcs/boundary_funcs.pyx
index 7c7240c..e52f089 100644
--- a/source/pysph/sph/funcs/boundary_funcs.pyx
+++ b/source/pysph/sph/funcs/boundary_funcs.pyx
@@ -40,8 +40,17 @@ cdef class MonaghanBoundaryForce(SPHFunctionParticle):
         self.s_ny = self.source.get_carray("ny")
         self.s_nz = self.source.get_carray("nz")
 
+        
+
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','m','rho'] )
         self.src_reads.extend( ['tx','ty','tz','nx','ny','nz','cs'] )
 
+        self.dst_reads.extend( ['x','y','z','h','m','cs','tag'] )
+
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
                        KernelBase kernel, double *nr):
         """ Perform the boundary force computation """
@@ -61,10 +70,12 @@ cdef class MonaghanBoundaryForce(SPHFunctionParticle):
         self._dst.y = self.d_y.data[dest_pid]
         self._dst.z = self.d_z.data[dest_pid]
             
-        cdef cPoint norm = cPoint(self.s_nx.data[source_pid], self.s_ny.data[source_pid],
+        cdef cPoint norm = cPoint(self.s_nx.data[source_pid],
+                                  self.s_ny.data[source_pid],
                                   self.s_nz.data[source_pid])
         
-        cdef cPoint tang = cPoint(self.s_tx.data[source_pid], self.s_ty.data[source_pid],
+        cdef cPoint tang = cPoint(self.s_tx.data[source_pid],
+                                  self.s_ty.data[source_pid],
                                   self.s_tz.data[source_pid])
 
         cs = self.d_cs.data[dest_pid]
@@ -125,6 +136,9 @@ cdef class BeckerBoundaryForce(SPHFunctionParticle):
         self.cl_kernel_src_file = "boundary_funcs.cl"
         self.cl_kernel_function_name = "BeckerBoundaryForce"
 
+    def set_src_dst_reads(self):
+        pass
+
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
                        KernelBase kernel, double *nr):
         """
@@ -213,6 +227,9 @@ cdef class LennardJonesForce(SPHFunctionParticle):
         self.cl_kernel_src_file = "boundary_funcs.cl"
         self.cl_kernel_function_name = "LennardJonesForce"
 
+    def set_src_dst_reads(self):
+        pass
+
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
                        KernelBase kernel, double *nr):
         """
diff --git a/source/pysph/sph/funcs/density_funcs.pyx b/source/pysph/sph/funcs/density_funcs.pyx
index ce33c51..fe657a5 100644
--- a/source/pysph/sph/funcs/density_funcs.pyx
+++ b/source/pysph/sph/funcs/density_funcs.pyx
@@ -27,6 +27,13 @@ cdef class SPHRho(CSPHFunctionParticle):
         self.cl_kernel_src_file = "density_funcs.clt"
         self.cl_kernel_function_name = "SPHRho"
 
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','m'] )
+        self.dst_reads.extend( ['x','y','z','h','tag'] )
+
     cdef void eval_nbr_csph(self, size_t source_pid, size_t dest_pid,
                             KernelBase kernel, double *nr, double *dnr):
         """ Compute the contribution from source_pid on dest_pid. """
@@ -89,12 +96,19 @@ cdef class SPHDensityRate(SPHFunctionParticle):
         self.id = 'densityrate'
         self.tag = "density"
 
-        self.src_reads.extend( ['u','v','w'] )
-        self.dst_reads.extend( ['u','v','w'] )
-
         self.cl_kernel_src_file = "density_funcs.cl"
         self.cl_kernel_function_name = "SPHDensityRate"
 
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','m'] )
+        self.dst_reads.extend( ['x','y','z','h','tag'] )
+
+        self.src_reads.extend( ['u','v','w'] )
+        self.dst_reads.extend( ['u','v','w'] )
+
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
                        KernelBase kernel, double *nr):
         """ Compute the contribution of particle at source_pid on particle at
diff --git a/source/pysph/sph/funcs/energy_funcs.pyx b/source/pysph/sph/funcs/energy_funcs.pyx
index 48f532f..709ec7f 100644
--- a/source/pysph/sph/funcs/energy_funcs.pyx
+++ b/source/pysph/sph/funcs/energy_funcs.pyx
@@ -19,11 +19,18 @@ cdef class EnergyEquationNoVisc(SPHFunctionParticle):
         self.id = 'energyeqn'
         self.tag = "energy"
 
-        self.src_reads.extend( ['u','v','w','p'] )
-        self.dst_reads.extend( ['u','v','w','p'] )
-
         self.cl_kernel_src_file = "energy_funcs.cl"
         self.cl_kernel_function_name = "EnergyEquationNoVisc"
+
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','m','rho'] )
+        self.dst_reads.extend( ['x','y','z','h','rho','tag'] )
+
+        self.src_reads.extend( ['u','v','w','p'] )
+        self.dst_reads.extend( ['u','v','w','p'] )
     
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
                        KernelBase kernel, double *nr):
@@ -107,11 +114,18 @@ cdef class EnergyEquationAVisc(SPHFunctionParticle):
         self.id = 'energyavisc'
         self.tag = "energy"
 
-        self.src_reads.extend( ['u','v','w','p','cs'] )
-        self.dst_reads.extend( ['u','v','w','p','cs'] )
-
         self.cl_kernel_src_file = "energy_funcs.cl"
         self.cl_kernel_function_name = "EnergyEquationAVisc"
+
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','m','rho'] )
+        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'] )
     
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
                        KernelBase kernel, double *nr):
@@ -222,15 +236,21 @@ cdef class EnergyEquation(SPHFunctionParticle):
         self.id = 'energyequation'
         self.tag = "energy"
 
-        self.src_reads.extend( ['u','v','w','p','cs'] )
-        self.dst_reads.extend( ['u','v','w','p','cs'] )
-
         self.cl_kernel_src_file = "energy_funcs.cl"
         self.cl_kernel_function_name = "EnergyEquation"
+
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','h','m','rho'] )
+        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'] )        
         
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
-                       KernelBase kernel, double *nr):
-    
+                       KernelBase kernel, double *nr):    
     
         cdef cPoint vab, rab
         cdef double Pa, Pb, rhoa, rhob, rhoab, mb
diff --git a/source/pysph/sph/funcs/eos_funcs.pyx b/source/pysph/sph/funcs/eos_funcs.pyx
index dc59342..0b9220a 100644
--- a/source/pysph/sph/funcs/eos_funcs.pyx
+++ b/source/pysph/sph/funcs/eos_funcs.pyx
@@ -20,14 +20,17 @@ cdef class IdealGasEquation(SPHFunction):
         self.id = 'idealgas'
         self.tag = "state"
 
-        self.dst_reads.extend( ['e','rho', 'p'] )
-
         self.cl_kernel_src_file = "eos_funcs.cl"
 
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.dst_reads.extend( ['e','rho'] )
+
     cdef void eval_single(self, size_t dest_pid, KernelBase kernel,
                           double* result):
         
-        cdef double Pa = self.d_p.data[dest_pid]
         cdef double ea = self.d_e.data[dest_pid]
         cdef double rhoa = self.d_rho.data[dest_pid]
         cdef double gamma = self.gamma
@@ -75,10 +78,14 @@ cdef class TaitEquation(SPHFunction):
         self.id = 'tait'
         self.tag = "state"
 
-        self.dst_reads.append('rho')
-
         self.cl_kernel_src_file = "eos_funcs.cl"
 
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.dst_reads.extend( ['rho'] )
+
     cdef void eval_single(self, size_t dest_pid, KernelBase kernel,
                           double* result):
 
diff --git a/source/pysph/sph/funcs/external_force.clt b/source/pysph/sph/funcs/external_force.clt
index a43423e..925f9a6 100644
--- a/source/pysph/sph/funcs/external_force.clt
+++ b/source/pysph/sph/funcs/external_force.clt
@@ -28,7 +28,7 @@ __kernel void NBodyForce(%(kernel_args)s)
 
     // The term `dest_id` will be suitably defined at this point.
 
-    REAL4 pa = (REAL4)( d_x[dest_id],d_y[dest_id], d_z[dest_id], d_h[dest_id] );
+    REAL4 pa = (REAL4)( d_x[dest_id],d_y[dest_id], d_z[dest_id], 0.0F );
     REAL4 rba;
     REAL invr, force_mag;
 
@@ -37,7 +37,7 @@ __kernel void NBodyForce(%(kernel_args)s)
         // SPH innermost loop code goes here.  The index `src_id` will
         // be available and looped over, this index.
 
-        REAL4 pb = (REAL4)( s_x[src_id],s_y[src_id],s_z[src_id],s_h[src_id] );
+        REAL4 pb = (REAL4)( s_x[src_id],s_y[src_id],s_z[src_id], 0.0F );
         rba = pb - pa;
       
         invr = 1.0F/( length(rba) + eps );
diff --git a/source/pysph/sph/funcs/external_force.pyx b/source/pysph/sph/funcs/external_force.pyx
index 4f0937c..142614d 100644
--- a/source/pysph/sph/funcs/external_force.pyx
+++ b/source/pysph/sph/funcs/external_force.pyx
@@ -33,6 +33,9 @@ cdef class GravityForce(SPHFunction):
         self.cl_kernel_src_file = "external_force.clt"
         self.cl_kernel_function_name = "GravityForce"
 
+    def set_src_dst_reads(self):
+        pass
+
     cdef void eval_single(self, size_t dest_pid,
                           KernelBase kernel, double *result):
         """ Perform the gravity force computation """
@@ -79,6 +82,9 @@ cdef class VectorForce(SPHFunction):
         self.cl_kernel_src_file = "external_force.cl"
         self.cl_kernel_function_name = "VectorForce"
 
+    def set_src_dst_reads(self):
+        pass
+
     cdef void eval_single(self, size_t dest_pid,
                           KernelBase kernel, double *result):
         """ Perform the force computation """
@@ -104,6 +110,9 @@ cdef class MoveCircleX(SPHFunction):
         self.id = 'circlex'
         self.tag = "position"
 
+    def set_src_dst_reads(self):
+        pass        
+
     cdef void eval_single(self, size_t dest_pid,
                           KernelBase kernel, double *result):
         cdef cPoint p = cPoint(self.d_x.data[dest_pid],
@@ -136,6 +145,9 @@ cdef class MoveCircleY(SPHFunction):
         self.id = 'circley'
         self.tag = "position"
 
+    def set_src_dst_reads(self):
+        pass        
+
     cdef void eval_single(self, size_t dest_pid,
                           KernelBase kernel, double *result):
         cdef cPoint p = cPoint(self.d_x.data[dest_pid],
@@ -177,6 +189,13 @@ cdef class NBodyForce(SPHFunctionParticle):
         self.cl_kernel_src_file = "external_force.clt"
         self.cl_kernel_function_name = "NBodyForce"
 
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.src_reads.extend( ['x','y','z','m'] )
+        self.dst_reads.extend( ['x','y','z','m','tag'] )
+
     cdef void eval_single(self, size_t dest_pid, KernelBase kernel,
                           double * result):
         """ Neighbors for the NBody example are by default all neighbors
diff --git a/source/pysph/sph/funcs/position_funcs.pyx b/source/pysph/sph/funcs/position_funcs.pyx
index 0d15f67..68b2099 100644
--- a/source/pysph/sph/funcs/position_funcs.pyx
+++ b/source/pysph/sph/funcs/position_funcs.pyx
@@ -19,10 +19,14 @@ cdef class PositionStepping(SPHFunction):
         self.id = 'positionstepper'
         self.tag = "position"
 
-        self.dst_reads.extend( ['u','v','w'] )
-
         self.cl_kernel_src_file = "position_funcs.clt"
         self.cl_kernel_function_name = "PositionStepping"
+
+    def set_src_dst_reads(self):
+        self.src_reads = []
+        self.dst_reads = []
+
+        self.dst_reads.extend( ['u','v','w'] )        
     
     cpdef eval(self, KernelBase kernel, DoubleArray output1,
                DoubleArray output2, DoubleArray output3):
diff --git a/source/pysph/sph/funcs/pressure_funcs.pyx b/source/pysph/sph/funcs/pressure_funcs.pyx
index 25657ea..0fd64f9 100644
--- a/source/pysph/sph/funcs/pressure_funcs.pyx
+++ b/source/pysph/sph/funcs/pressure_funcs.pyx
@@ -24,12 +24,13 @@ cdef class SPHPressureGradient(SPHFunctionParticle):
         self.id = 'pgrad'
         self.tag = "velocity"
 
-        self.dst_reads.extend( ['p', 'rho'] )
-        self.src_reads.append( 'p' )
-
         self.cl_kernel_src_file = "pressure_funcs.cl"
         self.cl_kernel_function_name = "SPHPressureGradient"
 
+    def set_src_dst_reads(self):
+        self.src_reads = ['x','y','z','h','m','rho','p']
+        self.dst_reads = ['x','y','z','h','rho','p','tag']
+
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid,
                    KernelBase kernel, double *nr):
         cdef double mb = self.s_m.data[source_pid]
@@ -113,11 +114,13 @@ cdef class MomentumEquation(SPHFunctionParticle):
         self.id = 'momentumequation'
         self.tag = "velocity"
 
-        self.src_reads.extend( ['u','v','w','p','cs'] )
-        self.src_reads.extend( ['u','v','w','p','cs','rho'] )
-
         self.cl_kernel_src_file = "pressure_funcs.cl"
         self.cl_kernel_function_name = "MomentumEquation"
+
+    def set_src_dst_reads(self):
+        self.src_reads = ['x','y','z','h','m','rho','p','u','v','w','cs']
+        self.dst_reads = ['x','y','z','h','rho','p',
+                          'u','v','w','cs','tag']
         
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid,
                        KernelBase kernel, double *nr):
diff --git a/source/pysph/sph/funcs/viscosity_funcs.pyx b/source/pysph/sph/funcs/viscosity_funcs.pyx
index 85e13de..dfb8e03 100644
--- a/source/pysph/sph/funcs/viscosity_funcs.pyx
+++ b/source/pysph/sph/funcs/viscosity_funcs.pyx
@@ -30,12 +30,14 @@ cdef class MonaghanArtificialVsicosity(SPHFunctionParticle):
         self.id = 'momavisc'
         self.tag = "velocity"
 
-        self.src_reads.extend( ['u','v','w','cs'] )
-        self.dst_reads.extend( ['u','v','w','cs','rho'] )
-
         self.cl_kernel_src_file = "viscosity_funcs.cl"
         self.cl_kernel_function_name = "MonaghanArtificialVsicosity"
 
+    def set_src_dst_reads(self):
+        self.src_reads = ['x','y','z','h','m','rho','u','v','w','cs']
+        self.dst_reads = ['x','y','z','h','p',
+                          'u','v','w','cs','rho','tag']
+
     cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, 
                        KernelBase kernel, double *nr):
         cdef cPoint va, vb, vab, rab
@@ -150,6 +152,10 @@ cdef class MorrisViscosity(SPHFunctionParticle):
         self.cl_kernel_src_file = "viscosity_funcs.cl"
         self.cl_kernel_function_name = "MorrisViscosity"
 
+    def set_src_dst_reads(self):
+        self.src_reads = ['x','y','z','h','m','rho','u','v','w',self.mu]
+        self.dst_reads = ['x','y','z','h','rho','u','v','w','tag',self.mu]
+
     cpdef setup_arrays(self):
         """
         """
diff --git a/source/pysph/sph/funcs/xsph_funcs.pyx b/source/pysph/sph/funcs/xsph_funcs.pyx
index f85defb..387a787 100644
--- a/source/pysph/sph/funcs/xsph_funcs.pyx
+++ b/source/pysph/sph/funcs/xsph_funcs.pyx
@@ -20,12 +20,13 @@ cdef class XSPHCorrection(CSPHFunctionParticle):
         self.id = 'xsph'
         self.tag = "position"
 
-        self.src_reads.extend( ['u','v','w'] )
-        self.dst_reads.extend( ['u','v','w','rho'] )
-
         self.cl_kernel_src_file = "xsph_funcs.cl"
         self.cl_kernel_function_name = "XSPHCorrection"
 
+    def set_src_dst_reads(self):
+        self.src_reads = ['x','y','z','h','m','rho','u','v','w']
+        self.dst_reads = ['x','y','z','h','rho','u','v','w']
+
     cdef void eval_nbr_csph(self, size_t source_pid, size_t dest_pid,
                             KernelBase kernel, double *nr, double *dnr):
         """
diff --git a/source/pysph/sph/sph_func.pyx b/source/pysph/sph/sph_func.pyx
index 87a376b..3c1290a 100755
--- a/source/pysph/sph/sph_func.pyx
+++ b/source/pysph/sph/sph_func.pyx
@@ -125,8 +125,8 @@ cdef class SPHFunction:
         
         self.num_outputs = 3
 
-        self.src_reads = ['x','y','z','h','m', 'rho']
-        self.dst_reads = ['x','y','z','h','tag']
+        self.src_reads = []
+        self.dst_reads = []
 
         self.kernel = None
 
@@ -140,6 +140,9 @@ cdef class SPHFunction:
 
         self.global_sizes = (self.dest.get_number_of_particles(), 1, 1)
         self.local_sizes = (1,1,1)
+
+        # setup the source and destination reads
+        self.set_src_dst_reads()
         
         if setup_arrays:
             self.setup_arrays()
@@ -300,7 +303,21 @@ cdef class SPHFunction:
 
     def get_cl_workgroup_code(self):
         return """unsigned int work_dim = get_work_dim();
-    unsigned int dest_id = get_gid(work_dim); """        
+    unsigned int dest_id = get_gid(work_dim); """
+
+    def set_src_dst_reads(self):
+        """ Populate the read requirements for the Function
+
+        The read requirements specify which particle properties will
+        be required by this function. Properties read from the source
+        particle array are appended to the list `src_reads` and those
+        from the destinatio particle array are appended to `dst_reads`
+
+        These read requirements are used to construct the OpenCL
+        kernel arguments at program creation time.
+
+        """
+        raise NotImplementedError("SPHFunction set_src_dst_reads called!")
 
 ################################################################################
 # `SPHFunctionParticle` class.