Repository: devitocodes/devito
Issue: Potentially wrong `shape_with_halo` with HALO > DOMAIN

This might cause NaNs or segfaults by essentially breaking the contract between python-land and generated code

Reproducer here: https://github.com/devitocodes/devito/compare/fix-shape-with-halo

For an Operator that produces NaNs, take [this test](https://github.com/devitocodes/devito/blob/consolidate-aliases-3/tests/test_dse.py#L582
) and change the grid such that, in at least one decomposed Dimension, it's smaller than the halo size, for example (3, 3, 3)
it will produce garbage values -- but you have to run it many times until it produces NaNs, thus causing np.close to fail ! or make the test "more reasonable" such that it produces a given expected value against which you can assert deterministically 

Changed files:
devito/types/dimension.py, devito/types/dense.py, devito/data/decomposition.py, tests/test_data.py, devito/mpi/distributed.py

diff --git a/devito/data/decomposition.py b/devito/data/decomposition.py
index bb4d4343f..3f291f85a 100644
--- a/devito/data/decomposition.py
+++ b/devito/data/decomposition.py
@@ -144,10 +144,10 @@ class Decomposition(tuple):
             There are three possible cases:
             * int. Given ``I``, a global index, return the corresponding
               relative local index if ``I`` belongs to the local subdomain,
-              None otherwise.
+              ``None`` otherwise.
             * int, DataSide. Given ``O`` and ``S``, respectively a global
               offset and a side, return the relative local offset. This
-              can be 0 if the local subdomain doesn't intersect with the
+              can be ``None`` if the local subdomain doesn't intersect with the
               region defined by the given global offset.
             * (int, int).  Given global ``(min, max)``, return ``(min', max')``
               representing the corresponding relative local min/max. If the
@@ -313,16 +313,20 @@ class Decomposition(tuple):
         elif len(args) == 2:
             # index_glb_to_loc(offset, side)
             if self.loc_empty:
-                return 0
-            rel_ofs, side = args
+                return None
+            abs_ofs, side = args
             if side is LEFT:
-                abs_ofs = self.glb_min + rel_ofs
-                size = self.loc_abs_max - base + 1
-                return min(abs_ofs - base, size) if abs_ofs > base else 0
+                rel_ofs = self.glb_min + abs_ofs - base
+                if abs_ofs >= base and abs_ofs <= top:
+                    return rel_ofs
+                else:
+                    return None
             else:
-                abs_ofs = self.glb_max - rel_ofs
-                size = top - self.loc_abs_min + 1
-                return min(top - abs_ofs, size) if abs_ofs < top else 0
+                rel_ofs = abs_ofs - (self.glb_max - top)
+                if abs_ofs >= self.glb_max - top and abs_ofs <= self.glb_max - base:
+                    return rel_ofs
+                else:
+                    return None
         else:
             raise TypeError("Expected 1 or 2 arguments, found %d" % len(args))
 
diff --git a/devito/mpi/distributed.py b/devito/mpi/distributed.py
index 5c0f0d02b..4f75d9443 100644
--- a/devito/mpi/distributed.py
+++ b/devito/mpi/distributed.py
@@ -245,6 +245,12 @@ class Distributor(AbstractDistributor):
     def topology(self):
         return self._topology
 
+    @cached_property
+    def is_boundary_rank(self):
+        """ MPI rank interfaces with the boundary of the domain. """
+        return any([True if i == 0 or i == j-1 else False for i, j in
+                   zip(self.mycoords, self.topology)])
+
     @cached_property
     def all_coords(self):
         """
diff --git a/devito/types/dense.py b/devito/types/dense.py
index 33e608070..02e0764ef 100644
--- a/devito/types/dense.py
+++ b/devito/types/dense.py
@@ -308,16 +308,29 @@ class DiscreteFunction(AbstractFunction, ArgProvider):
     @cached_property
     def _size_outhalo(self):
         """Number of points in the outer halo region."""
+
         if self._distributor is None:
+            # Computational domain is not distributed and hence the outhalo
+            # and inhalo correspond
             return self._size_inhalo
 
-        left = [self._distributor.glb_to_loc(d, i, LEFT, strict=False)
-                for d, i in zip(self.dimensions, self._size_inhalo.left)]
-        right = [self._distributor.glb_to_loc(d, i, RIGHT, strict=False)
-                 for d, i in zip(self.dimensions, self._size_inhalo.right)]
+        left = [abs(min(i.loc_abs_min-i.glb_min-j, 0)) if i and not i.loc_empty else 0
+                for i, j in zip(self._decomposition, self._size_inhalo.left)]
+        right = [max(i.loc_abs_max+j-i.glb_max, 0) if i and not i.loc_empty else 0
+                 for i, j in zip(self._decomposition, self._size_inhalo.right)]
 
         sizes = tuple(Size(i, j) for i, j in zip(left, right))
 
+        if self._distributor.is_parallel and any(left) > 0 or any(right) > 0:
+            try:
+                if not self._distributor.is_boundary_rank:
+                    warning("A space order of %d  and a halo size of %d has been set "
+                            "but the current rank (%d) has a domain size of only %d" %
+                            (self._space_order, max(self._size_inhalo),
+                             self._distributor.myrank, min(self.grid.shape_local)))
+            except AttributeError:
+                pass
+
         return DimensionTuple(*sizes, getters=self.dimensions, left=left, right=right)
 
     @property
diff --git a/devito/types/dimension.py b/devito/types/dimension.py
index 8e2f0b1f7..1c710322f 100644
--- a/devito/types/dimension.py
+++ b/devito/types/dimension.py
@@ -615,8 +615,10 @@ class SubDimension(DerivedDimension):
     def _arg_defaults(self, grid=None, **kwargs):
         if grid is not None and grid.is_distributed(self.root):
             # Get local thickness
-            ltkn = grid.distributor.glb_to_loc(self.root, self.thickness.left[1], LEFT)
-            rtkn = grid.distributor.glb_to_loc(self.root, self.thickness.right[1], RIGHT)
+            ltkn = grid.distributor.glb_to_loc(self.root,
+                                               self.thickness.left[1], LEFT) or 0
+            rtkn = grid.distributor.glb_to_loc(self.root,
+                                               self.thickness.right[1], RIGHT) or 0
             return {i.name: v for i, v in zip(self._thickness_map, (ltkn, rtkn))}
         else:
             return {k.name: v for k, v in self.thickness}
diff --git a/tests/test_data.py b/tests/test_data.py
index 4304c8241..05bc32159 100644
--- a/tests/test_data.py
+++ b/tests/test_data.py
@@ -258,6 +258,16 @@ class TestMetaData(object):
         assert tuple(i + j*2 for i, j in zip(u.shape, u._size_halo.left)) ==\\
             u.shape_with_halo
 
+        # Try with different grid shape and space_order
+        grid2 = Grid(shape=(3, 3, 3))
+        u2 = Function(name='u2', grid=grid2, space_order=4, padding=0)
+        assert u2.shape == (3, 3, 3)
+        assert u2._offset_domain == (4, 4, 4)
+        assert u2._offset_halo == ((0, 7), (0, 7), (0, 7))
+        assert tuple(i + j*2 for i, j in zip(u2.shape, u2._size_halo.left)) ==\\
+            u2.shape_with_halo
+        assert u2.shape_with_halo == (11, 11, 11)
+
     def test_wo_halo_w_padding(self):
         grid = Grid(shape=(4, 4, 4))
         u = Function(name='u', grid=grid, space_order=2, padding=((1, 1), (3, 3), (4, 4)))
@@ -342,6 +352,19 @@ class TestDecomposition(object):
         assert d.index_glb_to_loc((1, 6), rel=False) == (5, 6)
         assert d.index_glb_to_loc((None, None), rel=False) == (5, 7)
 
+    def test_glb_to_loc_w_side(self):
+        d = Decomposition([[0, 1, 2], [3, 4], [5, 6, 7], [8, 9, 10, 11]], 2)
+
+        # A global index as single argument
+        assert d.index_glb_to_loc(5, LEFT) == 0
+        assert d.index_glb_to_loc(6, RIGHT) == 2
+        assert d.index_glb_to_loc(7, LEFT) == 2
+        assert d.index_glb_to_loc(4, RIGHT) == 0
+        assert d.index_glb_to_loc(6, LEFT) == 1
+        assert d.index_glb_to_loc(5, RIGHT) == 1
+        assert d.index_glb_to_loc(2, LEFT) is None
+        assert d.index_glb_to_loc(3, RIGHT) is None
+
     def test_loc_to_glb_index_conversions(self):
         d = Decomposition([[0, 1, 2], [3, 4], [5, 6, 7], [8, 9, 10, 11]], 2)