Merge pull request #194 from adityapb/master

 Add update_minmax_cl method to calculate min and max in single scan

pysph/base/device_helper.py

@@ -1,7 +1,8 @@
 from __future__ import print_function
 import logging
 import numpy as np
-from pytools import memoize_method
+from pytools import memoize, memoize_method
+import mako.template as mkt
 
 from compyle.config import get_config
 from compyle.array import get_backend, wrap_array, Array
@@ -16,6 +17,85 @@
 logger = logging.getLogger()
 
 
+minmax_tpl = """//CL//
+
+    ${dtype} mmc_neutral()
+    {
+        ${dtype} result;
+
+        % for prop in prop_names:
+        % if not only_max:
+        result.cur_min_${prop} = INFINITY;
+        % endif
+        % if not only_min:
+        result.cur_max_${prop} = -INFINITY;
+        % endif
+        % endfor
+
+        return result;
+    }
+
+    ${dtype} mmc_from_scalar(${args})
+    {
+        ${dtype} result;
+
+        % for prop in prop_names:
+        % if not only_max:
+        result.cur_min_${prop} = ${prop};
+        % endif
+        % if not only_min:
+        result.cur_max_${prop} = ${prop};
+        % endif
+        % endfor
+
+        return result;
+    }
+
+    ${dtype} agg_mmc(${dtype} a, ${dtype} b)
+    {
+        ${dtype} result = a;
+
+        % for prop in prop_names:
+        % if not only_max:
+        if (b.cur_min_${prop} < result.cur_min_${prop})
+            result.cur_min_${prop} = b.cur_min_${prop};
+        % endif
+        % if not only_min:
+        if (b.cur_max_${prop} > result.cur_max_${prop})
+            result.cur_max_${prop} = b.cur_max_${prop};
+        % endif
+        % endfor
+
+        return result;
+    }
+
+    """
+
+
+def minmax_collector_key(device, dtype, props, name, *args):
+    return (device, dtype, tuple(props), name)
+
+
+@memoize(key=minmax_collector_key)
+def make_collector_dtype(device, dtype, props, name, only_min, only_max):
+    fields = [("pad", np.int32)]
+
+    for prop in props:
+        if not only_min:
+            fields.append(("cur_max_%s" % prop, dtype))
+        if not only_max:
+            fields.append(("cur_min_%s" % prop, dtype))
+
+    custom_dtype = np.dtype(fields)
+
+    from pyopencl.tools import get_or_register_dtype, match_dtype_to_c_struct
+
+    custom_dtype, c_decl = match_dtype_to_c_struct(device, name, custom_dtype)
+    custom_dtype = get_or_register_dtype(name, custom_dtype)
+
+    return custom_dtype, c_decl
+
+
 class ExtractParticles(Template):
     def __init__(self, name, prop_names):
         super(ExtractParticles, self).__init__(name=name)
@@ -54,8 +134,8 @@ class DeviceHelper(object):
     def __init__(self, particle_array, backend=None):
         self.backend = get_backend(backend)
         self._particle_array = pa = particle_array
-        use_double = get_config().use_double
-        self._dtype = np.float64 if use_double else np.float32
+        self.use_double = get_config().use_double
+        self._dtype = np.float64 if self.use_double else np.float32
         self.num_real_particles = pa.num_real_particles
         self._data = {}
         self.properties = []
@@ -156,6 +236,74 @@ def max(self, arg):
         return float(array.maximum(getattr(self, arg),
                                    backend=self.backend))
 
+    @memoize(key=lambda *args: (args[-2], args[-1]))
+    def _get_minmax_kernel(self, ctx, dtype, mmc_dtype, prop_names,
+                           only_min, only_max, name, mmc_c_decl):
+        tpl_args = ", ".join(
+            ["%(dtype)s %(prop)s" % {'dtype': dtype, 'prop': prop}
+                for prop in prop_names]
+        )
+
+        mmc_preamble = mmc_c_decl + minmax_tpl
+        preamble = mkt.Template(text=mmc_preamble).render(
+            args=tpl_args, prop_names=prop_names, dtype=name,
+            only_min=only_min, only_max=only_max
+        )
+
+        knl_args = ", ".join(
+            ["__global %(dtype)s *%(prop)s" % {'dtype': dtype, 'prop': prop}
+                for prop in prop_names]
+        )
+
+        map_args = ", ".join(
+            ["%(prop)s[i]" % {'dtype': dtype, 'prop': prop}
+                for prop in prop_names]
+        )
+
+        from pyopencl.reduction import ReductionKernel
+
+        knl = ReductionKernel(
+            ctx, mmc_dtype, neutral="mmc_neutral()",
+            reduce_expr="agg_mmc(a, b)",
+            map_expr="mmc_from_scalar(%s)" % map_args,
+            arguments=knl_args,
+            preamble=preamble
+        )
+
+        return knl
+
+    def update_minmax_cl(self, props, only_min=False, only_max=False):
+        if self.backend != 'opencl':
+            raise ValueError('''Optimized minmax update only supported
+                             using opencl backend''')
+        if only_min and only_max:
+            raise ValueError("Only one of only_min and only_max can be True")
+
+        dtype = 'double' if self.use_double else 'float'
+        op = 'min' if not only_max else ''
+        op += 'max' if not only_min else ''
+        name = "%s_collector_%s" % (op, ''.join([dtype] + props))
+
+        from compyle.opencl import get_context
+        ctx = get_context()
+
+        mmc_dtype, mmc_c_decl = make_collector_dtype(ctx.devices[0],
+                                                     self._dtype, props, name,
+                                                     only_min, only_max)
+        knl = self._get_minmax_kernel(ctx, dtype, mmc_dtype, props,
+                                      only_min, only_max, name, mmc_c_decl)
+
+        args = [getattr(self, prop).dev for prop in props]
+
+        result = knl(*args).get()
+
+        for prop in props:
+            proparr = self._data[prop]
+            if not only_max:
+                proparr.minimum = result["cur_min_%s" % prop]
+            if not only_min:
+                proparr.maximum = result["cur_max_%s" % prop]
+
     def update_min_max(self, props=None):
         """Update the min,max values of all properties """
         props = props if props else self.properties
@@ -583,7 +731,8 @@ def empty_clone(self, props=None):
         result_array.set_output_arrays(output_arrays)
         return result_array
 
-    def extract_particles(self, indices, dest_array=None, align=True, props=None):
+    def extract_particles(self, indices, dest_array=None,
+                          align=True, props=None):
         """Create new particle array for particles with given indices
 
         Parameters

pysph/base/gpu_domain_manager.py

@@ -71,7 +71,7 @@ def _compute_cell_size_for_binning(self):
 
         for pa_wrapper in self.pa_wrappers:
             h = pa_wrapper.pa.gpu.get_device_array('h')
-            h.update_min_max()
+            pa_wrapper.pa.gpu.update_minmax_cl(['h'])
 
             _hmax = h.maximum
             _hmin = h.minimum

pysph/base/gpu_nnps_base.pyx

@@ -271,9 +271,8 @@ cdef class GPUNNPS(NNPSBase):
             y = pa_wrapper.pa.gpu.get_device_array('y')
             z = pa_wrapper.pa.gpu.get_device_array('z')
 
-            x.update_min_max()
-            y.update_min_max()
-            z.update_min_max()
+            pa_wrapper.pa.gpu.update_minmax_cl(['x', 'y', 'z'])
+
             # find min and max of variables
             xmax = np.maximum(x.maximum, xmax)
             ymax = np.maximum(y.maximum, ymax)

pysph/base/tests/test_device_helper.py

@@ -394,3 +394,76 @@ def test_extract_particles(self, backend):
 
         # Then
         assert result_pa.gpu.get_number_of_particles() == 2
+
+    def test_update_minmax_cl(self):
+        backend = 'opencl'
+        check_import(backend)
+        self.setup()
+
+        # Given
+        x = [0.0, -1.0, 2.0, 3.0]
+        y = [0.0, 1.0, -2.0, 3.0]
+        z = [0.0, 1.0, 2.0, -3.0]
+        h = [4.0, 1.0, 2.0, 3.0]
+
+        pa = get_particle_array(x=x, y=y, z=z, h=h)
+        h = DeviceHelper(pa, backend=backend)
+        pa.set_device_helper(h)
+
+        # When
+        h.update_minmax_cl(['x', 'y', 'z', 'h'])
+
+        # Then
+        assert h.x.minimum == -1.0
+        assert h.x.maximum == 3.0
+
+        assert h.y.minimum == -2.0
+        assert h.y.maximum == 3.0
+
+        assert h.z.minimum == -3.0
+        assert h.z.maximum == 2.0
+
+        assert h.h.minimum == 1.0
+        assert h.h.maximum == 4.0
+
+        # When
+        h.x.maximum, h.x.minimum = 100., 100.
+        h.y.maximum, h.y.minimum = 100., 100.
+        h.z.maximum, h.z.minimum = 100., 100.
+        h.h.maximum, h.h.minimum = 100., 100.
+
+        h.update_minmax_cl(['x', 'y', 'z', 'h'], only_min=True)
+
+        # Then
+        assert h.x.minimum == -1.0
+        assert h.x.maximum == 100.0
+
+        assert h.y.minimum == -2.0
+        assert h.y.maximum == 100.0
+
+        assert h.z.minimum == -3.0
+        assert h.z.maximum == 100.0
+
+        assert h.h.minimum == 1.0
+        assert h.h.maximum == 100.0
+
+        # When
+        h.x.maximum, h.x.minimum = 100., 100.
+        h.y.maximum, h.y.minimum = 100., 100.
+        h.z.maximum, h.z.minimum = 100., 100.
+        h.h.maximum, h.h.minimum = 100., 100.
+
+        h.update_minmax_cl(['x', 'y', 'z', 'h'], only_max=True)
+
+        # Then
+        assert h.x.minimum == 100.0
+        assert h.x.maximum == 3.0
+
+        assert h.y.minimum == 100.0
+        assert h.y.maximum == 3.0
+
+        assert h.z.minimum == 100.0
+        assert h.z.maximum == 2.0
+
+        assert h.h.minimum == 100.0
+        assert h.h.maximum == 4.0

pysph/sph/integrator.py

@@ -54,13 +54,20 @@ def _get_dt_adapt_factors(self):
         a_eval = self.acceleration_evals[0]
         factors = [-1.0, -1.0, -1.0]
         for pa in a_eval.particle_arrays:
+            prop_names = []
             for i, name in enumerate(('dt_cfl', 'dt_force', 'dt_visc')):
                 if name in pa.properties:
                     if pa.gpu:
-                        max_val = pa.gpu.max(name)
+                        prop_names.append(name)
                     else:
                         max_val = np.max(pa.get(name))
-                    factors[i] = max(factors[i], max_val)
+                        factors[i] = max(factors[i], max_val)
+            if pa.gpu:
+                pa.gpu.update_minmax_cl(prop_names, only_max=True)
+                for i, name in enumerate(('dt_cfl', 'dt_force', 'dt_visc')):
+                    if name in pa.properties:
+                        max_val = getattr(pa.gpu, name).maximum
+                        factors[i] = max(factors[i], max_val)
         cfl_f, force_f, visc_f = factors
         return cfl_f, force_f, visc_f
 
@@ -101,8 +108,6 @@ def compute_h_minimum(self):
             else:
                 h = pa.get_carray('h')
 
-            h.update_min_max()
-
             if h.minimum < hmin:
                 hmin = h.minimum