Merge pull request #195 from psnbaba/interpolator_fix

added gradient interpolation

pysph/sph/tests/test_linalg.py

@@ -7,7 +7,7 @@
 
 def gj_solve_helper(a, b, n):
     m = np.zeros((n, n+1)).ravel().tolist()
-    augmented_matrix(a, b, n, 1, m)
+    augmented_matrix(a, b, n, 1, n, m)
     result = [0.0]*n
     is_singular = gj_solve(m, n, 1, result)
     return is_singular, result
@@ -30,11 +30,40 @@ def test_augmented_matrix(self):
         expect[:, :3] = a
         expect[:, 3:] = b
         # When
-        augmented_matrix(a.ravel(), b.ravel(), 3, 2, res)
+        augmented_matrix(a.ravel(), b.ravel(), 3, 2, 3, res)
         res = self._to_array(res, (3, 5))
         # Then
         np.testing.assert_array_almost_equal(res, expect)
 
+    def test_augmented_matrix_with_lower_dimension(self):
+        # Given
+        a = np.random.random((3, 3))
+        b = np.random.random((3, 2))
+        res = np.zeros((3, 5)).ravel().tolist()
+        expect = np.zeros((2, 4))
+        expect[:, :2] = a[:2, :2]
+        expect[:, 2:] = b[:2, :]
+        expect.resize(3, 5)
+        # When
+        augmented_matrix(a.ravel(), b.ravel(), 2, 2, 3, res)
+        res = self._to_array(res, (3, 5))
+        # Then
+        np.testing.assert_array_almost_equal(res, expect)
+
+    def test_augmented_matrix_with_gjsolve_with_lower_dimension(self):
+        # Given
+        nmax = 3
+        mat = np.array([[7., 4., 2.], [8., 9., 4.], [1., 4., 10.]])
+        b = np.array([5., 4., 2.])
+        expect = np.linalg.solve(mat[:2, :2], b[:2])
+        augmat = np.zeros((3, 4)).ravel().tolist()
+        res = np.zeros(2).ravel().tolist()
+        # When
+        augmented_matrix(mat.ravel(), b.ravel(), 2, 1, nmax, augmat)
+        gj_solve(augmat, 2, 1, res)
+        # Then
+        np.testing.assert_array_almost_equal(res, expect)
+
     def test_general_matrix(self):
         # Test Gauss Jordan solve.
         """
@@ -118,7 +147,7 @@ def test_inverse(self):
         mat = [[1.0, 2.0, 2.5], [2.5, 1.0, 0.0], [0.0, 0.0, 1.0]]
         b = np.identity(3).ravel().tolist()
         A = np.zeros((3, 6)).ravel().tolist()
-        augmented_matrix(np.ravel(mat), b, 3, 3, A)
+        augmented_matrix(np.ravel(mat), b, 3, 3, 3, A)
         result = np.zeros((3, 3)).ravel().tolist()
 
         # When

pysph/sph/wc/crksph.py

@@ -102,7 +102,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,
                         grad_m2[d2*gam + d*alp + bet] += V * temp2
 
         identity(m2inv, d)
-        augmented_matrix(m2, m2inv, d, d, temp_aug_m2)
+        augmented_matrix(m2, m2inv, d, d, d, temp_aug_m2)
 
         # If is_singular > 0 then matrix was singular
         is_singular = gj_solve(temp_aug_m2, d, d, m2inv)

pysph/sph/wc/density_correction.py

@@ -122,7 +122,7 @@ def loop_all(self, d_idx, d_rho, d_x, d_y, s_x, s_y, d_h, s_h, s_m,
         res[0] = 1.0
         res[1] = 0.0
         res[2] = 0.0
-        augmented_matrix(amls, res, 3, 1, aug_mls)
+        augmented_matrix(amls, res, 3, 1, 3, aug_mls)
         gj_solve(aug_mls, n, 1, res)
         b0 = res[0]
         b1 = res[1]

pysph/sph/wc/linalg.py

@@ -62,26 +62,31 @@ def mat_vec_mult(a=[1.0, 0.0], b=[1.0, 0.0], n=3, result=[0.0, 0.0]):
         result[i] = s
 
 
-def augmented_matrix(A=[0.0, 0.0], b=[0.0, 0.0], n=3, na=1, result=[0.0, 0.0]):
+def augmented_matrix(A=[0.0, 0.0], b=[0.0, 0.0], n=3, na=1, nmax=3,
+                     result=[0.0, 0.0]):
     """Create augmented matrix.
 
-    Given flattened matrix, `A`, and flattened columns `b` with `n` rols and
-    `na` additional columns, put these in result. Result must be already
-    allocated and be flattened.
+    Given flattened matrix, `A` of max rows/columns `nmax`, and flattened
+    columns `b` with `n` rows of interest and `na` additional columns, put
+    these in `result`. Result must be already allocated and be flattened.
+    The `result` will contain `(n + na)*n` first entries as the
+    augmented_matrix.
+
 
     Parameters
     ----------
     A: list: given matrix.
     b: list: additional columns to be augmented.
-    n: int : number of rows/columns in `A`.
+    n: int : number of rows/columns to use from `A`.
     na: int: number of added columns in `b`.
-    result: list: must have size of (n + na)*n.
+    nmax: int: the maximum dimension 'A'
+    result: list: must have size of (nmax + na)*n.
     """
     i, j, nt = declare('int', 3)
     nt = n + na
     for i in range(n):
         for j in range(n):
-            result[nt*i + j] = A[n * i + j]
+            result[nt*i + j] = A[nmax * i + j]
         for j in range(na):
             result[nt*i + n + j] = b[na*i + j]
 
@@ -97,7 +102,7 @@ def gj_solve(m=[1., 0.], n=3, nb=1, result=[0.0, 0.0]):
     ----------
 
     m : list: a flattened list representing the augmented matrix [A|b].
-    n : int: number of columns/rows of A.
+    n : int: number of columns/rows used from A in augmented_matrix.
     nb: int: number of columns added to A.
     result: list: with size n*nb
 

pysph/tools/interpolator.py

@@ -8,9 +8,12 @@
 from pysph.base.utils import get_particle_array
 from pysph.base.kernels import Gaussian
 from pysph.base.nnps import LinkedListNNPS as NNPS
-from pysph.sph.equation import Equation
+from pysph.sph.equation import Equation, Group
 from pysph.sph.acceleration_eval import AccelerationEval
 from pysph.sph.sph_compiler import SPHCompiler
+from compyle.api import declare
+from pysph.sph.wc.linalg import gj_solve, augmented_matrix
+from pysph.sph.basic_equations import SummationDensity
 
 
 class InterpolateFunction(Equation):
@@ -35,6 +38,117 @@ def loop(self, d_idx, s_idx, s_rho, s_m, s_temp_prop, d_prop, WIJ):
         d_prop[d_idx] += s_m[s_idx]/s_rho[s_idx]*WIJ*s_temp_prop[s_idx]
 
 
+class SPHFirstOrderApproximationPreStep(Equation):
+    def __init__(self, dest, sources, dim=1):
+        self.dim = dim
+
+        super(SPHFirstOrderApproximationPreStep, self).__init__(dest, sources)
+
+    def initialize(self, d_idx, d_moment):
+        i, j = declare('int', 2)
+
+        for i in range(4):
+            for j in range(4):
+                d_moment[16*d_idx + j+4*i] = 0.0
+
+    def loop(self, d_idx, s_idx, d_h, s_h, s_x, s_y, s_z, d_x, d_y, d_z, s_rho,
+             s_m, WIJ, XIJ, DWIJ, d_moment):
+        Vj = s_m[s_idx] / s_rho[s_idx]
+        i16 = declare('int')
+        i16 = 16*d_idx
+
+        d_moment[i16+0] += WIJ * Vj
+
+        d_moment[i16+1] += -XIJ[0] * WIJ * Vj
+        d_moment[i16+2] += -XIJ[1] * WIJ * Vj
+        d_moment[i16+3] += -XIJ[2] * WIJ * Vj
+
+        d_moment[i16+4] += DWIJ[0] * Vj
+        d_moment[i16+8] += DWIJ[1] * Vj
+        d_moment[i16+12] += DWIJ[2] * Vj
+
+        d_moment[i16+5] += -XIJ[0] * DWIJ[0] * Vj
+        d_moment[i16+6] += -XIJ[1] * DWIJ[0] * Vj
+        d_moment[i16+7] += -XIJ[2] * DWIJ[0] * Vj
+
+        d_moment[i16+9] += - XIJ[0] * DWIJ[1] * Vj
+        d_moment[i16+10] += -XIJ[1] * DWIJ[1] * Vj
+        d_moment[i16+11] += -XIJ[2] * DWIJ[1] * Vj
+
+        d_moment[i16+13] += -XIJ[0] * DWIJ[2] * Vj
+        d_moment[i16+14] += -XIJ[1] * DWIJ[2] * Vj
+        d_moment[i16+15] += -XIJ[2] * DWIJ[2] * Vj
+
+
+class SPHFirstOrderApproximation(Equation):
+    """ First order SPH approximation
+        -------------
+        The method used to solve the linear system in this function is not same
+        as in the reference. In the function `Ax=b` is solved where `A :=
+        moment` (Moment matrix) and `b := p_sph` (Property calculated using
+        basic SPH).
+        The calculation need the `moment` to be evaluated before this step
+        which is done in `SPHFirstOrderApproximationPreStep`
+
+        References
+        ----------
+
+        .. [Liu2006] M.B. Liu, G.R. Liu, "Restoring particle consistency in
+           smoothed particle hydrodynamics", Applied Numerical Mathematics
+           Volume 56, Issue 1 2006, Pages 19-36, ISSN 0168-9274
+    """
+    def _get_helpers_(self):
+        return [gj_solve, augmented_matrix]
+
+    def __init__(self, dest, sources, dim=1):
+        self.dim = dim
+
+        super(SPHFirstOrderApproximation, self).__init__(dest, sources)
+
+    def initialize(self, d_idx, d_prop, d_p_sph):
+        i = declare('int')
+
+        for i in range(3):
+            d_prop[4*d_idx+i] = 0.0
+            d_p_sph[4*d_idx+i] = 0.0
+
+    def loop(self, d_idx, d_h, s_h, s_x, s_y, s_z, d_x, d_y, d_z, s_rho,
+             s_m, WIJ, DWIJ, s_temp_prop, d_p_sph, s_idx):
+        i4 = declare('int')
+        Vj = s_m[s_idx] / s_rho[s_idx]
+        pj = s_temp_prop[s_idx]
+        i4 = 4*d_idx
+
+        d_p_sph[i4+0] += pj * WIJ * Vj
+        d_p_sph[i4+1] += pj * DWIJ[0] * Vj
+        d_p_sph[i4+2] += pj * DWIJ[1] * Vj
+        d_p_sph[i4+3] += pj * DWIJ[2] * Vj
+
+    def post_loop(self, d_idx, d_moment, d_prop, d_p_sph):
+
+        a_mat = declare('matrix(16)')
+        aug_mat = declare('matrix(20)')
+        b = declare('matrix(4)')
+        res = declare('matrix(4)')
+
+        i, n, i16, i4 = declare('int', 4)
+        i16 = 16*d_idx
+        i4 = 4*d_idx
+        for i in range(16):
+            a_mat[i] = d_moment[i16+i]
+        for i in range(20):
+            aug_mat[i] = 0.0
+        for i in range(4):
+            b[i] = d_p_sph[4*d_idx+i]
+            res[i] = 0.0
+
+        n = self.dim + 1
+        augmented_matrix(a_mat, b, n, 1, 4, aug_mat)
+        gj_solve(aug_mat, n, 1, res)
+        for i in range(4):
+            d_prop[i4+i] = res[i]
+
+
 def get_bounding_box(particle_arrays, tight=False, stretch=0.05):
     """Find the size of the domain given a sequence of particle arrays.
 
@@ -93,7 +207,7 @@ class Interpolator(object):
 
     def __init__(self, particle_arrays, num_points=125000, kernel=None,
                  x=None, y=None, z=None, domain_manager=None,
-                 equations=None, use_shepard=True):
+                 equations=None, method='shepard'):
         """
         The x, y, z coordinates need not be specified, and if they are not,
         the bounds of the interpolated domain is automatically computed and
@@ -119,9 +233,9 @@ def __init__(self, particle_arrays, num_points=125000, kernel=None,
         equations: sequence
             A sequence of equations or groups.  Defaults to None.  This is
             used only if the default interpolation equations are inadequate.
-        use_shepard: bool
-            Use Shepard interpolation for the interpolation when no equations
-            are specified. If False, a simple SPH interpolation is performed.
+        method : str
+            String with the following allowed methods: 'shepard', 'sph',
+            'order1'
         """
         self._set_particle_arrays(particle_arrays)
         bounds = get_bounding_box(self.particle_arrays)
@@ -138,7 +252,9 @@ def __init__(self, particle_arrays, num_points=125000, kernel=None,
         self.equations = equations
         self.func_eval = None
         self.domain_manager = domain_manager
-        self.use_shepard = use_shepard
+        self.method = method
+        if method not in ['sph', 'shepard', 'order1']:
+            raise RuntimeError('%s method is not implemented' % (method))
         if x is None and y is None and z is None:
             self.set_domain(bounds, shape)
         else:
@@ -200,7 +316,7 @@ def set_domain(self, bounds, shape):
         x, y, z = self._create_default_points(self.bounds, self.shape)
         self.set_interpolation_points(x, y, z)
 
-    def interpolate(self, prop, gradient=False):
+    def interpolate(self, prop, comp=0):
         """Interpolate given property.
 
         Parameters
@@ -209,19 +325,29 @@ def interpolate(self, prop, gradient=False):
         prop: str
             The name of the property to interpolate.
 
-        gradient: bool
-            Evaluate gradient and not function.
+        comp: int
+            The component of the gradient required
 
         Returns
         -------
         A numpy array suitably shaped with the property interpolated.
         """
+        assert isinstance(comp, int), 'Error: only interger value is allowed'
         for array in self.particle_arrays:
             data = array.get(prop, only_real_particles=False)
             array.get('temp_prop', only_real_particles=False)[:] = data
 
         self.func_eval.compute(0.0, 0.1)  # These are junk arguments.
-        result = self.pa.prop.copy()
+        if comp and (self.method in ['sph', 'shepard']):
+            raise RuntimeError("Error: use 'order1' method to evaluate"
+                               "gradient")
+        elif self.method in ['sph', 'shepard']:
+            result = self.pa.prop.copy()
+        else:
+            if comp > 3:
+                raise RuntimeError("Error: Only 0, 1, 2, 3 allowed")
+            result = self.pa.prop[comp::4].copy()
+
         result.shape = self.shape
         return result.squeeze()
 
@@ -279,27 +405,50 @@ def _create_particle_array(self, x, y, z):
 
         hmax = self._get_max_h_in_arrays()
         h = hmax*np.ones_like(xr)
-        prop = np.zeros_like(xr)
         pa = get_particle_array(
             name='interpolate',
             x=xr, y=yr, z=zr, h=h,
-            number_density=np.zeros_like(xr),
-            prop=prop,
-            grad_x=np.zeros_like(xr),
-            grad_y=np.zeros_like(xr),
-            grad_z=np.zeros_like(xr)
+            number_density=np.zeros_like(xr)
         )
+        if self.method in ['sph', 'shepard']:
+            pa.add_property('prop')
+        else:
+            pa.add_property('moment', stride=16)
+            pa.add_property('p_sph', stride=4)
+            pa.add_property('prop', stride=4)
+
         return pa
 
     def _compile_acceleration_eval(self, arrays):
         names = [x.name for x in self.particle_arrays]
         if self.equations is None:
-            if self.use_shepard:
-                equations = [InterpolateFunction(dest='interpolate',
-                                                 sources=names)]
+            if self.method == 'shepard':
+                equations = [
+                    InterpolateFunction(dest='interpolate',
+                                        sources=names)
+                ]
+            elif self.method == 'sph':
+                equations = [
+                    InterpolateSPH(dest='interpolate',
+                                        sources=names)
+                ]
             else:
-                equations = [InterpolateSPH(dest='interpolate',
-                                            sources=names)]
+                equations = [
+                    Group(equations=[
+                        SummationDensity(dest=name, sources=names)
+                        for name in names],
+                          real=False),
+                    Group(equations=[
+                        SPHFirstOrderApproximationPreStep(dest='interpolate',
+                                                          sources=names,
+                                                          dim=self.dim)],
+                          real=True),
+                    Group(equations=[
+                        SPHFirstOrderApproximation(dest='interpolate',
+                                                   sources=names,
+                                                   dim=self.dim)],
+                          real=True)
+                ]
         else:
             equations = self.equations
         self.func_eval = AccelerationEval(arrays, equations, self.kernel)