From 6e5ad90db812c99e0fd10c214fcb523670519a2b Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Thu, 14 Jun 2018 15:53:44 +0200
Subject: [PATCH] Second attempt, added a parameter to disable diagonal ghost
 exchange when not desired (ghost accumulation)

---
 hysop/constants.py                       |   2 +
 hysop/constants.py.in                    |   2 +
 hysop/fields/cartesian_discrete_field.py |  50 +++--
 hysop/fields/discrete_field.py           |   6 +-
 hysop/fields/ghost_exchangers.py         |  62 ++++---
 hysop/fields/tests/test_cartesian.py     | 225 ++++++++++++-----------
 hysop/mesh/cartesian_mesh.py             |  72 ++++----
 7 files changed, 245 insertions(+), 174 deletions(-)

diff --git a/hysop/constants.py b/hysop/constants.py
index e93a47331..1ed7b2c72 100644
--- a/hysop/constants.py
+++ b/hysop/constants.py
@@ -102,6 +102,8 @@ def implementation_to_backend(implementation):
 
 GhostOperation = EnumFactory.create('GhostOperation',
         ['EXCHANGE', 'ACCUMULATE'])
+GhostMask = EnumFactory.create('GhostMask',
+        ['CROSS', 'FULL'])
 ExchangeMethod = EnumFactory.create('ExchangeMethod',
         ['ISEND_IRECV', 'NEIGHBOR_ALL_TO_ALL_V', 'NEIGHBOR_ALL_TO_ALL_W'])
 
diff --git a/hysop/constants.py.in b/hysop/constants.py.in
index ef2f9bc7b..d12240e72 100644
--- a/hysop/constants.py.in
+++ b/hysop/constants.py.in
@@ -102,6 +102,8 @@ def implementation_to_backend(implementation):
 
 GhostOperation = EnumFactory.create('GhostOperation',
         ['EXCHANGE', 'ACCUMULATE'])
+GhostMask = EnumFactory.create('GhostMask',
+        ['CROSS', 'FULL'])
 ExchangeMethod = EnumFactory.create('ExchangeMethod',
         ['ISEND_IRECV', 'NEIGHBOR_ALL_TO_ALL_V', 'NEIGHBOR_ALL_TO_ALL_W'])
 
diff --git a/hysop/fields/cartesian_discrete_field.py b/hysop/fields/cartesian_discrete_field.py
index 999e595d0..21304c2a6 100644
--- a/hysop/fields/cartesian_discrete_field.py
+++ b/hysop/fields/cartesian_discrete_field.py
@@ -9,7 +9,7 @@ Documentation and examples can be found in :ref:`fields`.
 from hysop import vprint, dprint, MPI
 from hysop.deps import np, hashlib
 from hysop.core.arrays.all import HostArray, OpenClArray
-from hysop.constants import Backend, DirectionLabels, GhostOperation, ExchangeMethod
+from hysop.constants import Backend, DirectionLabels, GhostOperation, GhostMask, ExchangeMethod
 from hysop.topology.topology import Topology, TopologyView
 from hysop.tools.decorators import debug, static_vars
 from hysop.tools.types import check_instance, to_tuple, first_not_None, to_list
@@ -437,8 +437,8 @@ CartesianDiscreteFieldView (id={}, tag={})
         for d in xrange(self.nb_components):
             self.backend.rand(out=self.data[d], **kwds)
    
-    def initialize(self, formula, vectorize=False, 
-            only_finite=True, exchange_ghosts=True, 
+    def initialize(self, formula, vectorize=False, only_finite=True, 
+            exchange_ghosts=True, exchange_kwds=None,
             quiet=False, reorder=None, **kwds):
         """
         Initialize the field components
@@ -455,7 +455,12 @@ CartesianDiscreteFieldView (id={}, tag={})
         only_finite: bool, optional
             Check that initialized values are not +inf, -inf or NaN.
             Defaults to True.
-        reorder: str or tuple of str
+        exchange_ghosts: bool, optional, defaults to True
+            Should we exchange ghosts after initialization ?
+        exchange_kwds: dict, optional,
+            Extra exchange keyword arguments passed to ghost exchange.
+            Only used if exchange_ghosts is set to True.
+        reorder: str or tuple of str, optional
             Reorder all kwd according to current topology state.
             ie. kwd should be contained in kwds and kwds[kwd] should
                 be of length self.domain.dim.
@@ -576,7 +581,8 @@ CartesianDiscreteFieldView (id={}, tag={})
                 vprint(msg)
         
         if exchange_ghosts:
-            evt = self.exchange_ghosts()
+            exchange_kwds = first_not_None(exchange_kwds, {})
+            evt = self.exchange_ghosts(**exchange_kwds)
             if (evt is not None):
                 evt.wait()
 
@@ -710,26 +716,41 @@ CartesianDiscreteFieldView (id={}, tag={})
         with printoptions(**_print_opts):
             print strarr
     
+    def accumulate_ghosts(self, **kwds):
+        """
+        Exchange ghosts using cached ghost exchangers which are built at first use.
+        ie. Exchange every ghosts components of self.data using current topology state.
+        Specialization for ghost summation.
+
+        Defaults to ghost accumulation excluding diagonals.
+        """
+        return self.exchange_ghosts(ghost_op=GhostOperation.ACCUMULATE, 
+                                    ghost_mask=GhostMask.CROSS,
+                                    **kwds)
+    
     def exchange_ghosts(self, components=None, directions=None,
-            ghosts=None, ghost_op=None, exchange_method=None,
+            ghosts=None, ghost_op=None, ghost_mask=None, exchange_method=None, 
             evt=None, build_launcher=False, **kwds):
         """
         Exchange ghosts using cached ghost exchangers which are built at first use.
         ie. Exchange every ghosts components of self.data using current topology state.
-
-        Defaults to ghosts exchange (ie. overwrite operation).
+        
+        Defaults to full ghosts exchange, including diagonals (ie. overwrite operation).
         """
         ghosts          = to_tuple(first_not_None(ghosts, self.ghosts))
         components      = to_tuple(first_not_None(components, range(len(self.data))))
         directions      = to_tuple(first_not_None(directions, range(self.dim)))
         ghost_op        = first_not_None(ghost_op, GhostOperation.EXCHANGE)
+        ghost_mask      = first_not_None(ghost_mask, GhostMask.FULL)
         exchange_method = first_not_None(exchange_method, ExchangeMethod.ISEND_IRECV)
         if self.has_ghosts() and any(mg>0 for mg in ghosts):
             topology_state = self.topology_state
-            key = (topology_state, ghosts, ghost_op, exchange_method, components, directions)
+            key = (topology_state, ghosts, ghost_op, ghost_mask, 
+                            exchange_method, components, directions)
             if (key not in self._dfield._ghost_exchangers):
                 self._dfield._ghost_exchangers[key] = self.build_ghost_exchanger(ghosts=ghosts, 
-                        ghost_op=ghost_op, components=components, directions=directions,
+                        ghost_op=ghost_op, ghost_mask=ghost_mask,
+                        components=components, directions=directions,
                         exchange_method=exchange_method)
             ghost_exchanger = self._dfield._ghost_exchangers[key]
             if build_launcher:
@@ -746,9 +767,10 @@ CartesianDiscreteFieldView (id={}, tag={})
         return first_not_None(new_evt, evt)
 
     def build_ghost_exchanger(self, name=None, components=None, directions=None, 
-            data=None, ghosts=None, ghost_op=None, exchange_method=None):
+            data=None, ghosts=None, ghost_op=None, ghost_mask=None, exchange_method=None):
         """Build a ghost exchanger for cartesian discrete fields, possibly on different data."""
-        ghost_op   = first_not_None(ghost_op, GhostOperation.EXCHANGE)
+        ghost_op        = first_not_None(ghost_op, GhostOperation.EXCHANGE)
+        ghost_mask      = first_not_None(ghost_mask, GhostMask.FULL)
         exchange_method = first_not_None(exchange_method, ExchangeMethod.ISEND_IRECV)
         check_instance(ghost_op, GhostOperation)
         check_instance(exchange_method, ExchangeMethod)
@@ -757,7 +779,7 @@ CartesianDiscreteFieldView (id={}, tag={})
         components = to_tuple(first_not_None(components, range(len(data))))
         directions = to_tuple(first_not_None(directions, range(self.dim)))
         data       = tuple(data[i] for i in components)
-        ghosts = first_not_None(ghosts, self.ghosts)
+        ghosts     = first_not_None(ghosts, self.ghosts)
 
         d0 = data[0]
         if isinstance(data[0], (np.ndarray, HostArray)):
@@ -771,7 +793,7 @@ CartesianDiscreteFieldView (id={}, tag={})
         return CartesianDiscreteFieldGhostExchanger(name=name,
                 topology=self.topology, data=data, kind=kind, 
                 ghosts=ghosts, directions=directions,
-                ghost_op=ghost_op, exchange_method=exchange_method)
+                ghost_op=ghost_op, ghost_mask=ghost_mask, exchange_method=exchange_method)
 
 
     size               = property(_get_size)
diff --git a/hysop/fields/discrete_field.py b/hysop/fields/discrete_field.py
index d3392f838..fdc6b35c4 100644
--- a/hysop/fields/discrete_field.py
+++ b/hysop/fields/discrete_field.py
@@ -128,9 +128,11 @@ class DiscreteFieldView(TaggedObjectView, VariableTag):
         ie. Exchange every ghosts components of self.data using current topology state.
         Specialization for ghost summation.
         """
-        return self.exchange_ghosts(ghost_op=GhostOperation.ACCUMULATE, **kwds)
+        return self.exchange_ghosts(ghost_op=GhostOperation.ACCUMULATE, 
+                                    **kwds)
     
-    def exchange_ghosts(self, ghosts=None, ghost_op=GhostOperation.EXCHANGE, 
+    def exchange_ghosts(self, ghosts=None, 
+            ghost_op=GhostOperation.EXCHANGE, 
             evt=None, **kwds):
         """
         Exchange ghosts using cached ghost exchangers which are built at first use.
diff --git a/hysop/fields/ghost_exchangers.py b/hysop/fields/ghost_exchangers.py
index e71f4e506..3b8cbd793 100644
--- a/hysop/fields/ghost_exchangers.py
+++ b/hysop/fields/ghost_exchangers.py
@@ -2,7 +2,8 @@ from hysop.deps import np, hashlib
 from hysop.tools.types import check_instance, to_tuple, first_not_None
 from hysop.tools.numerics import default_invalid_value
 from hysop.tools.mpi_utils import iter_mpi_requests, dtype_to_mpi_type
-from hysop.constants import GhostOperation, Backend, DirectionLabels, ExchangeMethod
+from hysop.constants import GhostOperation, GhostMask, \
+        Backend, DirectionLabels, ExchangeMethod
 from hysop.topology.topology import Topology, TopologyView
 from hysop.topology.cartesian_topology import CartesianTopologyView
 from hysop.core.arrays.all import HostArray, OpenClArray
@@ -12,11 +13,14 @@ from hysop.backend.device.opencl import cl, clArray
 from hysop.backend.device.opencl.opencl_kernel_launcher import HostLauncherI
 
 class GhostExchanger(object):
-    def __init__(self, name, topology, data, ghost_op):
+    def __init__(self, name, topology, data, 
+            exchange_method, ghost_op, ghost_mask):
         check_instance(name, str)
         check_instance(topology, TopologyView)
         check_instance(data, tuple, minsize=1)
+        check_instance(exchange_method, ExchangeMethod)
         check_instance(ghost_op, GhostOperation)
+        check_instance(ghost_mask, GhostMask)
 
         dtype = data[0].dtype
         for d in data[1:]:
@@ -33,7 +37,9 @@ class GhostExchanger(object):
         self.backend = data[0].backend if hasattr(data[0], 'backend') else topology.backend
         self.host_backend = self.backend.host_array_backend
         self.base_tag = int(hashlib.sha1(name).hexdigest(), 16) % (104729)
+        self.exchange_method = exchange_method
         self.ghost_op = ghost_op
+        self.ghost_mask = ghost_mask
         self.base_dtype = base_dtype
         self.base_mpi_type = base_mpi_type
         self.name = name
@@ -42,29 +48,39 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
 
     def __init__(self, name, topology, data,
             kind=None, directions=None,
-            ghosts=None, ghost_op=None,
+            ghosts=None, ghost_op=None, ghost_mask=None,
             exchange_method=None):
         """
-        /!\ Ghosts are not exchanged in the diagonals /!\
-            Just Cartesian Communicator neighbours are
-            considered here.
-
-                 X   P   X
-                     ^
-                 P < P > P
-                     v
-                 X   P   X
-
-        Diagonal ghosts are set to NAN to ensure they are not used.
+        By default, we exchange all ghosts, including diagonals.
+        This is done by setting ghost_mask to GhostMask.FULL.
+        
+        Just Cartesian Communicator neighbours are
+        considered here (there is no direct communication
+        between diagonal processes).
+
+             X   P   X
+                 ^
+             P < P > P
+                 v
+             X   P   X
+        
+        If ghost_mask is set to GhostMask.CROSS, diagonal ghosts 
+        are set to NAN to ensure they are not used.
         """
 
         ghost_op = first_not_None(ghost_op, GhostOperation.EXCHANGE)
+        check_instance(ghost_op, GhostOperation)
+
+        ghost_mask = first_not_None(ghost_mask, GhostMask.FULL)
+        check_instance(ghost_mask, GhostMask)
 
         exchange_method = first_not_None(exchange_method,
                                 ExchangeMethod.NEIGHBOR_ALL_TO_ALL_V)
-
+        check_instance(exchange_method, ExchangeMethod)
+        
         super(CartesianDiscreteFieldGhostExchanger, self).__init__(topology=topology,
-                data=data, name=name, ghost_op=ghost_op)
+                data=data, name=name, exchange_method=exchange_method, 
+                ghost_op=ghost_op, ghost_mask=ghost_mask)
 
         mesh = self.mesh
         dim  = mesh.dim
@@ -84,13 +100,14 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
             ghosts*=dim
 
         self.directions   = directions
-        self.outer_ghosts = mesh.get_local_outer_ghost_slices(ghosts=ghosts)
-        self.inner_ghosts = mesh.get_local_inner_ghost_slices(ghosts=ghosts)
+        self.outer_ghosts = mesh.get_local_outer_ghost_slices(ghosts=ghosts, 
+                                                              ghost_mask=ghost_mask)
+        self.inner_ghosts = mesh.get_local_inner_ghost_slices(ghosts=ghosts,
+                                                              ghost_mask=ghost_mask)
         self.all_inner_ghost_slices = mesh.get_all_local_inner_ghost_slices(ghosts=ghosts)
         self.all_outer_ghost_slices = mesh.get_all_local_outer_ghost_slices(ghosts=ghosts)
-        self.dim          = dim
-        self.ghosts       = ghosts
-        self.exchange_method = exchange_method
+        self.dim             = dim
+        self.ghosts          = ghosts
 
         kind = first_not_None(kind, topology.backend.kind)
         if (kind == Backend.HOST):
@@ -198,7 +215,8 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
             lp.w_send_requests = {}
             lp.w_recv_requests = {}
 
-            if (ghost_op is GhostOperation.EXCHANGE) and (dim>1) and not self.mesh.isExchInclDiag:
+            if (ghost_op is GhostOperation.EXCHANGE) and (dim>1) and \
+                    (self.ghost_mask is GhostMask.CROSS):
                 # set all outer diagonal ghosts to NAN
                 for directions in all_outer_ghost_slices[dim]:
                     if all(ghosts[d]>0 for d in directions):
diff --git a/hysop/fields/tests/test_cartesian.py b/hysop/fields/tests/test_cartesian.py
index 16e2d8ae8..0e5dfae3d 100644
--- a/hysop/fields/tests/test_cartesian.py
+++ b/hysop/fields/tests/test_cartesian.py
@@ -1,7 +1,7 @@
 import os, subprocess, sys, time
 from hysop import __ENABLE_LONG_TESTS__
 from hysop.deps import it, np
-from hysop.constants import Backend, ExchangeMethod, GhostOperation, DirectionLabels
+from hysop.constants import Backend, ExchangeMethod, GhostOperation, GhostMask, DirectionLabels
 from hysop.tools.parameters import Discretization
 from hysop.tools.numerics import is_integer, is_fp
 from hysop.tools.numpywrappers import npw
@@ -85,45 +85,50 @@ def test_serial_initialization_2d():
         print '  Topo {}::{} | {}'.format(topo.full_tag, topo.backend.kind, topo.backend)
         dF0 = F0.discretize(topo)
         dF1 = F1.discretize(topo)
-        dF0.initialize(__random_init)
-        dF1.initialize(__random_init)
 
-        data = dF0.data + dF1.data
-        for (i,d) in enumerate(data):
-            print '    *buffer {}'.format(i)
-            buf = d.get().handle
-            slo_x = (slice(0,nghosts[0]),
-                     slice(nghosts[0],nghosts[0]+npts[0]-1),
-                     slice(nghosts[0]+npts[0]-1,None))
-            slo_y = (slice(0,nghosts[1]),
-                     slice(nghosts[1],nghosts[1]+npts[1]-1),
-                     slice(nghosts[1]+npts[1]-1,None))
-            sli_x = (slice(nghosts[0], 2*nghosts[0]),
-                     slice(2*nghosts[0],npts[0]-1),
-                     slice(npts[0]-1, nghosts[0]+npts[0]-1))
-            sli_y = (slice(nghosts[1], 2*nghosts[1]),
-                     slice(2*nghosts[1],npts[1]-1),
-                     slice(npts[1]-1, nghosts[1]+npts[1]-1))
-            assert buf.shape==(npts[0]+2*nghosts[0]-1,npts[1]+2*nghosts[1]-1)
-            assert buf[slo_x[0],slo_y[0]].shape == nghosts
-            assert buf[slo_x[2],slo_y[2]].shape == nghosts
-            assert buf[slo_x[1],slo_y[1]].shape == tuple([_-1 for _ in npts])
-            assert buf[sli_x[0],sli_y[0]].shape == nghosts
-            assert buf[sli_x[2],sli_y[2]].shape == nghosts
-            assert np.all(buf[slo_x[0],slo_y[1]] == buf[sli_x[2],slo_y[1]])
-            assert np.all(buf[slo_x[2],slo_y[1]] == buf[sli_x[0],slo_y[1]])
-            assert np.all(buf[slo_x[1],slo_y[0]] == buf[slo_x[1],sli_y[2]])
-            assert np.all(buf[slo_x[1],slo_y[2]] == buf[slo_x[1],sli_y[0]])
-            if topo.mesh.isExchInclDiag:
-                assert np.all(buf[slo_x[0],slo_y[0]]==buf[sli_x[2],sli_y[2]])
-                assert np.all(buf[slo_x[2],slo_y[0]]==buf[sli_x[0],sli_y[2]])
-                assert np.all(buf[slo_x[2],slo_y[2]]==buf[sli_x[0],sli_y[0]])
-                assert np.all(buf[slo_x[0],slo_y[2]]==buf[sli_x[2],sli_y[0]])
-            else:
-                assert np.all(np.isnan(buf[slo_x[0],slo_y[0]]))
-                assert np.all(np.isnan(buf[slo_x[2],slo_y[0]]))
-                assert np.all(np.isnan(buf[slo_x[2],slo_y[2]]))
-                assert np.all(np.isnan(buf[slo_x[0],slo_y[2]]))
+        for ghost_mask in GhostMask.all:
+            dF0.initialize(__random_init, exchange_kwds=dict(ghost_mask=ghost_mask))
+            dF1.initialize(__random_init, exchange_kwds=dict(ghost_mask=ghost_mask))
+
+            data = dF0.data + dF1.data
+            for (i,d) in enumerate(data):
+                print '    *buffer {}'.format(i)
+                buf = d.get().handle
+                so_x = (slice(0,nghosts[0]),
+                         slice(nghosts[0],nghosts[0]+npts[0]-1),
+                         slice(nghosts[0]+npts[0]-1,None))
+                so_y = (slice(0,nghosts[1]),
+                         slice(nghosts[1],nghosts[1]+npts[1]-1),
+                         slice(nghosts[1]+npts[1]-1,None))
+                si_x = (slice(nghosts[0], 2*nghosts[0]),
+                         slice(2*nghosts[0],npts[0]-1),
+                         slice(npts[0]-1, nghosts[0]+npts[0]-1))
+                si_y = (slice(nghosts[1], 2*nghosts[1]),
+                         slice(2*nghosts[1],npts[1]-1),
+                         slice(npts[1]-1, nghosts[1]+npts[1]-1))
+                assert buf.shape==(npts[0]+2*nghosts[0]-1,npts[1]+2*nghosts[1]-1)
+                assert buf[so_x[0],so_y[0]].shape == nghosts
+                assert buf[so_x[2],so_y[2]].shape == nghosts
+                assert buf[so_x[1],so_y[1]].shape == tuple([_-1 for _ in npts])
+                assert buf[si_x[0],si_y[0]].shape == nghosts
+                assert buf[si_x[2],si_y[2]].shape == nghosts
+                assert np.all(buf[so_x[0],so_y[1]] == buf[si_x[2],so_y[1]])
+                assert np.all(buf[so_x[2],so_y[1]] == buf[si_x[0],so_y[1]])
+                assert np.all(buf[so_x[1],so_y[0]] == buf[so_x[1],si_y[2]])
+                assert np.all(buf[so_x[1],so_y[2]] == buf[so_x[1],si_y[0]])
+                if (ghost_mask is GhostMask.FULL):
+                    assert np.all(buf[so_x[0],so_y[0]]==buf[si_x[2],si_y[2]])
+                    assert np.all(buf[so_x[2],so_y[0]]==buf[si_x[0],si_y[2]])
+                    assert np.all(buf[so_x[2],so_y[2]]==buf[si_x[0],si_y[0]])
+                    assert np.all(buf[so_x[0],so_y[2]]==buf[si_x[2],si_y[0]])
+                elif (ghost_mask is GhostMask.CROSS):
+                    assert np.all(np.isnan(buf[so_x[0],so_y[0]]))
+                    assert np.all(np.isnan(buf[so_x[2],so_y[0]]))
+                    assert np.all(np.isnan(buf[so_x[2],so_y[2]]))
+                    assert np.all(np.isnan(buf[so_x[0],so_y[2]]))
+                else:
+                    msg='Unknown ghost mask {}.'.format(ghost_mask)
+                    raise NotImplementedError(msg)
 
 
 def test_serial_initialization_3d():
@@ -146,66 +151,70 @@ def test_serial_initialization_3d():
         print '  Topo {}::{} | {}'.format(topo.full_tag, topo.backend.kind, topo.backend)
         dF0 = F0.discretize(topo)
         dF1 = F1.discretize(topo)
-        dF0.initialize(__random_init)
-        dF1.initialize(__random_init)
-
-        data = dF0.data + dF1.data
-        for (i,d) in enumerate(data):
-            print '    *buffer {}'.format(i)
-            buf = d.get().handle
-            slo_x = (slice(0,nghosts[0]),
-                     slice(nghosts[0],nghosts[0]+npts[0]-1),
-                     slice(nghosts[0]+npts[0]-1,None))
-            slo_y = (slice(0,nghosts[1]),
-                     slice(nghosts[1],nghosts[1]+npts[1]-1),
-                     slice(nghosts[1]+npts[1]-1,None))
-            slo_z = (slice(0,nghosts[2]),
-                     slice(nghosts[2],nghosts[2]+npts[2]-1),
-                     slice(nghosts[2]+npts[2]-1,None))
-            sli_x = (slice(nghosts[0], 2*nghosts[0]),
-                     slice(2*nghosts[0],npts[0]-1),
-                     slice(npts[0]-1, nghosts[0]+npts[0]-1))
-            sli_y = (slice(nghosts[1], 2*nghosts[1]),
-                     slice(2*nghosts[1],npts[1]-1),
-                     slice(npts[1]-1, nghosts[1]+npts[1]-1))
-            sli_z = (slice(nghosts[2], 2*nghosts[2]),
-                     slice(2*nghosts[2],npts[2]-1),
-                     slice(npts[2]-1, nghosts[2]+npts[2]-1))
-            assert buf.shape==(npts[0]+2*nghosts[0]-1,npts[1]+2*nghosts[1]-1,npts[2]
-                                +2*nghosts[2]-1)
-            assert buf[slo_x[0],slo_y[0],slo_z[0]].shape == nghosts
-            assert buf[slo_x[2],slo_y[2],slo_z[0]].shape == nghosts
-            assert buf[slo_x[0],slo_y[0],slo_z[2]].shape == nghosts
-            assert buf[slo_x[2],slo_y[2],slo_z[2]].shape == nghosts
-            assert buf[slo_x[1],slo_y[1],slo_z[1]].shape == tuple([_-1 for _ in npts])
-            assert buf[sli_x[0],sli_y[0],sli_z[0]].shape == nghosts
-            assert buf[sli_x[2],sli_y[2],sli_z[0]].shape == nghosts
-            assert buf[sli_x[0],sli_y[0],sli_z[2]].shape == nghosts
-            assert buf[sli_x[2],sli_y[2],sli_z[2]].shape == nghosts
-            assert np.all(buf[slo_x[0],slo_y[1],slo_z[1]] == buf[sli_x[2],slo_y[1],slo_z[1]])
-            assert np.all(buf[slo_x[2],slo_y[1],slo_z[1]] == buf[sli_x[0],slo_y[1],slo_z[1]])
-            assert np.all(buf[slo_x[1],slo_y[0],slo_z[1]] == buf[slo_x[1],sli_y[2],slo_z[1]])
-            assert np.all(buf[slo_x[1],slo_y[2],slo_z[1]] == buf[slo_x[1],sli_y[0],slo_z[1]])
-            assert np.all(buf[slo_x[1],slo_y[1],slo_z[0]] == buf[slo_x[1],slo_y[1],sli_z[2]])
-            assert np.all(buf[slo_x[1],slo_y[1],slo_z[2]] == buf[slo_x[1],slo_y[1],sli_z[0]])
-            if topo.mesh.isExchInclDiag:
-                assert np.all(buf[slo_x[0],slo_y[0],slo_z[0]]==buf[sli_x[2],sli_y[2],sli_z[2]])
-                assert np.all(buf[slo_x[2],slo_y[0],slo_z[0]]==buf[sli_x[0],sli_y[2],sli_z[2]])
-                assert np.all(buf[slo_x[2],slo_y[2],slo_z[0]]==buf[sli_x[0],sli_y[0],sli_z[2]])
-                assert np.all(buf[slo_x[0],slo_y[2],slo_z[0]]==buf[sli_x[2],sli_y[0],sli_z[2]])
-                assert np.all(buf[slo_x[0],slo_y[0],slo_z[2]]==buf[sli_x[2],sli_y[2],sli_z[0]])
-                assert np.all(buf[slo_x[2],slo_y[0],slo_z[2]]==buf[sli_x[0],sli_y[2],sli_z[0]])
-                assert np.all(buf[slo_x[2],slo_y[2],slo_z[2]]==buf[sli_x[0],sli_y[0],sli_z[0]])
-                assert np.all(buf[slo_x[0],slo_y[2],slo_z[2]]==buf[sli_x[2],sli_y[0],sli_z[0]])
-            else:
-                assert np.all(np.isnan(buf[slo_x[0],slo_y[0],slo_z[0]]))
-                assert np.all(np.isnan(buf[slo_x[2],slo_y[0],slo_z[0]]))
-                assert np.all(np.isnan(buf[slo_x[2],slo_y[2],slo_z[0]]))
-                assert np.all(np.isnan(buf[slo_x[0],slo_y[2],slo_z[0]]))
-                assert np.all(np.isnan(buf[slo_x[0],slo_y[0],slo_z[2]]))
-                assert np.all(np.isnan(buf[slo_x[2],slo_y[0],slo_z[2]]))
-                assert np.all(np.isnan(buf[slo_x[2],slo_y[2],slo_z[2]]))
-                assert np.all(np.isnan(buf[slo_x[0],slo_y[2],slo_z[2]]))
+        for ghost_mask in GhostMask.all:
+            dF0.initialize(__random_init, exchange_kwds=dict(ghost_mask=ghost_mask))
+            dF1.initialize(__random_init, exchange_kwds=dict(ghost_mask=ghost_mask))
+
+            data = dF0.data + dF1.data
+            for (i,d) in enumerate(data):
+                print '    *buffer {}'.format(i)
+                buf = d.get().handle
+                so_x = (slice(0,nghosts[0]),
+                         slice(nghosts[0],nghosts[0]+npts[0]-1),
+                         slice(nghosts[0]+npts[0]-1,None))
+                so_y = (slice(0,nghosts[1]),
+                         slice(nghosts[1],nghosts[1]+npts[1]-1),
+                         slice(nghosts[1]+npts[1]-1,None))
+                so_z = (slice(0,nghosts[2]),
+                         slice(nghosts[2],nghosts[2]+npts[2]-1),
+                         slice(nghosts[2]+npts[2]-1,None))
+                si_x = (slice(nghosts[0], 2*nghosts[0]),
+                         slice(2*nghosts[0],npts[0]-1),
+                         slice(npts[0]-1, nghosts[0]+npts[0]-1))
+                si_y = (slice(nghosts[1], 2*nghosts[1]),
+                         slice(2*nghosts[1],npts[1]-1),
+                         slice(npts[1]-1, nghosts[1]+npts[1]-1))
+                si_z = (slice(nghosts[2], 2*nghosts[2]),
+                         slice(2*nghosts[2],npts[2]-1),
+                         slice(npts[2]-1, nghosts[2]+npts[2]-1))
+                assert buf.shape==(npts[0]+2*nghosts[0]-1,npts[1]+2*nghosts[1]-1,npts[2]
+                                    +2*nghosts[2]-1)
+                assert buf[so_x[0],so_y[0],so_z[0]].shape == nghosts
+                assert buf[so_x[2],so_y[2],so_z[0]].shape == nghosts
+                assert buf[so_x[0],so_y[0],so_z[2]].shape == nghosts
+                assert buf[so_x[2],so_y[2],so_z[2]].shape == nghosts
+                assert buf[so_x[1],so_y[1],so_z[1]].shape == tuple([_-1 for _ in npts])
+                assert buf[si_x[0],si_y[0],si_z[0]].shape == nghosts
+                assert buf[si_x[2],si_y[2],si_z[0]].shape == nghosts
+                assert buf[si_x[0],si_y[0],si_z[2]].shape == nghosts
+                assert buf[si_x[2],si_y[2],si_z[2]].shape == nghosts
+                assert np.all(buf[so_x[0],so_y[1],so_z[1]] == buf[si_x[2],so_y[1],so_z[1]])
+                assert np.all(buf[so_x[2],so_y[1],so_z[1]] == buf[si_x[0],so_y[1],so_z[1]])
+                assert np.all(buf[so_x[1],so_y[0],so_z[1]] == buf[so_x[1],si_y[2],so_z[1]])
+                assert np.all(buf[so_x[1],so_y[2],so_z[1]] == buf[so_x[1],si_y[0],so_z[1]])
+                assert np.all(buf[so_x[1],so_y[1],so_z[0]] == buf[so_x[1],so_y[1],si_z[2]])
+                assert np.all(buf[so_x[1],so_y[1],so_z[2]] == buf[so_x[1],so_y[1],si_z[0]])
+                if (ghost_mask is GhostMask.FULL):
+                    assert np.all(buf[so_x[0],so_y[0],so_z[0]]==buf[si_x[2],si_y[2],si_z[2]])
+                    assert np.all(buf[so_x[2],so_y[0],so_z[0]]==buf[si_x[0],si_y[2],si_z[2]])
+                    assert np.all(buf[so_x[2],so_y[2],so_z[0]]==buf[si_x[0],si_y[0],si_z[2]])
+                    assert np.all(buf[so_x[0],so_y[2],so_z[0]]==buf[si_x[2],si_y[0],si_z[2]])
+                    assert np.all(buf[so_x[0],so_y[0],so_z[2]]==buf[si_x[2],si_y[2],si_z[0]])
+                    assert np.all(buf[so_x[2],so_y[0],so_z[2]]==buf[si_x[0],si_y[2],si_z[0]])
+                    assert np.all(buf[so_x[2],so_y[2],so_z[2]]==buf[si_x[0],si_y[0],si_z[0]])
+                    assert np.all(buf[so_x[0],so_y[2],so_z[2]]==buf[si_x[2],si_y[0],si_z[0]])
+                elif (ghost_mask is GhostMask.CROSS):
+                    assert np.all(np.isnan(buf[so_x[0],so_y[0],so_z[0]]))
+                    assert np.all(np.isnan(buf[so_x[2],so_y[0],so_z[0]]))
+                    assert np.all(np.isnan(buf[so_x[2],so_y[2],so_z[0]]))
+                    assert np.all(np.isnan(buf[so_x[0],so_y[2],so_z[0]]))
+                    assert np.all(np.isnan(buf[so_x[0],so_y[0],so_z[2]]))
+                    assert np.all(np.isnan(buf[so_x[2],so_y[0],so_z[2]]))
+                    assert np.all(np.isnan(buf[so_x[2],so_y[2],so_z[2]]))
+                    assert np.all(np.isnan(buf[so_x[0],so_y[2],so_z[2]]))
+                else:
+                    msg='Unknown ghost mask {}.'.format(ghost_mask)
+                    raise NotImplementedError(msg)
 
 
 def iter_backends():
@@ -263,7 +272,8 @@ def test_mpi_ghost_exchange(comm):
                     def ghost_vals(shape,dtype,i,d,rank,local_dir):
                         if (shape is None) or (len(shape)==0):
                             raise ValueError('Shape is None or an empty tuple: {}'.format(shape))
-                        base = np.full(shape=shape, dtype=dtype, fill_value=ghost_base(i,d,rank,local_dir))
+                        base = np.full(shape=shape, dtype=dtype, fill_value=ghost_base(i, d, rank,
+                                                                                local_dir))
                         I = np.ix_(*tuple(np.arange(shape[direction], dtype=dtype) \
                                                 for direction in xrange(len(shape))))
                         return base + I[d]
@@ -305,12 +315,16 @@ def test_mpi_ghost_exchange(comm):
                                     rdata = data[rghosts]
                                     if (ldata is not None):
                                         ldata = np.atleast_1d(ldata.get())
-                                        target_vals = ghost_vals(ldata.shape, dtype,i,d,left_rank,1)
-                                        assert np.allclose(ldata, target_vals), (rank, target_vals)
+                                        target_vals = ghost_vals(ldata.shape, dtype, i, d,
+                                                left_rank,1)
+                                        assert np.allclose(ldata, target_vals), (rank, 
+                                                                                    target_vals)
                                     if (rdata is not None):
                                         rdata = np.atleast_1d(rdata.get())
-                                        target_vals = ghost_vals(rdata.shape, dtype, i,d,right_rank,0)
-                                        assert np.allclose(rdata, target_vals), (rank, target_vals)
+                                        target_vals = ghost_vals(rdata.shape, dtype, i, d, 
+                                                                                right_rank, 0)
+                                        assert np.allclose(rdata, target_vals), (rank, 
+                                                                                    target_vals)
                             if rank==0:
                                 print
 
@@ -421,7 +435,8 @@ def test_mpi_ghost_accumulate(comm):
 
                                     dF.exchange_ghosts(directions=directions,
                                                        exchange_method=exchange_method,
-                                                       ghost_op=GhostOperation.ACCUMULATE)
+                                                       ghost_op=GhostOperation.ACCUMULATE,
+                                                       ghost_mask=GhostMask.CROSS)
 
                                     for (i,data) in enumerate(dF.data):
                                         for displacements in all_displacements:
diff --git a/hysop/mesh/cartesian_mesh.py b/hysop/mesh/cartesian_mesh.py
index 2fc95059f..ffac96902 100644
--- a/hysop/mesh/cartesian_mesh.py
+++ b/hysop/mesh/cartesian_mesh.py
@@ -10,7 +10,7 @@ See also
 * :ref:`topologies` in HySoP user guide.
 """
 
-from hysop.constants import np, BoundaryCondition, HYSOP_INTEGER, HYSOP_REAL
+from hysop.constants import np, BoundaryCondition, GhostMask, HYSOP_INTEGER, HYSOP_REAL
 from hysop.deps import it
 from hysop.tools.decorators import debug
 from hysop.tools.numpywrappers import npw
@@ -22,8 +22,7 @@ class CartesianMeshView(MeshView):
     """Mesh views on cartesian meshes. """
 
     __slots__ = ('_mesh', '_topology_state')
-    isExchInclDiag = True  # Is ghosts exchange includes diagonal ghosts
-
+  
     @debug
     def __new__(cls, mesh, topology_state, **kwds):
         """Create a CartesianMeshView."""
@@ -215,29 +214,37 @@ class CartesianMeshView(MeshView):
         """
         return self.__get_transposed_mesh_attr('_local_compute_slices')
 
-    def get_local_inner_ghost_slices(self, ghosts=None):
+    def get_local_inner_ghost_slices(self, ghosts=None, ghost_mask=GhostMask.CROSS):
         """
-        Return a list of slices defining the ghosts in this arrays as local indices.
-        Those slices corresponds to neighbour processes overlap on local
-        compute slices (depending on isExchInclDiag parameter it INCLUDES
-        or EXCLUDE diagonal processes).
-        For each direction, a pair of tuples of slices is returned, one
-        for left and one for the right ghosts.
+        Return a list of slices defining the ghosts in this arrays a local indices.
+        Those slices corresponds to neighbour processes overlap on local compute slices.
+        Depending on the ghost_mask parameter, one can include or exclude diagonal ghosts: 
+            GhostMask.FULL:   INCLUDES diagonal ghosts
+            GhostMask.CROSS:  EXCLUDES diagonal ghosts
+        For each direction, a pair of tuples of slices is returned, one for left and one 
+        for the right ghosts.
         """
         dim = self.dim
         local_start = self.local_start
         local_stop  = self.local_stop
         compute_resolution = self.compute_resolution
-        local_resolution = self.local_resolution
+        local_resolution   = self.local_resolution
 
         ghosts = to_tuple(ghosts or self.ghosts)
         ghosts = ghosts*dim if len(ghosts)==1 else ghosts
         assert len(ghosts)==dim
-
-        if self.isExchInclDiag:
+        
+        if (ghost_mask is GhostMask.FULL):
             gh_slices = lambda j: slice(None)
-        else:
+            resolution = local_resolution
+        elif (ghost_mask is GhostMask.CROSS):
             gh_slices = lambda j: slice(ghosts[j], compute_resolution[j]+ghosts[j])
+            resolution = compute_resolution
+        else:
+            msg_mask='Unknown ghost ghost_mask configuration {}.'
+            msg_mask=msg_mask.format(ghost_mask)
+            raise NotImplementedError(msg_mask)
+
         local_inner_ghost_slices = []
         for i in xrange(dim):
             inner_lslices = tuple(
@@ -248,37 +255,43 @@ class CartesianMeshView(MeshView):
                 gh_slices(j) if j!=i else
                 slice(local_stop[i]-ghosts[i], local_stop[i])
                 for j in xrange(dim) )
-            if self.isExchInclDiag:
-                inner_shape = local_resolution.copy()
-            else:
-                inner_shape = compute_resolution.copy()
+            inner_shape = resolution.copy()
             inner_shape[i] = ghosts[i]
             local_inner_ghost_slices.append( (inner_lslices, inner_rslices, inner_shape) )
 
         return tuple(local_inner_ghost_slices)
 
-    def get_local_outer_ghost_slices(self, ghosts=None):
+    def get_local_outer_ghost_slices(self, ghosts=None, ghost_mask=GhostMask.CROSS):
         """
-        Return a list of slices defining the ghosts in this arrays as local indices.
-        Those slices corresponds to local to process ghosts (depending
-        on isExchInclDiag it INCLUDES or EXCLUDES diagonal processes).
-        For each direction, a pair of tuples of slices is returned, one
-        for left and one for the right ghosts.
+        Return a list of slices defining the ghosts in this arrays a local indices.
+        Those slices corresponds to local to process ghost slices.
+        Depending on the ghost_mask parameter, one can include or exclude diagonal ghosts: 
+            GhostMask.FULL:   INCLUDES diagonal ghosts
+            GhostMask.CROSS:  EXCLUDES diagonal ghosts
+        For each direction, a pair of tuples of slices is returned, one for left and one 
+        for the right ghosts.
         """
         dim = self.dim
         local_start = self.local_start
         local_stop  = self.local_stop
         compute_resolution = self.compute_resolution
-        local_resolution = self.local_resolution
+        local_resolution   = self.local_resolution
 
         ghosts = to_tuple(ghosts or self.ghosts)
         ghosts = ghosts*dim if len(ghosts)==1 else ghosts
         assert len(ghosts)==dim
 
-        if self.isExchInclDiag:
+        if (ghost_mask is GhostMask.FULL):
             gh_slices = lambda j: slice(None)
-        else:
+            resolution = local_resolution
+        elif (ghost_mask is GhostMask.CROSS):
             gh_slices = lambda j: slice(ghosts[j], compute_resolution[j]+ghosts[j])
+            resolution = compute_resolution
+        else:
+            msg_mask='Unknown ghost ghost_mask configuration {}.'
+            msg_mask=msg_mask.format(ghost_mask)
+            raise NotImplementedError(msg_mask)
+
         local_outer_ghost_slices = []
         for i in xrange(dim):
             outer_lslices = tuple(
@@ -289,10 +302,7 @@ class CartesianMeshView(MeshView):
                 gh_slices(j) if j!=i else
                 slice(local_stop[i], local_stop[i]+ghosts[i])
                 for j in xrange(dim) )
-            if self.isExchInclDiag:
-                outer_shape = local_resolution.copy()
-            else:
-                outer_shape    = compute_resolution.copy()
+            outer_shape = resolution.copy()
             outer_shape[i] = ghosts[i]
             local_outer_ghost_slices.append( (outer_lslices, outer_rslices, outer_shape) )
 
-- 
GitLab