diff --git a/devito/__init__.py b/devito/__init__.py
index 535b93ecbb..e91a4a9726 100644
--- a/devito/__init__.py
+++ b/devito/__init__.py
@@ -96,6 +96,15 @@ def reinit_compiler(val):
 # optimisations.
 configuration.add('safe-math', 0, [0, 1], preprocessor=bool, callback=reinit_compiler)
 
+# By default Devito allows some extra halo space to satisfy derivatives
+# which may result in OOB accesses. However, in most of cases these extra halo points
+# are redundant. In single-node simulations the space overhead is usually minimal
+# compared to the computational domain, however in multi-node simulations, this may
+# result in redundant ghost-cell communications, which may have high a communication
+# cost impact.
+preprocessor = lambda i: {0: 'None', 1: 'HALF', 2: 'CUSTOM'}.get(i, i)
+configuration.add('safe-halo', 0, [0, 1, 2], preprocessor=preprocessor)
+
 # Enable/disable automatic padding for allocated data
 configuration.add('autopadding', False, [False, True])
 
diff --git a/devito/ir/equations/algorithms.py b/devito/ir/equations/algorithms.py
index c27ef5c54f..7350dabd41 100644
--- a/devito/ir/equations/algorithms.py
+++ b/devito/ir/equations/algorithms.py
@@ -121,8 +121,9 @@ def lower_exprs(expressions, **kwargs):
             f = i.function
 
             # Introduce shifting to align with the computational domain
-            indices = [(lower_exprs(a) + o) for a, o in
-                       zip(i.indices, f._size_nodomain.left)]
+            indices = []
+            for a, o in zip(i.indices, f._size_nodomain.left):
+                indices.append(lower_exprs(a) + o)
 
             # Substitute spacing (spacing only used in own dimension)
             indices = [i.xreplace({d.spacing: 1, -d.spacing: -1})
diff --git a/devito/parameters.py b/devito/parameters.py
index 885d8a56e4..5edbfb4d6e 100644
--- a/devito/parameters.py
+++ b/devito/parameters.py
@@ -152,7 +152,8 @@ def _signature_items(self):
     'DEVITO_FIRST_TOUCH': 'first-touch',
     'DEVITO_JIT_BACKDOOR': 'jit-backdoor',
     'DEVITO_IGNORE_UNKNOWN_PARAMS': 'ignore-unknowns',
-    'DEVITO_SAFE_MATH': 'safe-math'
+    'DEVITO_SAFE_MATH': 'safe-math',
+    'DEVITO_SAFE_HALO': 'safe-halo',
 }
 
 env_vars_deprecated = {
diff --git a/devito/types/dense.py b/devito/types/dense.py
index 08146eaf83..b8f18b6ea1 100644
--- a/devito/types/dense.py
+++ b/devito/types/dense.py
@@ -1113,7 +1113,12 @@ def __halo_setup__(self, **kwargs):
         else:
             space_order = kwargs.get('space_order', 1)
             if isinstance(space_order, int):
-                halo = (space_order, space_order)
+                if configuration['safe-halo'] == 'HALF':
+                    halo = (int(space_order/2), int(space_order/2))
+                elif configuration['safe-halo'] == 'CUSTOM':
+                    halo = (int(space_order/2) + 1, int(space_order/2) + 1)
+                else:
+                    halo = (space_order, space_order)
             elif isinstance(space_order, tuple) and len(space_order) == 3:
                 _, left_points, right_points = space_order
                 halo = (left_points, right_points)
diff --git a/examples/seismic/acoustic/acoustic_example.py b/examples/seismic/acoustic/acoustic_example.py
index 7e922cbb6e..a07717660c 100644
--- a/examples/seismic/acoustic/acoustic_example.py
+++ b/examples/seismic/acoustic/acoustic_example.py
@@ -38,6 +38,7 @@ def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
     save = full_run and not checkpointing
     # Define receiver geometry (spread across x, just below surface)
     rec, u, summary = solver.forward(save=save, autotune=autotune)
+    print("Norm u:", norm(u))
 
     if preset == 'constant-isotropic':
         # With a new m as Constant
diff --git a/examples/seismic/elastic/elastic_example.py b/examples/seismic/elastic/elastic_example.py
index bd1845baaa..54072a2129 100644
--- a/examples/seismic/elastic/elastic_example.py
+++ b/examples/seismic/elastic/elastic_example.py
@@ -29,6 +29,7 @@ def run(shape=(50, 50), spacing=(20.0, 20.0), tn=1000.0,
     info("Applying Forward")
     # Define receiver geometry (spread across x, just below surface)
     rec1, rec2, v, tau, summary = solver.forward(autotune=autotune)
+    print("Norm v:", norm(v[0]))
     return (summary.gflopss, summary.oi, summary.timings,
             [rec1, rec2, v, tau])