Fix CRKSPH OpenCL kernel compilation

pysph/sph/gas_dynamics/magma2.py

@@ -29,7 +29,7 @@
 from pysph.sph.integrator_step import IntegratorStep
 from pysph.sph.scheme import Scheme, add_bool_argument
 from pysph.sph.wc.linalg import (augmented_matrix, gj_solve, identity,
-                                 mat_mult, mat_vec_mult, dot)
+                                 mat_mult, mat_vec_mult, dot_product)
 
 GHOST_TAG = get_ghost_tag()
 
@@ -1119,7 +1119,7 @@ def __init__(self, dest, sources, dim, fkern, eta_crit=0.3, eta_fold=0.2,
         super().__init__(dest, sources)
 
     def _get_helpers_(self):
-        return [dot]
+        return [dot_product]
 
     def initialize(self, d_idx, d_au, d_av, d_aw, d_ae):
         d_au[d_idx] = 0.0
@@ -1155,8 +1155,8 @@ def loop(self, d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au,
             etaj[row] = XIJ[row] / hj
 
         # Limiter
-        etaisq = dot(etai, etai, dim)
-        etajsq = dot(etaj, etaj, dim)
+        etaisq = dot_product(etai, etai, dim)
+        etajsq = dot_product(etaj, etaj, dim)
         etaij = sqrt(min(etaisq, etajsq))
 
         aanum = 0.0
@@ -1217,8 +1217,8 @@ def loop(self, d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au,
 
         # Artificial viscosity
         vsigng = sqrt(abs(d_p[d_idx] - s_p[s_idx]) * RHOIJ1)
-        mui = min(0.0, dot(vij, etai, dim) / (etaisq + epssq))
-        muj = min(0.0, dot(vij, etaj, dim) / (etajsq + epssq))
+        mui = min(0.0, dot_product(vij, etai, dim) / (etaisq + epssq))
+        muj = min(0.0, dot_product(vij, etaj, dim) / (etajsq + epssq))
         qi = d_rho[d_idx] * mui * (-d_alpha[d_idx] * d_cs[d_idx] + beta * mui)
         qj = s_rho[s_idx] * muj * (-s_alpha[s_idx] * s_cs[s_idx] + beta * muj)
         pi = d_p[d_idx] + qi
@@ -1233,7 +1233,7 @@ def loop(self, d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au,
         d_aw[d_idx] -= mjpibyrhoisq * DWI[2] + mjpjbyrhojsq * DWJ[2]
 
         # Accelerations for the thermal energy
-        vijdotdwi = dot(VIJ, DWI, dim)
+        vijdotdwi = dot_product(VIJ, DWI, dim)
         d_ae[d_idx] += mjpibyrhoisq * vijdotdwi
 
         # artificial conduction
@@ -1273,8 +1273,8 @@ def loop(self, d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au,
             etaj[row] = XIJ[row] / hj
 
         # Limiter
-        etaisq = dot(etai, etai, dim)
-        etajsq = dot(etaj, etaj, dim)
+        etaisq = dot_product(etai, etai, dim)
+        etajsq = dot_product(etaj, etaj, dim)
         etaij = sqrt(min(etaisq, etajsq))
 
         aanum = 0.0
@@ -1340,8 +1340,8 @@ def loop(self, d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au,
 
         # Artificial viscosity
         vsigng = sqrt(abs(d_p[d_idx] - s_p[s_idx]) * RHOIJ1)
-        mui = min(0.0, dot(vij, etai, dim) / (etaisq + epssq))
-        muj = min(0.0, dot(vij, etaj, dim) / (etajsq + epssq))
+        mui = min(0.0, dot_product(vij, etai, dim) / (etaisq + epssq))
+        muj = min(0.0, dot_product(vij, etaj, dim) / (etajsq + epssq))
         qi = d_rho[d_idx] * mui * (-d_alpha[d_idx] * d_cs[d_idx] + beta * mui)
         qj = s_rho[s_idx] * muj * (-s_alpha[s_idx] * s_cs[s_idx] + beta * muj)
         pi = d_p[d_idx] + qi
@@ -1356,7 +1356,7 @@ def loop(self, d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au,
         d_aw[d_idx] -= mjpibyrhoisq * gmi[2] + mjpjbyrhojsq * gmj[2]
 
         # Accelerations for the thermal energy
-        vijdotdwi = dot(VIJ, gmi, dim)
+        vijdotdwi = dot_product(VIJ, gmi, dim)
         d_ae[d_idx] += mjpibyrhoisq * vijdotdwi
 
         # artificial conduction
@@ -1395,8 +1395,8 @@ def loop(self, d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au,
             gmij[row] = 0.0
 
         # Limiter
-        etaisq = dot(etai, etai, dim)
-        etajsq = dot(etaj, etaj, dim)
+        etaisq = dot_product(etai, etai, dim)
+        etajsq = dot_product(etaj, etaj, dim)
         etaij = sqrt(min(etaisq, etajsq))
 
         aanum = 0.0
@@ -1464,8 +1464,8 @@ def loop(self, d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au,
 
         # Artificial viscosity
         vsigng = sqrt(abs(d_p[d_idx] - s_p[s_idx]) * RHOIJ1)
-        mui = min(0.0, dot(vij, etai, dim) / (etaisq + epssq))
-        muj = min(0.0, dot(vij, etaj, dim) / (etajsq + epssq))
+        mui = min(0.0, dot_product(vij, etai, dim) / (etaisq + epssq))
+        muj = min(0.0, dot_product(vij, etaj, dim) / (etajsq + epssq))
         qi = d_rho[d_idx] * mui * (-d_alpha[d_idx] * d_cs[d_idx] + beta * mui)
         qj = s_rho[s_idx] * muj * (-s_alpha[s_idx] * s_cs[s_idx] + beta * muj)
         pi = d_p[d_idx] + qi
@@ -1480,7 +1480,7 @@ def loop(self, d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au,
         d_aw[d_idx] -= comn * gmij[2]
 
         # Accelerations for the thermal energy
-        vijdotgmij = dot(VIJ, gmij, dim)
+        vijdotgmij = dot_product(VIJ, gmij, dim)
         d_ae[d_idx] -= self.alphac * s_m[
             s_idx] * vsigng * eij * normgmij * RHOIJ1
         d_ae[d_idx] += s_m[s_idx] * pi * invrhosq * vijdotgmij
@@ -1496,14 +1496,14 @@ def __init__(self, dest, sources, dim):
         super().__init__(dest, sources)
 
     def _get_helpers_(self):
-        return [dot]
+        return [dot_product]
 
     def initialize(self, d_idx, d_tilmu):
         d_tilmu[d_idx] = -INFINITY
 
     def loop(self, d_tilmu, d_idx, d_h, VIJ, XIJ, R2IJ):
         d_tilmu[d_idx] = max(d_tilmu[d_idx],
-                             d_h[d_idx] * dot(VIJ, XIJ, self.dim) / (
+                             d_h[d_idx] * dot_product(VIJ, XIJ, self.dim) / (
                                      R2IJ + 0.01))
 
 

pysph/sph/tests/test_crksph.py

@@ -0,0 +1,74 @@
+import numpy as np
+import pytest
+
+from compyle.config import use_config
+from pysph.base.kernels import QuinticSpline
+from pysph.base.nnps import LinkedListNNPS
+from pysph.base.utils import get_particle_array
+from pysph.sph.acceleration_eval import make_acceleration_evals
+from pysph.sph.sph_compiler import SPHCompiler
+from pysph.sph.wc.crksph import CRKSPHScheme
+
+
+def evaluate_crksph(backend):
+    x, y = np.mgrid[-0.6:0.6:7j, -0.6:0.6:7j]
+    x, y = x.ravel(), y.ravel()
+    pa = get_particle_array(
+        name='fluid', x=x, y=y, h=0.3, m=0.04,
+        rho=1.0 + 0.1*x, e=2.0 + 0.2*y,
+        u=-0.3*x + 0.1*y, v=0.2*x - 0.4*y
+    )
+    scheme = CRKSPHScheme(
+        fluids=['fluid'], dim=2, rho0=1.0, c0=1.0, nu=0.0,
+        h0=0.3, p0=1.0, gamma=1.4
+    )
+    scheme.setup_properties([pa])
+    pa.u0[:] = pa.u
+    pa.v0[:] = pa.v
+    kernel = QuinticSpline(dim=2)
+    evaluators = make_acceleration_evals(
+        [pa], scheme.get_equations(), kernel, backend=backend
+    )
+    SPHCompiler(evaluators, integrator=None).compile()
+    if backend == 'opencl':
+        from pysph.base.gpu_nnps import ZOrderGPUNNPS
+        nnps = ZOrderGPUNNPS(
+            dim=2, particles=[pa], radius_scale=kernel.radius_scale,
+            cache=True, backend=backend
+        )
+    else:
+        nnps = LinkedListNNPS(
+            dim=2, particles=[pa], radius_scale=kernel.radius_scale,
+            cache=True
+        )
+    for evaluator in evaluators:
+        evaluator.set_nnps(nnps)
+        evaluator.compute(0.0, 0.001)
+    props = ('rho', 'p', 'cs', 'ai', 'bi', 'gradai', 'gradbi',
+             'gradv', 'au', 'av', 'ae')
+    if backend == 'opencl':
+        pa.gpu.pull(*props)
+    return {name: pa.get(name).copy() for name in props}
+
+
+@pytest.mark.parametrize('use_double', [False, True])
+def test_crksph_opencl_matches_cython(use_double):
+    pytest.importorskip('pyopencl')
+    pytest.importorskip('pysph.base.gpu_nnps')
+    from compyle.opencl import get_context
+    devices = get_context().devices
+    if use_double and not all(d.double_fp_config for d in devices):
+        pytest.skip('OpenCL device does not support double precision')
+
+    # Exercise both stages, including the correction solve, viscosity limiter,
+    # energy equation and speed of sound used by the Kelvin-Helmholtz example.
+    with use_config(use_opencl=False, use_double=use_double):
+        expected = evaluate_crksph('cython')
+        actual = evaluate_crksph('opencl')
+    for name in expected:
+        assert np.isfinite(actual[name]).all(), name
+        rtol, atol = (1e-9, 1e-10) if use_double else (2e-4, 2e-5)
+        np.testing.assert_allclose(
+            actual[name], expected[name], rtol=rtol, atol=atol,
+            err_msg=name
+        )

pysph/sph/tests/test_linalg.py

@@ -1,5 +1,5 @@
 from pysph.sph.wc.linalg import (
-    augmented_matrix, gj_solve, mat_mult, mat_vec_mult
+    augmented_matrix, dot, dot_product, gj_solve, mat_mult, mat_vec_mult
 )
 import numpy as np
 import unittest
@@ -15,6 +15,14 @@ def gj_solve_helper(a, b, n):
 
 class TestLinalg(unittest.TestCase):
 
+    def test_dot_product(self):
+        a = [1.5, -2.0, 3.0]
+        b = [-4.0, 0.5, 2.0]
+        for n in (1, 2, 3):
+            expected = np.dot(a[:n], b[:n])
+            self.assertAlmostEqual(dot_product(a, b, n), expected)
+            self.assertAlmostEqual(dot(a, b, n), expected)
+
     def _to_array(self, x, shape=None):
         x = np.asarray(x)
         if shape:

pysph/sph/wc/crksph.py

@@ -17,7 +17,7 @@
 from compyle.api import declare
 from pysph.sph.equation import Equation, Group, MultiStageEquations
 from pysph.sph.wc.linalg import (
-    augmented_matrix, dot, gj_solve, identity, mat_vec_mult
+    augmented_matrix, dot_product, gj_solve, identity, mat_vec_mult
 )
 from pysph.sph.scheme import Scheme
 from pysph.base.utils import get_particle_array
@@ -35,7 +35,8 @@ def __init__(self, dest, sources, dim=2):
         super(CRKSPHPreStep, self).__init__(dest, sources)
 
     def _get_helpers_(self):
-        return [augmented_matrix, gj_solve, identity, dot, mat_vec_mult]
+        return [augmented_matrix, gj_solve, identity, dot_product,
+                mat_vec_mult]
 
     def loop_all(self, d_idx, d_x, d_y, d_z, d_h, s_x, s_y, s_z, s_h, s_m,
                  s_rho, SPH_KERNEL, NBRS, N_NBRS, d_ai, d_gradai, d_bi, s_V,
@@ -118,7 +119,7 @@ def loop_all(self, d_idx, d_x, d_y, d_z, d_h, s_x, s_y, s_z, s_h, s_m,
             mat_vec_mult(m2inv, m1, d, temp_vec)
 
             # Eq. 12.
-            ai = 1.0/(m0 - dot(temp_vec, m1, d))
+            ai = 1.0/(m0 - dot_product(temp_vec, m1, d))
             # Eq. 13.
             mat_vec_mult(m2inv, m1, d, bi)
             for gam in range(d):
@@ -548,7 +549,7 @@ def __init__(self, dest, sources, dim, gx=0.0, gy=0.0, gz=0.0, cl=2, cq=1,
         super(MomentumEquation, self).__init__(dest, sources)
 
     def _get_helpers_(self):
-        return [dot]
+        return [dot_product]
 
     def initialize(self, d_idx, d_au, d_av, d_aw):
         d_au[d_idx] = self.gx
@@ -593,10 +594,10 @@ def loop(self, d_idx, s_idx, d_m, s_m, d_rho, s_rho, d_p, s_p, d_cs, s_cs,
                 tmprj += s_gradv[d*d*s_idx + d*alp + bet] * XIJ[alp] * XIJ[bet]
         rij = tmpri/tmprj
 
-        tmprij = min(1, 4*rij/((1 + rij)*(1 + rij)))
-        phiij = max(0, tmprij)
+        tmprij = min(1.0, 4*rij/((1 + rij)*(1 + rij)))
+        phiij = max(0.0, tmprij)
 
-        tmpxij = dot(XIJ, XIJ, d)
+        tmpxij = dot_product(XIJ, XIJ, d)
         tmpxij2 = sqrt(tmpxij)
         etai_scalar = tmpxij2/hi
         etaj_scalar = tmpxij2/hj
@@ -617,11 +618,11 @@ def loop(self, d_idx, s_idx, d_m, s_m, d_rho, s_rho, d_p, s_p, d_cs, s_cs,
         uijhat[1] = d_v[d_idx] - s_v[s_idx] - 0.5*phiij * tmpdvxij[1]
         uijhat[2] = d_w[d_idx] - s_w[s_idx] - 0.5*phiij * tmpdvxij[2]
 
-        tmpmui = dot(uijhat, XIJ, d) / (tmpxij/hi + EPS * hi)
-        mui = min(0, tmpmui)
+        tmpmui = dot_product(uijhat, XIJ, d) / (tmpxij/hi + EPS * hi)
+        mui = min(0.0, tmpmui)
 
-        tmpmuj = dot(uijhat, XIJ, d) / (tmpxij/hi + EPS * hj)
-        muj = min(0, tmpmuj)
+        tmpmuj = dot_product(uijhat, XIJ, d) / (tmpxij/hi + EPS * hj)
+        muj = min(0.0, tmpmuj)
 
         Qi = rhoi * (-Cl*ci*mui + Cq*mui*mui)
         Qj = rhoj * (-Cl*cj*muj + Cq*muj*muj)
@@ -677,7 +678,7 @@ def __init__(self, dest, sources, dim, gamma, gx=0.0, gy=0.0, gz=0.0,
         super(EnergyEquation, self).__init__(dest, sources)
 
     def _get_helpers_(self):
-        return [dot]
+        return [dot_product]
 
     def initialize(self, d_idx, d_ae):
         d_ae[d_idx] = 0.0
@@ -720,10 +721,10 @@ def loop(self, d_idx, s_idx, d_au, d_av, d_aw, d_ae, s_au, s_av, s_aw,
                 tmprj += s_gradv[d*d*s_idx + d*alp + bet] * XIJ[alp] * XIJ[bet]
         rij = tmpri/tmprj
 
-        tmprij = min(1, 4*rij/((1 + rij)*(1 + rij)))
-        phiij = max(0, tmprij)
+        tmprij = min(1.0, 4*rij/((1 + rij)*(1 + rij)))
+        phiij = max(0.0, tmprij)
 
-        tmpxij = dot(XIJ, XIJ, d)
+        tmpxij = dot_product(XIJ, XIJ, d)
         tmpxij2 = sqrt(tmpxij)
         etai_scalar = tmpxij2/hi
         etaj_scalar = tmpxij2/hj
@@ -744,11 +745,11 @@ def loop(self, d_idx, s_idx, d_au, d_av, d_aw, d_ae, s_au, s_av, s_aw,
         uijhat[1] = d_v0[d_idx] - s_v0[s_idx] - 0.5*phiij * tmpdvxij[1]
         uijhat[2] = d_w0[d_idx] - s_w0[s_idx] - 0.5*phiij * tmpdvxij[2]
 
-        tmpmui = dot(uijhat, XIJ, d) / (tmpxij/hi + EPS * hi)
-        mui = min(0, tmpmui)
+        tmpmui = dot_product(uijhat, XIJ, d) / (tmpxij/hi + EPS * hi)
+        mui = min(0.0, tmpmui)
 
-        tmpmuj = dot(uijhat, XIJ, d) / (tmpxij/hi + EPS * hj)
-        muj = min(0, tmpmuj)
+        tmpmuj = dot_product(uijhat, XIJ, d) / (tmpxij/hi + EPS * hj)
+        muj = min(0.0, tmpmuj)
 
         Qi = rhoi * (-Cl*ci*mui + Cq*mui*mui)
         Qj = rhoj * (-Cl*cj*muj + Cq*muj*muj)
@@ -766,7 +767,7 @@ def loop(self, d_idx, s_idx, d_au, d_av, d_aw, d_ae, s_au, s_av, s_aw,
         delu[1] = s_v0[s_idx] + s_v[s_idx] - d_v0[d_idx] - d_v[d_idx]
         delu[2] = s_w0[s_idx] + s_w[s_idx] - d_w0[d_idx] - d_w[d_idx]
 
-        aeij = dot(delu, auij, d)
+        aeij = dot_product(delu, auij, d)
 
         si = d_p[d_idx]/((d_rho[d_idx])**gamma)
         sj = s_p[s_idx]/((s_rho[s_idx])**gamma)
@@ -807,7 +808,7 @@ def __init__(self, dest, sources=None, gamma=7.0):
         self.gamma = gamma
 
     def initialize(self, d_cs, d_idx, d_p, d_rho):
-        d_cs[d_idx] = (self.gamma * d_p[d_idx] / d_rho[d_idx])**0.5
+        d_cs[d_idx] = sqrt(self.gamma * d_p[d_idx] / d_rho[d_idx])
 
 
 class CRKSPHUpdateGhostProps(Equation):

pysph/sph/wc/linalg.py

@@ -13,14 +13,20 @@ def identity(a=[0.0, 0.0], n=3):
                 a[n*i + j] = 0.0
 
 
-def dot(a=[0.0, 0.0], b=[0.0, 0.0], n=3):
+def dot_product(a=[0.0, 0.0], b=[0.0, 0.0], n=3):
+    """Dot product of arrays, with a name distinct from OpenCL's builtin."""
     i = declare('int')
     result = 0.0
     for i in range(n):
         result += a[i]*b[i]
     return result
 
 
+# Backward-compatible name for Python callers. Transpiled equations should use
+# dot_product to avoid colliding with the OpenCL vector builtin.
+dot = dot_product
+
+
 def mat_mult(a=[1.0, 0.0], b=[1.0, 0.0], n=3, result=[0.0, 0.0]):
     """Multiply two square matrices (not element-wise).
 
@@ -131,10 +137,10 @@ def gj_solve(m=[1., 0.], n=3, nb=1, result=[0.0, 0.0]):
     rr, rrcol, rb, rbr, kup, kupr, kleft, kleftr = declare('int', 8)
     for rrcol in range(0, colrange):
         for rr in range(rrcol + 1, eqns):
-            dnr = float(m[nt*rrcol + rrcol])
+            dnr = m[nt*rrcol + rrcol]
             if abs(dnr) < 1e-12:
                 return 1.0
-            cc = -float(m[nt*rr + rrcol]) / dnr
+            cc = -m[nt*rr + rrcol] / dnr
             for j in range(augCol):
                 m[nt*rr + j] = m[nt*rr + j] + cc * m[nt*rrcol + j]