diff --git a/Makefile b/Makefile
index 64131f4..332e8b5 100644
--- a/Makefile
+++ b/Makefile
@@ -73,10 +73,10 @@ epydoc :
 doc :
 	cd docs; make html
 
-develop : $(DIRS)
+develop :
 	python setup.py develop
 
-install : $(DIRS)
+install :
 	python setup.py install
 
 clang :
diff --git a/source/pysph/base/particle_array.pyx b/source/pysph/base/particle_array.pyx
index efe76b5..7faa02a 100644
--- a/source/pysph/base/particle_array.pyx
+++ b/source/pysph/base/particle_array.pyx
@@ -202,13 +202,12 @@ cdef class ParticleArray:
         default_values = {}
 
         for prop, arr in self.properties.iteritems():
-            if not prop.startswith('_'):
-                pinfo = {}
-                pinfo['name'] = prop            
-                pinfo['type'] = arr.get_c_type()
-                pinfo['data'] = arr.get_npy_array()
-                pinfo['default'] = self.default_values[prop]
-                props[prop] = pinfo
+            pinfo = {}
+            pinfo['name'] = prop            
+            pinfo['type'] = arr.get_c_type()
+            pinfo['data'] = arr.get_npy_array()
+            pinfo['default'] = self.default_values[prop]
+            props[prop] = pinfo
         
         d['properties'] = props
 
diff --git a/source/pysph/base/particles.py b/source/pysph/base/particles.py
index 9b49307..fe99b28 100644
--- a/source/pysph/base/particles.py
+++ b/source/pysph/base/particles.py
@@ -312,8 +312,7 @@ def get_particle_array(cl_precision="double", **props):
 
     default_props = {'x':0.0, 'y':0.0, 'z':0.0, 'u':0.0, 'v':0.0 ,
                      'w':0.0, 'm':1.0, 'h':1.0, 'p':0.0,'e':0.0,
-                     'rho':1.0, 'cs':0.0, '_tmpx':0.0,
-                     '_tmpy':0.0, '_tmpz':0.0}
+                     'rho':1.0, 'cs':0.0}
     
     #Add the properties requested
     
diff --git a/source/pysph/parallel/tests/test_particle_array_pickling.py b/source/pysph/parallel/tests/particle_array_pickling.py
similarity index 81%
rename from source/pysph/parallel/tests/test_particle_array_pickling.py
rename to source/pysph/parallel/tests/particle_array_pickling.py
index b6110ad..23d82cc 100644
--- a/source/pysph/parallel/tests/test_particle_array_pickling.py
+++ b/source/pysph/parallel/tests/particle_array_pickling.py
@@ -32,9 +32,6 @@
     # make sure the sending data defines the underscored variables
     assert pa.properties.has_key('_rho')
     assert pa.properties.has_key('_p')
-    assert pa.properties.has_key('_tmpx')
-    assert pa.properties.has_key('_tmpy')
-    assert pa.properties.has_key('_tmpz')
 
     comm.send(obj=pa, dest=1)
 
@@ -45,9 +42,6 @@
     # make sure that the underscored variables are not there
     assert not pa2.properties.has_key('_rho')
     assert not pa2.properties.has_key('_p')
-    assert not pa2.properties.has_key('_tmpx')
-    assert not pa2.properties.has_key('_tmpy')
-    assert not pa2.properties.has_key('_tmpz')
 
     # now append the received array to the local array
     pa.append_parray(pa2)
diff --git a/source/pysph/parallel/tests/test_parallel_runs.py b/source/pysph/parallel/tests/test_parallel_runs.py
index e3a7d37..e6f6488 100644
--- a/source/pysph/parallel/tests/test_parallel_runs.py
+++ b/source/pysph/parallel/tests/test_parallel_runs.py
@@ -84,6 +84,10 @@ def test_parallel_cell_check3(self):
     def test_share_data(self):
         for i in range(1,6):
             run_mpi_script('share_data.py', i)
+
+    def test_particle_array_pickling(self):
+        for i in range(2,3):
+            run_mpi_script('particle_array_pickling.py', i)
         
 
 #for filename in tests:
diff --git a/source/pysph/solver/fluid_solver.py b/source/pysph/solver/fluid_solver.py
index 656f070..b534233 100644
--- a/source/pysph/solver/fluid_solver.py
+++ b/source/pysph/solver/fluid_solver.py
@@ -75,7 +75,7 @@ def setup_solver(self):
         self.add_operation(SPHIntegration(
             sph.MomentumEquation.withargs(alpha=0.01, beta=0.0, hks=False),
             from_types=[Fluids], on_types=[Fluids], 
-            updates=['u','v'], 
+            updates=['u','v','w'], 
             id='mom')
                            )
 
diff --git a/source/pysph/solver/integrator.py b/source/pysph/solver/integrator.py
index 49cd5eb..909b08f 100644
--- a/source/pysph/solver/integrator.py
+++ b/source/pysph/solver/integrator.py
@@ -222,8 +222,6 @@ def __init__(self, particles=None, calcs=[], pcalcs = []):
 
         self.rupdate_list = []
 
-        self.step_props = '_tmpx', '_tmpy', '_tmpz'
-
     def set_rupdate_list(self):
         for i in range(len(self.particles.arrays)):
             self.rupdate_list.append([])
diff --git a/source/pysph/solver/solver.py b/source/pysph/solver/solver.py
index 27e7a7b..d1f5081 100644
--- a/source/pysph/solver/solver.py
+++ b/source/pysph/solver/solver.py
@@ -141,8 +141,8 @@ def add_operation_step(self, types, xsph=False, eps=0.5):
         id = 'step'
         
         self.add_operation(SPHIntegration(
-            sph.PositionStepping, on_types=types, updates=updates, id=id,
-            kernel=kernel)
+                sph.PositionStepping.withargs(dim=self.dim), on_types=types,
+                updates=updates, id=id, kernel=kernel)
                            )
 
     def add_operation_xsph(self, eps, hks=False):
@@ -171,13 +171,12 @@ def add_operation_xsph(self, eps, hks=False):
         assert self.operation_dict.has_key('step'), err
 
         types = self.operation_dict['step'].on_types
-        from_types = self.operation_dict['step'].from_types
         updates = self.operation_dict['step'].updates
 
                            
         self.add_operation(SPHIntegration(
 
-            sph.XSPHCorrection.withargs(eps=eps, hks=hks), from_types=from_types,
+            sph.XSPHCorrection.withargs(eps=eps, hks=hks, dim=self.dim), from_types=types,
             on_types=types, updates=updates, id=id, kernel=self.default_kernel)
 
                            )        
diff --git a/source/pysph/sph/funcs/adke_funcs.pyx b/source/pysph/sph/funcs/adke_funcs.pyx
index c6d4840..d6b5c3c 100644
--- a/source/pysph/sph/funcs/adke_funcs.pyx
+++ b/source/pysph/sph/funcs/adke_funcs.pyx
@@ -27,6 +27,7 @@ cdef class ADKEPilotRho(CSPHFunctionParticle):
 
         self.cl_kernel_src_file = "adke_funcs.cl"
         self.cl_kernel_function_name = "ADKEPilotRho"
+        self.num_outputs = 1
 
     def set_src_dst_reads(self):
 
@@ -94,6 +95,7 @@ cdef class ADKESmoothingUpdate(SPHFunction):
 
         self.cl_kernel_src_file = "adke_funcs.cl"
         self.cl_kernel_function_name = "ADKESmoothingUpdate"
+        self.num_outputs = 1
 
     cpdef setup_arrays(self):
         """
@@ -237,6 +239,7 @@ cdef class ADKEConductionCoeffUpdate(SPHFunction):
         self.cl_kernel_src_file = "adke_funcs.cl"
 
         self.cl_kernel_function_name = "ADKEConductionCoeffUpdate"
+        self.num_outputs = 1
 
     def set_src_dst_reads(self):
         self.src_reads = []
diff --git a/source/pysph/sph/funcs/arithmetic_funcs.pxd b/source/pysph/sph/funcs/arithmetic_funcs.pxd
index e8017bc..d075252 100644
--- a/source/pysph/sph/funcs/arithmetic_funcs.pxd
+++ b/source/pysph/sph/funcs/arithmetic_funcs.pxd
@@ -9,6 +9,7 @@ cdef class PropertyGet(SPHFunction):
 cdef class PropertyAdd(SPHFunction):
     cdef list prop_names, d_props
     cdef int num_props
+    cdef public double constant
 
 cdef class PropertyNeg(SPHFunction):
     cdef list prop_names, d_props
@@ -16,6 +17,7 @@ cdef class PropertyNeg(SPHFunction):
 cdef class PropertyMul(SPHFunction):
     cdef list prop_names, d_props
     cdef int num_props
+    cdef public double constant
 
 cdef class PropertyInv(SPHFunction):
     cdef list prop_names, d_props
diff --git a/source/pysph/sph/funcs/arithmetic_funcs.pyx b/source/pysph/sph/funcs/arithmetic_funcs.pyx
index 12f95ea..8cf7a4f 100644
--- a/source/pysph/sph/funcs/arithmetic_funcs.pyx
+++ b/source/pysph/sph/funcs/arithmetic_funcs.pyx
@@ -37,7 +37,7 @@ cdef class PropertyGet(SPHFunction):
 cdef class PropertyAdd(SPHFunction):
     """ function to add property arrays (any number of arrays) """
     def __init__(self, ParticleArray source, ParticleArray dest=None,
-                 prop_names=['rho'], **kwargs):
+                 prop_names=['rho'], constant=0.0, **kwargs):
         """ Constructor for SPH
 
         Parameters:
@@ -45,8 +45,10 @@ cdef class PropertyAdd(SPHFunction):
         source -- The source particle array
         dest -- The destination particle array
         prop_names -- The properties to get the sum of
+        constant -- a constant value to add to the result
         """
         self.prop_names = prop_names
+        self.constant = constant
         SPHFunction.__init__(self, source, dest, setup_arrays=True)
         self.num_outputs = 1
         self.id = 'property_add'
@@ -55,7 +57,7 @@ cdef class PropertyAdd(SPHFunction):
         #Setup the basic properties like m, x rho etc.
         SPHFunction.setup_arrays(self)
         self.d_props = [self.dest.get_carray(i) for i in self.prop_names]
-        self.num_props = len(self.prop_names)
+        self.num_props = len(self.d_props)
     
     cpdef eval(self, KernelBase kernel, DoubleArray output1,
                DoubleArray output2, DoubleArray output3):
@@ -67,7 +69,7 @@ cdef class PropertyAdd(SPHFunction):
         for n in range(1, self.num_props):
             arr = self.d_props[n]
             for i in range(np):
-                output1.data[i] += arr.data[i]
+                output1.data[i] += arr.data[i] + self.constant
 
 cdef class PropertyNeg(SPHFunction):
     """ function to return the negative of upto 3 particle arrays """
@@ -105,7 +107,7 @@ cdef class PropertyNeg(SPHFunction):
 cdef class PropertyMul(SPHFunction):
     """ function to get product of property arrays (any number of arrays) """
     def __init__(self, ParticleArray source, ParticleArray dest=None,
-                 prop_names=['rho'], **kwargs):
+                 prop_names=['rho'], constant=0.0, **kwargs):
         """ Constructor for SPH
 
         Parameters:
@@ -113,8 +115,10 @@ cdef class PropertyMul(SPHFunction):
         source -- The source particle array.
         dest -- The destination particle array.
         prop_names -- The properties to get product of
+        constant -- a constant value to multiply to the result
         """
         self.prop_names = prop_names
+        self.constant = constant
         SPHFunction.__init__(self, source, dest, setup_arrays=True)
         self.num_outputs = 1
         self.id = 'property_mul'
@@ -123,7 +127,7 @@ cdef class PropertyMul(SPHFunction):
         #Setup the basic properties like m, x rho etc.
         SPHFunction.setup_arrays(self)
         self.d_props = [self.dest.get_carray(i) for i in self.prop_names]
-        self.num_props = len(self.prop_names)
+        self.num_props = len(self.d_props)
     
     cpdef eval(self, KernelBase kernel, DoubleArray output1,
                DoubleArray output2, DoubleArray output3):
@@ -135,7 +139,7 @@ cdef class PropertyMul(SPHFunction):
         for n in range(1, self.num_props):
             arr = self.d_props[n]
             for i in range(np):
-                output1.data[i] *= arr.data[i]
+                output1.data[i] *= arr.data[i] * self.constant
 
 cdef class PropertyInv(SPHFunction):
     """ function to return the inverse of upto 3 particle arrays """
diff --git a/source/pysph/sph/funcs/basic_funcs.pyx b/source/pysph/sph/funcs/basic_funcs.pyx
index 79cc29b..206bede 100644
--- a/source/pysph/sph/funcs/basic_funcs.pyx
+++ b/source/pysph/sph/funcs/basic_funcs.pyx
@@ -32,6 +32,7 @@ cdef class SPH(CSPHFunctionParticle):
 
         self.cl_kernel_src_file = "basic_funcs.clt"
         self.cl_kernel_function_name = "SPH"
+        self.num_outputs = 1
 
     def set_src_dst_reads(self):
         self.src_reads = []
@@ -337,6 +338,7 @@ cdef class SPHLaplacian(SPHFunctionParticle):
 
         self.cl_kernel_src_file = "basic_funcs.cl"
         self.cl_kernel_function_name = "SPHLaplacian"
+        self.num_outputs = 1
 
     def set_src_dst_reads(self):
         self.src_reads = []
diff --git a/source/pysph/sph/funcs/density_funcs.pyx b/source/pysph/sph/funcs/density_funcs.pyx
index 45a109c..0864fa4 100644
--- a/source/pysph/sph/funcs/density_funcs.pyx
+++ b/source/pysph/sph/funcs/density_funcs.pyx
@@ -26,6 +26,7 @@ cdef class SPHRho(CSPHFunctionParticle):
 
         self.cl_kernel_src_file = "density_funcs.clt"
         self.cl_kernel_function_name = "SPHRho"
+        self.num_outputs = 1
 
     def set_src_dst_reads(self):
         self.src_reads = []
@@ -98,6 +99,7 @@ cdef class SPHDensityRate(SPHFunctionParticle):
 
         self.cl_kernel_src_file = "density_funcs.cl"
         self.cl_kernel_function_name = "SPHDensityRate"
+        self.num_outputs = 1
 
     def set_src_dst_reads(self):
         self.src_reads = []
diff --git a/source/pysph/sph/funcs/energy_funcs.pyx b/source/pysph/sph/funcs/energy_funcs.pyx
index 845f912..d5428ad 100644
--- a/source/pysph/sph/funcs/energy_funcs.pyx
+++ b/source/pysph/sph/funcs/energy_funcs.pyx
@@ -23,6 +23,7 @@ cdef class EnergyEquationNoVisc(SPHFunctionParticle):
 
         self.cl_kernel_src_file = "energy_funcs.cl"
         self.cl_kernel_function_name = "EnergyEquationNoVisc"
+        self.num_outputs = 1
 
     def set_src_dst_reads(self):
         self.src_reads = []
@@ -128,6 +129,7 @@ cdef class EnergyEquationAVisc(SPHFunctionParticle):
 
         self.cl_kernel_src_file = "energy_funcs.cl"
         self.cl_kernel_function_name = "EnergyEquationAVisc"
+        self.num_outputs = 1
 
     def set_src_dst_reads(self):
         self.src_reads = []
@@ -251,6 +253,7 @@ cdef class EnergyEquation(SPHFunctionParticle):
 
         self.cl_kernel_src_file = "energy_funcs.cl"
         self.cl_kernel_function_name = "EnergyEquationWithVisc"
+        self.num_outputs = 1
 
     def set_src_dst_reads(self):
         self.src_reads = []
@@ -400,6 +403,7 @@ cdef class ArtificialHeat(SPHFunctionParticle):
 
         self.cl_kernel_src_file = "energy_funcs.cl"
         self.cl_kernel_function_name = "ArtificialHeat"
+        self.num_outputs = 1
 
     def set_src_dst_reads(self):
         self.src_reads = ['x','y','z','h','m','rho','e,','u','v','w','cs','q']
diff --git a/source/pysph/sph/funcs/eos_funcs.pyx b/source/pysph/sph/funcs/eos_funcs.pyx
index 0b3c653..8244f53 100644
--- a/source/pysph/sph/funcs/eos_funcs.pyx
+++ b/source/pysph/sph/funcs/eos_funcs.pyx
@@ -24,6 +24,7 @@ cdef class IdealGasEquation(SPHFunction):
 
         self.cl_kernel_src_file = "eos_funcs.cl"
         self.cl_kernel_function_name = "IdealGasEquation"
+        self.num_outputs = 2
 
     def set_src_dst_reads(self):
         self.src_reads = []
@@ -93,6 +94,7 @@ cdef class TaitEquation(SPHFunction):
         self.tag = "state"
 
         self.cl_kernel_src_file = "eos_funcs.cl"
+        self.num_outputs = 2
 
     def set_src_dst_reads(self):
         self.src_reads = []
diff --git a/source/pysph/sph/funcs/external_force.pyx b/source/pysph/sph/funcs/external_force.pyx
index 280435a..b58e3b2 100644
--- a/source/pysph/sph/funcs/external_force.pyx
+++ b/source/pysph/sph/funcs/external_force.pyx
@@ -109,6 +109,7 @@ cdef class MoveCircleX(SPHFunction):
 
         self.id = 'circlex'
         self.tag = "position"
+        self.num_outputs = 2
 
     def set_src_dst_reads(self):
         pass        
@@ -144,6 +145,7 @@ cdef class MoveCircleY(SPHFunction):
 
         self.id = 'circley'
         self.tag = "position"
+        self.num_outputs = 2
 
     def set_src_dst_reads(self):
         pass        
diff --git a/source/pysph/sph/funcs/position_funcs.pyx b/source/pysph/sph/funcs/position_funcs.pyx
index ea77910..5bdab42 100644
--- a/source/pysph/sph/funcs/position_funcs.pyx
+++ b/source/pysph/sph/funcs/position_funcs.pyx
@@ -12,7 +12,7 @@ cdef class PositionStepping(SPHFunction):
 
     #Defined in the .pxd file
     def __init__(self, ParticleArray source, ParticleArray dest, 
-                 bint setup_arrays=True):
+                 bint setup_arrays=True, int dim=3):
 
         SPHFunction.__init__(self, source, dest, setup_arrays)
         
@@ -21,12 +21,13 @@ cdef class PositionStepping(SPHFunction):
 
         self.cl_kernel_src_file = "position_funcs.clt"
         self.cl_kernel_function_name = "PositionStepping"
+        self.num_outputs = dim
 
     def set_src_dst_reads(self):
         self.src_reads = []
         self.dst_reads = []
 
-        self.dst_reads.extend( ['u','v','w'] )        
+        self.dst_reads.extend( ['u','v','w'][:self.num_outputs] )
     
     cpdef eval(self, KernelBase kernel, DoubleArray output1,
                DoubleArray output2, DoubleArray output3):
@@ -35,14 +36,29 @@ cdef class PositionStepping(SPHFunction):
 
         self.setup_iter_data()
         cdef size_t np = self.dest.get_number_of_particles()
-        
-        for i in range(np):
-            if tag_arr.data[i] == LocalReal:
-                output1[i] = self.d_u.data[i]
-                output2[i] = self.d_v.data[i]
-                output3[i] = self.d_w.data[i]
-            else:
-                output1[i] = output2[i] = output3[i] = 0
+
+        if self.num_outputs == 3:
+            for i in range(np):
+                if tag_arr.data[i] == LocalReal:
+                    output1.data[i] += self.d_u.data[i]
+                    output2.data[i] += self.d_v.data[i]
+                    output3.data[i] += self.d_w.data[i]
+                else:
+                    output1.data[i] = output2.data[i] = output3.data[i] = 0
+        elif self.num_outputs == 2:
+            for i in range(np):
+                if tag_arr.data[i] == LocalReal:
+                    output1.data[i] += self.d_u.data[i]
+                    output2.data[i] += self.d_v.data[i]
+                else:
+                    output1.data[i] = output2.data[i] = 0
+        elif self.num_outputs == 1:
+            for i in range(np):
+                if tag_arr.data[i] == LocalReal:
+                    output1.data[i] += self.d_u.data[i]
+                else:
+                    output1.data[i] = 0
+
 
     def _set_extra_cl_args(self):
         pass
diff --git a/source/pysph/sph/funcs/pressure_funcs.pyx b/source/pysph/sph/funcs/pressure_funcs.pyx
index c5e2274..bb11a39 100644
--- a/source/pysph/sph/funcs/pressure_funcs.pyx
+++ b/source/pysph/sph/funcs/pressure_funcs.pyx
@@ -18,7 +18,7 @@ cdef class SPHPressureGradient(SPHFunctionParticle):
     """
 
     def __init__(self, ParticleArray source, ParticleArray dest,
-                 bint setup_arrays=True, **kwargs):
+                 bint setup_arrays=True, int dim=3, **kwargs):
         
         SPHFunctionParticle.__init__(self, source, dest, setup_arrays,
                                      **kwargs)
@@ -28,6 +28,7 @@ cdef class SPHPressureGradient(SPHFunctionParticle):
 
         self.cl_kernel_src_file = "pressure_funcs.cl"
         self.cl_kernel_function_name = "SPHPressureGradient"
+        self.num_outputs = dim
 
     def set_src_dst_reads(self):
         self.src_reads = ['x','y','z','h','m','rho','p']
@@ -82,10 +83,12 @@ cdef class SPHPressureGradient(SPHFunctionParticle):
 
         if self.bonnet_and_lok_correction:
             self.bonnet_and_lok_gradient_correction(dest_pid, &grad)
-            
+        
         nr[0] += temp*grad.x
-        nr[1] += temp*grad.y
-        nr[2] += temp*grad.z
+        if self.num_outputs > 1:
+            nr[1] += temp*grad.y
+            if self.num_outputs > 2:
+                nr[2] += temp*grad.z
 
     def cl_eval(self, object queue, object context):
 
@@ -113,7 +116,7 @@ cdef class MomentumEquation(SPHFunctionParticle):
 
     def __init__(self, ParticleArray source, ParticleArray dest, 
                  bint setup_arrays=True, alpha=1, beta=1, gamma=1.4, 
-                 eta=0.1, **kwargs):
+                 eta=0.1, int dim=3, **kwargs):
 
         SPHFunctionParticle.__init__(self, source, dest, setup_arrays,
                                      **kwargs)
@@ -128,6 +131,7 @@ cdef class MomentumEquation(SPHFunctionParticle):
 
         self.cl_kernel_src_file = "pressure_funcs.cl"
         self.cl_kernel_function_name = "MomentumEquation"
+        self.num_outputs = dim
 
     def set_src_dst_reads(self):
         self.src_reads = ['x','y','z','h','m','rho','p','u','v','w','cs']
@@ -235,8 +239,10 @@ cdef class MomentumEquation(SPHFunctionParticle):
             self.bonnet_and_lok_gradient_correction(dest_pid, &grad)
 
         nr[0] += tmp*grad.x
-        nr[1] += tmp*grad.y
-        nr[2] += tmp*grad.z
+        if self.num_outputs > 1:
+            nr[1] += tmp*grad.y
+            if self.num_outputs > 2:
+                nr[2] += tmp*grad.z
 
     def cl_eval(self, object queue, object context):
 
diff --git a/source/pysph/sph/funcs/xsph_funcs.pyx b/source/pysph/sph/funcs/xsph_funcs.pyx
index 387a787..c5373a0 100644
--- a/source/pysph/sph/funcs/xsph_funcs.pyx
+++ b/source/pysph/sph/funcs/xsph_funcs.pyx
@@ -11,7 +11,7 @@ cdef class XSPHCorrection(CSPHFunctionParticle):
     #Defined in the .pxd file
 
     def __init__(self, ParticleArray source, ParticleArray dest, 
-                 bint setup_arrays=True, double eps = 0.5, **kwargs):
+                 bint setup_arrays=True, double eps = 0.5, int dim=3, **kwargs):
 
         CSPHFunctionParticle.__init__(self, source, dest, setup_arrays,
                                      **kwargs)
@@ -22,6 +22,7 @@ cdef class XSPHCorrection(CSPHFunctionParticle):
 
         self.cl_kernel_src_file = "xsph_funcs.cl"
         self.cl_kernel_function_name = "XSPHCorrection"
+        self.num_outputs = dim
 
     def set_src_dst_reads(self):
         self.src_reads = ['x','y','z','h','m','rho','u','v','w']
@@ -81,8 +82,10 @@ cdef class XSPHCorrection(CSPHFunctionParticle):
         temp = mb * w/rhoab
 
         nr[0] += temp*Vba.x*self.eps
-        nr[1] += temp*Vba.y*self.eps
-        nr[2] += temp*Vba.z*self.eps
+        if self.num_outputs > 1:
+            nr[1] += temp*Vba.y*self.eps
+            if self.num_outputs > 2:
+                nr[2] += temp*Vba.z*self.eps
         
 ##########################################################################
 
@@ -95,6 +98,11 @@ cdef class XSPHDensityRate(SPHFunctionParticle):
 
     #Defined in the .pxd file
     #cdef DoubleArray s_ubar, s_vbar, s_wbar
+    
+    def __init__(self, ParticleArray source, ParticleArray dest,
+                 bint setup_arrays=True, *args, **kwargs):
+        SPHFunctionParticle.__init__(self, source, dest, setup_arrays, **kwargs)
+        self.num_outputs = 1
 
     cpdef setup_arrays(self):
         """ Setup the arrays required to read data from source and dest. """
diff --git a/source/pysph/sph/sph_calc.pxd b/source/pysph/sph/sph_calc.pxd
index 4a84d0d..0839665 100644
--- a/source/pysph/sph/sph_calc.pxd
+++ b/source/pysph/sph/sph_calc.pxd
@@ -78,6 +78,5 @@ cdef class SPHCalc:
     cdef setup_internals(self)
     cpdef check_internals(self)
 
-    cdef reset_output_arrays(self, DoubleArray output1, DoubleArray output2,
-                             DoubleArray output3)
+    cdef reset_output_array(self, DoubleArray output)
 
diff --git a/source/pysph/sph/sph_calc.pyx b/source/pysph/sph/sph_calc.pyx
index 925ca35..7de97f3 100644
--- a/source/pysph/sph/sph_calc.pyx
+++ b/source/pysph/sph/sph_calc.pyx
@@ -5,7 +5,6 @@ This module provides the SPHCalc class, which does the actual SPH summation.
     
 """
 
-#include "stdlib.pxd"
 from libc.stdlib cimport *
 
 cimport numpy
@@ -172,7 +171,10 @@ cdef class SPHCalc:
                 msg = 'SPHFunction.source not same as'
                 msg += ' SPHCalc.sources[%d]'%(i)
                 raise ValueError, msg
-
+            if funcs[i].num_outputs > len(self.updates):
+                raise ValueError, ('Number of updates not same as num_outputs'
+                       '; required %d, got %d %s for func %s'%(funcs[i].num_outputs,
+                                       len(self.updates), self.updates, funcs[i]))
             # not valid for SPHFunction
             #if funcs[i].dest != self.dest:
             #    msg = 'SPHFunction.dest not same as'
@@ -228,8 +230,11 @@ cdef class SPHCalc:
         cdef DoubleArray output1 = self.dest.get_carray(output_array1)
         cdef DoubleArray output2 = self.dest.get_carray(output_array2)
         cdef DoubleArray output3 = self.dest.get_carray(output_array3)
+        
+        if output1 is not None: self.reset_output_array(output1)
+        if output2 is not None: self.reset_output_array(output2)
+        if output3 is not None: self.reset_output_array(output3)
 
-        self.reset_output_arrays(output1, output2, output3)
         self.sph_array(output1, output2, output3, exclude_self)
 
         # call an update on the particles if the destination pa is dirty
@@ -264,14 +269,11 @@ cdef class SPHCalc:
 
             func.eval(self.kernel, output1, output2, output3)
 
-    cdef reset_output_arrays(self, DoubleArray output1, DoubleArray output2,
-                             DoubleArray output3):
+    cdef reset_output_array(self, DoubleArray output):
 
         cdef int i
-        for i in range(output1.length):
-            output1.data[i] = 0.0
-            output2.data[i] = 0.0
-            output3.data[i] = 0.0
+        for i in range(output.length):
+            output.data[i] = 0.0
 
 #############################################################################
 
diff --git a/source/pysph/sph/sph_func.pxd b/source/pysph/sph/sph_func.pxd
index 3d9ec62..9575a77 100755
--- a/source/pysph/sph/sph_func.pxd
+++ b/source/pysph/sph/sph_func.pxd
@@ -13,7 +13,7 @@ from pysph.base.nnps cimport FixedDestNbrParticleLocator
 cdef class SPHFunction:
     cdef public ParticleArray source, dest
     cdef public FixedDestNbrParticleLocator nbr_locator
-    cdef readonly num_outputs
+    cdef public num_outputs
 
     cdef public str name, id
     cdef public str tag
@@ -65,8 +65,6 @@ cdef class SPHFunction:
     cpdef setup_arrays(self)
     cpdef setup_iter_data(self)
 
-    cpdef int output_fields(self) except -1
-    
     cpdef eval(self, KernelBase kernel, DoubleArray output1,
                DoubleArray output2, DoubleArray output3)
 
diff --git a/source/pysph/sph/sph_func.pyx b/source/pysph/sph/sph_func.pyx
index b5dd11b..bac3e9f 100755
--- a/source/pysph/sph/sph_func.pyx
+++ b/source/pysph/sph/sph_func.pyx
@@ -207,15 +207,31 @@ cdef class SPHFunction:
 
         self.setup_iter_data()
         cdef size_t np = self.dest.get_number_of_particles()
-
-        for i in range(np):
-            if tag_arr.data[i] == LocalReal:
-                self.eval_single(i, kernel, result)
-                output1.data[i] += result[0]
-                output2.data[i] += result[1]
-                output3.data[i] += result[2]
-            else:
-                output1.data[i] = output2.data[i] = output3.data[i] = 0
+        
+        if self.num_outputs == 3:
+            for i in range(np):
+                if tag_arr.data[i] == LocalReal:
+                    self.eval_single(i, kernel, result)
+                    output1.data[i] += result[0]
+                    output2.data[i] += result[1]
+                    output3.data[i] += result[2]
+                else:
+                    output1.data[i] = output2.data[i] = output3.data[i] = 0
+        elif self.num_outputs == 2:
+            for i in range(np):
+                if tag_arr.data[i] == LocalReal:
+                    self.eval_single(i, kernel, result)
+                    output1.data[i] += result[0]
+                    output2.data[i] += result[1]
+                else:
+                    output1.data[i] = output2.data[i] = 0
+        elif self.num_outputs == 1:
+            for i in range(np):
+                if tag_arr.data[i] == LocalReal:
+                    self.eval_single(i, kernel, result)
+                    output1.data[i] += result[0]
+                else:
+                    output1.data[i] = 0
 
     cdef void eval_single(self, size_t dest_pid, KernelBase kernel,
                           double * result):
@@ -226,9 +242,6 @@ cdef class SPHFunction:
         """
         raise NotImplementedError, 'SPHFunction.eval_single()'
     	
-    cpdef int output_fields(self) except - 1:
-        return self.num_outputs
-    
     cpdef setup_iter_data(self):
         """ setup operations performed in each iteration
         
@@ -311,13 +324,13 @@ cdef class SPHFunction:
         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`
+        from the destination 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!")
+        pass
 
 ################################################################################
 # `SPHFunctionParticle` class.
diff --git a/source/pysph/sph/tests/test_boundary_funcs.py b/source/pysph/sph/tests/test_boundary_funcs.py
index 0c95134..3d4131a 100644
--- a/source/pysph/sph/tests/test_boundary_funcs.py
+++ b/source/pysph/sph/tests/test_boundary_funcs.py
@@ -55,7 +55,7 @@ def setUp(self):
         rhof = numpy.ones_like(xf)
 
         self.fluid = base.get_particle_array(x=xf, y=yf, h=hf, m=mf, rho=rhof,
-                                             cs=cs,
+                                             cs=cs, tmpx=mf, tmpy=mf, tmpz=mf,
                                              name='fluid', type=Fluid)
 
 
@@ -76,7 +76,7 @@ def setUp(self):
                 
                 sph.MonaghanBoundaryForce.withargs(delp=dp),
                 from_types = [Solid], on_types=[Fluid],
-                updates=['u','v'], id='boundary')
+                updates=['u','v','w'], id='boundary')
                                   
                                   )
                              
diff --git a/source/pysph/sph/tests/test_sph_calc.py b/source/pysph/sph/tests/test_sph_calc.py
index cc7f9fa..1b162c9 100644
--- a/source/pysph/sph/tests/test_sph_calc.py
+++ b/source/pysph/sph/tests/test_sph_calc.py
@@ -26,7 +26,7 @@ def test_sph_calc():
     z = numpy.array([0,])
     h = numpy.ones_like(x)
 
-    pa = base.get_particle_array(name="test", x=x, y=y, z=z,h=h)
+    pa = base.get_particle_array(name="test", x=x, y=y, z=z,h=h, _tmpx=z, _tmpy=z, _tmpz=z)
     particles = base.Particles(arrays=[pa,])
     kernel = base.CubicSplineKernel(dim=1)
 
diff --git a/source/pysph/sph/tests/test_update_smoothing.py b/source/pysph/sph/tests/test_update_smoothing.py
index 8f90841..94a49ca 100644
--- a/source/pysph/sph/tests/test_update_smoothing.py
+++ b/source/pysph/sph/tests/test_update_smoothing.py
@@ -51,7 +51,7 @@ def setUp(self):
         h0 = dx
 
         pa = base.get_particle_array(name="test", type=Fluid,
-                                     x=x, h=h, rho=rho,  m=m)
+                                     x=x, h=h, rho=rho, m=m, tmp=m)
 
         particles = base.Particles(arrays=[pa], variable_h=True)
 
@@ -192,6 +192,10 @@ class TestUpdateSmoothing(sph.SPHFunctionParticle):
     The smoothing length of the fifth particle is doubled i.e. h[5] *= 2
 
     """
+    def __init__(self, source, dest, setup_arrays=True):
+        sph.SPHFunctionParticle.__init__(self, source, dest, setup_arrays)
+        self.num_outputs = 1
+        self.dst_reads = ['h']
 
     def eval(self, kernel, output1, output2, output3):
         for i in range(output1.length):
@@ -209,24 +213,22 @@ def setUp(self):
         self.h0 = dx
 
         self.pa = pa = base.get_particle_array(name="test", type=Fluid,
-                                     x=x, h=h, rho=rho,  m=m)
+                                     x=x, h=h, rho=rho, m=m, tmp=m)
 
         self.particles = base.Particles(arrays=[pa], variable_h=True)
 
         # smoothing length update operation
 
         self.adke = solver.SPHOperation(
-
-            TestUpdateSmoothing.withargs( h0=self.h0, k=1.0, eps=0.0 ),
+            TestUpdateSmoothing,
             from_types=[Fluid], on_types=[Fluid],
-            updates=["h"], id="pilotrho" )
-
+            updates=["tmp"], # updates='h' will first reset array h to 0
+            id="pilotrho" )
         self.kernel = base.CubicSplineKernel(dim=1)
     
     def test_solve(self):
 
-        for integrator in (solver.EulerIntegrator,solver.RK2Integrator,
-                           solver.RK4Integrator):
+        for integrator in [solver.EulerIntegrator]:
 
             self.solver = solver.Solver(self.kernel.dim, integrator)
             self.solver.add_operation(self.adke)
@@ -241,7 +243,7 @@ def test_solve(self):
             h_expect = self.pa.h[5] * 2**self.solver.integrator.nsteps
 
             self.solver.solve(False)
-            self.assertAlmostEqual(self.pa.h[5], h_expect)
+            self.assertAlmostEqual(self.pa.tmp[5], h_expect)
     
 if __name__ == "__main__":
     unittest.main()