From 8e815ad99ba89d4ea7a5538b11ffd63f9de4a131 Mon Sep 17 00:00:00 2001
From: Jean-Matthieu Etancelin <jean-matthieu.etancelin@univ-reims.fr>
Date: Tue, 18 Nov 2014 17:52:28 +0100
Subject: [PATCH] Improve GPU-diffusion profiling. Add a full update ghost. fix
 advection_dir topology parameters. fix tools parameters. Add update ghosts
 test.

---
 HySoP/hysop/gpu/gpu_diffusion.py              |  14 +-
 .../numerics/tests/test_update_ghosts.py      | 170 ++++++++++++++++++
 HySoP/hysop/numerics/update_ghosts.py         |  25 ++-
 HySoP/hysop/operator/advection_dir.py         |  13 +-
 HySoP/hysop/tools/parameters.py               |   2 +-
 5 files changed, 199 insertions(+), 25 deletions(-)
 create mode 100644 HySoP/hysop/numerics/tests/test_update_ghosts.py

diff --git a/HySoP/hysop/gpu/gpu_diffusion.py b/HySoP/hysop/gpu/gpu_diffusion.py
index 4131d9d28..e70f3ab4a 100644
--- a/HySoP/hysop/gpu/gpu_diffusion.py
+++ b/HySoP/hysop/gpu/gpu_diffusion.py
@@ -158,6 +158,8 @@ class GPUDiffusion(DiscreteOperator, GPUOperator):
                            nb_part_per_wi, tile_size)
         self.num_diffusion = KernelLauncher(
             prg.diffusion, self.cl_env.queue, gwi, lwi)
+        self.copy = KernelLauncher(cl.enqueue_copy,
+                                   self.cl_env.queue)
 
     def _compute_diffusion(self, simulation):
         assert self.field_tmp is not None
@@ -168,8 +170,10 @@ class GPUDiffusion(DiscreteOperator, GPUOperator):
             self.gpu_precision(self.viscosity * simulation.timeStep),
             self._mesh_size,
             wait_for=wait_evt)
-        c_evt = cl.enqueue_copy(self.cl_env.queue, self.field.gpu_data[0],
-                                self.field_tmp, wait_for=[d_evt])
+        c_evt = self.copy.launch_sizes_in_args(
+            self.field.gpu_data[0], self.field_tmp, wait_for=[d_evt])
+        #c_evt = cl.enqueue_copy(self.cl_env.queue, self.field.gpu_data[0],
+        #                        self.field_tmp, wait_for=[d_evt])
         self.field.events.append(c_evt)
 
     def set_field_tmp(self, field_tmp):
@@ -277,3 +281,9 @@ class GPUDiffusion(DiscreteOperator, GPUOperator):
 
     def apply(self, simulation):
         self._compute(simulation)
+
+    def get_profiling_info(self):
+        for k in [self.num_diffusion, self.copy]:
+            if k is not None:
+                for p in k.profile:
+                    self.profiler += p
diff --git a/HySoP/hysop/numerics/tests/test_update_ghosts.py b/HySoP/hysop/numerics/tests/test_update_ghosts.py
new file mode 100644
index 000000000..d050f43ff
--- /dev/null
+++ b/HySoP/hysop/numerics/tests/test_update_ghosts.py
@@ -0,0 +1,170 @@
+import parmepy.tools.numpywrappers as npw
+from parmepy.numerics.update_ghosts import UpdateGhosts, UpdateGhostsFull
+from parmepy.domain.box import Box
+from parmepy.fields.continuous import Field
+from parmepy.tools.parameters import Discretization
+import numpy as np
+
+
+def _setup(ghosts, topo_dim):
+    dom = Box()
+
+    f = Field(dom, isVector=True)
+    resolTopo = Discretization([33, 33, 33], ghosts=ghosts)
+    topo = dom.create_topology(resolTopo, dim=topo_dim)
+    f.discretize(topo)
+    df = f.discreteFields[topo]
+    for i, dfd in enumerate(df.data):
+        dfd[...] = 0.
+        dfd[topo.mesh.iCompute] = 1. * (i + 1)
+    return df, topo
+
+
+def verify(df, gh):
+    sli = (
+        (slice(0, gh[0]), slice(gh[1], -gh[1]), slice(gh[2], -gh[2])),
+        (slice(-gh[0], None), slice(gh[1], -gh[1]), slice(gh[2], -gh[2])),
+        (slice(gh[0], -gh[0]), slice(0, gh[1]), slice(gh[2], -gh[2])),
+        (slice(gh[0], -gh[0]), slice(-gh[1], None), slice(gh[2], -gh[2])),
+        (slice(gh[0], -gh[0]), slice(gh[1], -gh[1]), slice(0, gh[2])),
+        (slice(gh[0], -gh[0]), slice(gh[1], -gh[1]), slice(-gh[2], None)))
+    slni = (
+        (slice(0, gh[0]), slice(0, gh[1]), slice(None)),
+        (slice(0, gh[0]), slice(None), slice(0, gh[2])),
+        (slice(None), slice(0, gh[1]), slice(0, gh[2])),
+        (slice(-gh[0], None), slice(-gh[1], None), slice(None)),
+        (slice(-gh[0], None), slice(None), slice(-gh[2], None)),
+        (slice(None), slice(-gh[1], None), slice(-gh[2], None)))
+    for i, dfd in enumerate(df.data):
+        for s in sli:
+            assert np.max(dfd[s]) == np.min(dfd[s]) and \
+                np.max(dfd[s]) == 1. * (i + 1)
+        for s in slni:
+            assert np.max(dfd[s]) == np.min(dfd[s]) and \
+                np.max(dfd[s]) == 0.
+
+
+def verify_full(df, gh):
+    for i, dfd in enumerate(df.data):
+        for s in ((slice(0, gh[0]), slice(None), slice(None)),
+                  (slice(-gh[0], None), slice(None), slice(None)),
+                  (slice(None), slice(0, gh[1]), slice(None)),
+                  (slice(None), slice(-gh[1], None), slice(None)),
+                  (slice(None), slice(None), slice(0, gh[2])),
+                  (slice(None), slice(None), slice(-gh[2], None))):
+            assert np.max(dfd[s]) == np.min(dfd[s]) and \
+                np.max(dfd[s]) == 1. * (i + 1)
+        assert np.max(dfd) == np.min(dfd) and \
+            np.max(dfd) == 1. * (i + 1)
+
+
+def test_update_ghosts_simple_1D():
+    gh = [1, 1, 1]
+    df, topo, = _setup(gh, 1)
+    update = UpdateGhosts(topo, 3)
+    update(df.data)
+    verify(df, gh)
+
+
+def test_update_ghosts_simple_2D():
+    gh = [1, 1, 1]
+    df, topo, = _setup(gh, 2)
+    update = UpdateGhosts(topo, 3)
+    update(df.data)
+    verify(df, gh)
+
+
+def test_update_ghosts_simple_3D():
+    gh = [1, 1, 1]
+    df, topo, = _setup(gh, 3)
+    update = UpdateGhosts(topo, 3)
+    update(df.data)
+    verify(df, gh)
+
+
+def test_update_ghosts_1D():
+    gh = [2, 3, 4]
+    df, topo, = _setup(gh, 1)
+    update = UpdateGhosts(topo, 3)
+    update(df.data)
+    verify(df, gh)
+
+
+def test_update_ghosts_2D():
+    gh = [2, 3, 4]
+    df, topo, = _setup(gh, 2)
+    update = UpdateGhosts(topo, 3)
+    update(df.data)
+    verify(df, gh)
+
+
+def test_update_ghosts_3D():
+    gh = [2, 3, 4]
+    df, topo, = _setup(gh, 3)
+    update = UpdateGhosts(topo, 3)
+    update(df.data)
+    verify(df, gh)
+
+
+def test_update_ghosts_full_simple_1D():
+    gh = [1, 1, 1]
+    df, topo, = _setup(gh, 1)
+    update = UpdateGhostsFull(topo, 3)
+    update(df.data)
+    verify_full(df, gh)
+
+
+def test_update_ghosts_full_simple_2D():
+    gh = [1, 1, 1]
+    df, topo, = _setup(gh, 2)
+    update = UpdateGhostsFull(topo, 3)
+    update(df.data)
+    verify_full(df, gh)
+
+
+def test_update_ghosts_full_simple_3D():
+    gh = [1, 1, 1]
+    df, topo, = _setup(gh, 3)
+    update = UpdateGhostsFull(topo, 3)
+    update(df.data)
+    verify_full(df, gh)
+
+
+def test_update_ghosts_full_1D():
+    gh = [2, 3, 4]
+    df, topo, = _setup(gh, 1)
+    update = UpdateGhostsFull(topo, 3)
+    update(df.data)
+    verify_full(df, gh)
+
+
+def test_update_ghosts_full_2D():
+    gh = [2, 3, 4]
+    df, topo, = _setup(gh, 2)
+    update = UpdateGhostsFull(topo, 3)
+    update(df.data)
+    verify_full(df, gh)
+
+
+def test_update_ghosts_full_3D():
+    gh = [2, 3, 4]
+    df, topo, = _setup(gh, 3)
+    update = UpdateGhostsFull(topo, 3)
+    update(df.data)
+    verify_full(df, gh)
+
+
+if __name__ == '__main__':
+    for f in (test_update_ghosts_simple_1D,
+              test_update_ghosts_simple_2D,
+              test_update_ghosts_simple_3D,
+              test_update_ghosts_1D,
+              test_update_ghosts_2D,
+              test_update_ghosts_3D,
+              test_update_ghosts_full_simple_1D,
+              test_update_ghosts_full_simple_2D,
+              test_update_ghosts_full_simple_3D,
+              test_update_ghosts_full_1D,
+              test_update_ghosts_full_2D,
+              test_update_ghosts_full_3D):
+        f()
diff --git a/HySoP/hysop/numerics/update_ghosts.py b/HySoP/hysop/numerics/update_ghosts.py
index f18dce221..8e9c9a997 100644
--- a/HySoP/hysop/numerics/update_ghosts.py
+++ b/HySoP/hysop/numerics/update_ghosts.py
@@ -46,7 +46,7 @@ class UpdateGhosts(object):
         self._recvbuffer = []
         # domain dimension
         self._dim = self.topology.domain.dimension
-        
+
         self._setup_slices()
 
         ## shape of numpy arrays to be updated.
@@ -123,7 +123,7 @@ class UpdateGhosts(object):
                 if self.topology.shape[d] == 1]
         for d in dirs:
             self._applyBC_in_dir(variables, d)
-                
+
     def _applyBC_in_dir(self, variables, d):
         """Apply periodic boundary condition in direction d."""
         for v in variables:
@@ -151,7 +151,7 @@ class UpdateGhosts(object):
         self.applyBC(variables)
 
     def _apply_in_dir(self, variables, d, i):
-        """Communicate ghosts values in direction d for neighbour 
+        """Communicate ghosts values in direction d for neighbour
         in direction i of the topology"""
         comm = self.topology.comm
         rank = self.topology.rank
@@ -219,7 +219,7 @@ class UpdateGhostsFull(UpdateGhosts):
         at each call.
         nbElements and memshape will be used to allocate memory for local
         buffers used for send-recv process.
-        This version differs from UpdateGhosts by computing also ghosts values 
+        This version differs from UpdateGhosts by computing also ghosts values
         in edges and corners of the domain. The directions are computed in reversed order.
         """
         super(UpdateGhostsFull, self).__init__(topo, nbElements)
@@ -244,13 +244,14 @@ class UpdateGhostsFull(UpdateGhosts):
                 self._g_tonext.append(list(defslice))
                 self._g_tonext[d][d] = slice(-2 * self.ghosts[d],
                                              -self.ghosts[d])
-                
-                ## Slices for other directions corresmonding to directions not
-                ## yet exchanged : x < d. For directions x > d, the slices is a 
-                ## full slice that includes ghosts. This assumes that directions 
-                ## x > d have already been computed (by communications or local
+
+                ## Slices for other directions corresponding to directions not
+                ## yet exchanged : x > d. For directions x < d, the slices is a
+                ## full slice that includes ghosts. This assumes that directions
+                ## x < d have already been computed (by communications or local
                 ## exchanges)
-                otherDim = [x for x in xrange(self._dim) if x < d]
+                #                if d == 0:
+                otherDim = [x for x in xrange(self._dim) if x > d]
                 for d2 in otherDim:
                     self._g_fromprevious[d][d2] = slice(self.ghosts[d2],
                                                         -self.ghosts[d2])
@@ -285,7 +286,7 @@ class UpdateGhostsFull(UpdateGhosts):
         assert len(exchange_dir) + len(localBC_dir) == self._dim
 
         i = 0
-        for d in xrange(self._dim-1, -1, -1):
+        for d in xrange(self._dim):
             if d in localBC_dir:
                 self._applyBC_in_dir(variables, d)
             elif d in exchange_dir:
@@ -293,5 +294,3 @@ class UpdateGhostsFull(UpdateGhosts):
                 # update index in neighbours list
                 i += 1
                 # End of loop through send/recv directions.
-
-
diff --git a/HySoP/hysop/operator/advection_dir.py b/HySoP/hysop/operator/advection_dir.py
index d08deff01..75c004af2 100644
--- a/HySoP/hysop/operator/advection_dir.py
+++ b/HySoP/hysop/operator/advection_dir.py
@@ -30,7 +30,7 @@ class AdvectionDir(Computational):
 
     @debug
     def __init__(self, velocity, direction, advected_fields=None,
-                 name_suffix='', cutdir=None, **kwds):
+                 name_suffix='', **kwds):
         ## Get the other arguments to pass to the discrete operators
         self._kwds = kwds.copy()
         self._kwds.pop('discretization')
@@ -78,8 +78,6 @@ class AdvectionDir(Computational):
         ## Positions of the particles
         self.particle_positions = None
 
-        self._default_cutdir = cutdir
-
     @debug
     def discretize(self):
         if self._is_discretized:
@@ -93,12 +91,9 @@ class AdvectionDir(Computational):
                 raise ValueError("Multiscale advection is not yet supported "
                                  "in pure Python, use Scales or GPU.")
 
-        ## Topology cutdir for parallel advection
-        if self._default_cutdir is None:
-            cutdir = [False] * self.domain.dimension
-            cutdir[-1] = True
-        else:
-            cutdir = self._default_cutdir
+        ## Default topology cutdir for parallel advection
+        cutdir = [False] * self.domain.dimension
+        cutdir[-1] = True
 
         if self._single_topo:
             # One topo for all fields ...
diff --git a/HySoP/hysop/tools/parameters.py b/HySoP/hysop/tools/parameters.py
index 53db939b5..b24a6f9cf 100644
--- a/HySoP/hysop/tools/parameters.py
+++ b/HySoP/hysop/tools/parameters.py
@@ -61,7 +61,7 @@ class Discretization(namedtuple("Discretization", ['resolution', 'ghosts'])):
             assert ghosts.size == resolution.size, msg
             assert all(ghosts >= 0)
         else:
-            ghosts = npw.dim_zeros(resolution.size)
+            ghosts = npw.int_zeros(resolution.size)
         return super(Discretization, cls).__new__(cls, resolution, ghosts)
 
     def __eq__(self, other):
-- 
GitLab