diff --git a/devito/data/decomposition.py b/devito/data/decomposition.py
index 7fad455d33..b111fe8217 100644
--- a/devito/data/decomposition.py
+++ b/devito/data/decomposition.py
@@ -320,7 +320,7 @@ def index_glb_to_loc(self, *args, rel=True):
             if self.loc_empty:
                 return None
             abs_ofs, side = args
-            if side is LEFT:
+            if side == LEFT:
                 rel_ofs = self.glb_min + abs_ofs - base
                 if abs_ofs >= base and abs_ofs <= top:
                     return rel_ofs
diff --git a/devito/ir/equations/algorithms.py b/devito/ir/equations/algorithms.py
index c8f207855b..e0281e1139 100644
--- a/devito/ir/equations/algorithms.py
+++ b/devito/ir/equations/algorithms.py
@@ -8,7 +8,7 @@
                           ConditionalDimension)
 from devito.types.array import Array
 from devito.types.basic import AbstractFunction
-from devito.types.dimension import MultiSubDimension
+from devito.types.dimension import MultiSubDimension, Thickness
 from devito.data.allocators import DataReference
 from devito.logger import warning
 
@@ -213,18 +213,24 @@ def _(expr, mapper, rebuilt, sregistry):
     _concretize_subdims(expr.implicit_dims, mapper, rebuilt, sregistry)
 
 
+@_concretize_subdims.register(Thickness)
+def _(tkn, mapper, rebuilt, sregistry):
+    if tkn in mapper:
+        # Already have a substitution for this thickness
+        return
+
+    mapper[tkn] = tkn._rebuild(name=sregistry.make_name(prefix=tkn.name))
+
+
 @_concretize_subdims.register(SubDimension)
 def _(d, mapper, rebuilt, sregistry):
     if d in mapper:
         # Already have a substitution for this dimension
         return
 
-    tkns_subs = {tkn: tkn._rebuild(name=sregistry.make_name(prefix=tkn.name))
-                 for tkn in d.tkns}
-    left, right = [mM.subs(tkns_subs) for mM in (d.symbolic_min, d.symbolic_max)]
-    thickness = tuple((v, d._thickness_map[k]) for k, v in tkns_subs.items())
-
-    mapper[d] = d._rebuild(symbolic_min=left, symbolic_max=right, thickness=thickness)
+    tkns = tuple(t._rebuild(name=sregistry.make_name(prefix=t.name)) for t in d.tkns)
+    mapper.update({tkn0: tkn1 for tkn0, tkn1 in zip(d.tkns, tkns)})
+    mapper[d] = d._rebuild(thickness=tkns)
 
 
 @_concretize_subdims.register(ConditionalDimension)
@@ -261,11 +267,9 @@ def _(d, mapper, rebuilt, sregistry):
         # Already have a substitution for this dimension
         return
 
-    abstract_tkns = MultiSubDimension._symbolic_thickness(d.parent.name)
-    concrete_tkns = tuple(tkn._rebuild(name=sregistry.make_name(prefix=tkn.name))
-                          for tkn in abstract_tkns)
-
-    kwargs = {'thickness': concrete_tkns}
+    tkns = tuple(tkn._rebuild(name=sregistry.make_name(prefix=tkn.name))
+                 for tkn in d.thickness)
+    kwargs = {'thickness': tkns}
     fkwargs = {}
 
     idim0 = d.implicit_dimension
diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py
index 3d6bb3ba60..bb558c8d58 100644
--- a/devito/ir/support/basic.py
+++ b/devito/ir/support/basic.py
@@ -1376,7 +1376,7 @@ def disjoint_test(e0, e1, d, it):
     if d.is_Custom:
         subs = {}
     elif d.is_Sub and d.is_left:
-        subs = {d.root.symbolic_min: 0, **dict([d.thickness.left])}
+        subs = {d.root.symbolic_min: 0, d.ltkn: d.ltkn.value}
     else:
         return False
 
diff --git a/devito/mpi/halo_scheme.py b/devito/mpi/halo_scheme.py
index 0062b5b32f..7074b4c5d5 100644
--- a/devito/mpi/halo_scheme.py
+++ b/devito/mpi/halo_scheme.py
@@ -89,8 +89,7 @@ def __init__(self, exprs, ispace):
                              if d._defines & self.dimensions])
         subdims = [d for d in dims if d.is_Sub and not d.local]
         for i in subdims:
-            ltk, _ = i.thickness.left
-            rtk, _ = i.thickness.right
+            ltk, rtk = i.tkns
             self._honored[i.root] = frozenset([(ltk, rtk)])
         self._honored = frozendict(self._honored)
 
diff --git a/devito/operator/operator.py b/devito/operator/operator.py
index 24258ad671..9bfdcdc255 100644
--- a/devito/operator/operator.py
+++ b/devito/operator/operator.py
@@ -31,6 +31,8 @@
                           split, timed_pass, timed_region, contains_val)
 from devito.types import (Buffer, Grid, Evaluable, host_layer, device_layer,
                           disk_layer)
+from devito.types.dimension import Thickness
+
 
 __all__ = ['Operator']
 
@@ -640,6 +642,11 @@ def _prepare_arguments(self, autotune=None, **kwargs):
         for d in reversed(toposort):
             args.update(d._arg_values(self._dspace[d], grid, **kwargs))
 
+        # Process Thicknesses
+        for p in self.parameters:
+            if isinstance(p, Thickness):
+                args.update(p._arg_values(grid=grid, **kwargs))
+
         # Process Objects
         for o in self.objects:
             args.update(o._arg_values(grid=grid, **kwargs))
diff --git a/devito/types/dimension.py b/devito/types/dimension.py
index d8fbeff959..7e1f7fc811 100644
--- a/devito/types/dimension.py
+++ b/devito/types/dimension.py
@@ -9,7 +9,7 @@
 from devito.data import LEFT, RIGHT
 from devito.exceptions import InvalidArgument
 from devito.logger import debug
-from devito.tools import Pickable, is_integer, flatten
+from devito.tools import Pickable, is_integer, memoized_meth
 from devito.types.args import ArgProvider
 from devito.types.basic import Symbol, DataSymbol, Scalar
 from devito.types.constant import Constant
@@ -22,8 +22,7 @@
            'VirtualDimension', 'Spacing', 'dimensions']
 
 
-Thickness = namedtuple('Thickness', 'left right')
-SubDimensionOffset = namedtuple('SubDimensionOffset', 'value extreme thickness')
+SubDimensionThickness = namedtuple('SubDimensionThickness', 'left right')
 
 
 class Dimension(ArgProvider):
@@ -536,6 +535,74 @@ def _arg_check(self, *args, **kwargs):
 # The Dimensions below are exposed in the user API. They can only be created by
 # the user
 
+class Thickness(DataSymbol):
+    """A DataSymbol to represent a thickness of a SubDimension"""
+
+    __rkwargs__ = DataSymbol.__rkwargs__ + ('root', 'side', 'local', 'value')
+
+    def __new__(cls, *args, root=None, side=None, local=False, **kwargs):
+        newobj = super().__new__(cls, *args, **kwargs)
+        newobj._root = root
+        newobj._side = side
+        newobj._local = local
+
+        return newobj
+
+    def __init_finalize__(self, *args, **kwargs):
+        self._value = kwargs.pop('value', None)
+
+        kwargs.setdefault('is_const', True)
+        super().__init_finalize__(*args, **kwargs)
+
+    @property
+    def root(self):
+        return self._root
+
+    @property
+    def side(self):
+        return self._side
+
+    @property
+    def local(self):
+        return self._local
+
+    @property
+    def value(self):
+        return self._value
+
+    def _arg_check(self, *args, **kwargs):
+        pass
+
+    def _arg_values(self, grid=None, **kwargs):
+        # Allow override of thickness values to disable BCs
+        # However, arguments from the user are considered global
+        # So overriding the thickness to a nonzero value should not cause
+        # boundaries to exist between ranks where they did not before
+        rtkn = kwargs.get(self.name, self.value)
+        if grid is not None and grid.is_distributed(self.root):
+            # Get local thickness
+            if self.local:
+                # Dimension is of type `left`/`right` - compute the offset
+                # and then add 1 to get the appropriate thickness
+                if self.value is not None:
+                    tkn = grid.distributor.glb_to_loc(self.root, rtkn-1, self.side)
+                    tkn = tkn+1 if tkn is not None else 0
+                else:
+                    tkn = 0
+            else:
+                # Dimension is of type `middle`
+                tkn = grid.distributor.glb_to_loc(self.root, rtkn, self.side) or 0
+        else:
+            tkn = rtkn or 0
+
+        return {self.name: tkn}
+
+    def _arg_finalize(self, *args, **kwargs):
+        return {}
+
+    def _arg_apply(self, *args, **kwargs):
+        pass
+
 
 class AbstractSubDimension(DerivedDimension):
 
@@ -550,21 +617,31 @@ class AbstractSubDimension(DerivedDimension):
 
     is_AbstractSub = True
 
-    __rargs__ = (DerivedDimension.__rargs__ +
-                 ('symbolic_min', 'symbolic_max', 'thickness'))
+    __rargs__ = DerivedDimension.__rargs__ + ('thickness',)
     __rkwargs__ = ()
 
-    def __init_finalize__(self, name, parent, left, right, thickness, **kwargs):
+    _thickness_type = Symbol
+
+    def __init_finalize__(self, name, parent, thickness, **kwargs):
         super().__init_finalize__(name, parent)
-        self._interval = sympy.Interval(left, right)
-        self._thickness = Thickness(*thickness)
+        thickness = thickness or (None, None)
+        if any(isinstance(tkn, self._thickness_type) for tkn in thickness):
+            self._thickness = SubDimensionThickness(*thickness)
+        else:
+            self._thickness = self._symbolic_thickness(thickness=thickness)
 
-    @classmethod
-    def _symbolic_thickness(cls, name, stype=Scalar):
-        return (stype(name="%s_ltkn" % name, dtype=np.int32,
-                      is_const=True, nonnegative=True),
-                stype(name="%s_rtkn" % name, dtype=np.int32,
-                      is_const=True, nonnegative=True))
+    @cached_property
+    def _interval(self):
+        left = self.parent.symbolic_min + self.ltkn
+        right = self.parent.symbolic_max - self.rtkn
+        return sympy.Interval(left, right)
+
+    @memoized_meth
+    def _symbolic_thickness(self, **kwargs):
+        kwargs = {'dtype': np.int32, 'is_const': True, 'nonnegative': True}
+
+        names = ["%s_%stkn" % (self.parent.name, s) for s in ('l', 'r')]
+        return SubDimensionThickness(*[Symbol(name=n, **kwargs) for n in names])
 
     @cached_property
     def symbolic_min(self):
@@ -583,28 +660,20 @@ def symbolic_size(self):
     def thickness(self):
         return self._thickness
 
-    @cached_property
-    def _thickness_map(self):
-        return dict(self.thickness)
-
-    @property
-    def is_abstract(self):
-        return all(i is None for i in flatten(self.thickness))
+    tkns = thickness  # Shortcut for thickness
 
     @property
     def ltkn(self):
         # Shortcut for the left thickness symbol
-        return self.thickness.left[0]
+        return self.thickness.left
 
     @property
     def rtkn(self):
         # Shortcut for the right thickness symbol
-        return self.thickness.right[0]
+        return self.thickness.right
 
-    @property
-    def tkns(self):
-        # Shortcut for both thickness symbols
-        return self.ltkn, self.rtkn
+    def __hash__(self):
+        return id(self)
 
 
 class SubDimension(AbstractSubDimension):
@@ -659,37 +728,47 @@ class SubDimension(AbstractSubDimension):
 
     __rargs__ = AbstractSubDimension.__rargs__ + ('local',)
 
-    def __init_finalize__(self, name, parent, left, right, thickness, local,
+    _thickness_type = Thickness
+
+    def __init_finalize__(self, name, parent, thickness, local,
                           **kwargs):
-        super().__init_finalize__(name, parent, left, right, thickness)
         self._local = local
+        super().__init_finalize__(name, parent, thickness)
 
     @classmethod
     def left(cls, name, parent, thickness, local=True):
-        lst, rst = cls._symbolic_thickness(parent.name)
-        return cls(name, parent,
-                   left=parent.symbolic_min,
-                   right=parent.symbolic_min+lst-1,
-                   thickness=((lst, thickness), (rst, None)),
-                   local=local)
+        return cls(name, parent, thickness=(thickness, None), local=local)
 
     @classmethod
     def right(cls, name, parent, thickness, local=True):
-        lst, rst = cls._symbolic_thickness(parent.name)
-        return cls(name, parent,
-                   left=parent.symbolic_max-rst+1,
-                   right=parent.symbolic_max,
-                   thickness=((lst, None), (rst, thickness)),
-                   local=local)
+        return cls(name, parent, thickness=(None, thickness), local=local)
 
     @classmethod
     def middle(cls, name, parent, thickness_left, thickness_right, local=False):
-        lst, rst = cls._symbolic_thickness(parent.name)
-        return cls(name, parent,
-                   left=parent.symbolic_min+lst,
-                   right=parent.symbolic_max-rst,
-                   thickness=((lst, thickness_left), (rst, thickness_right)),
-                   local=local)
+        return cls(name, parent, thickness=(thickness_left, thickness_right), local=local)
+
+    @memoized_meth
+    def _symbolic_thickness(self, thickness=None):
+        kwargs = {'dtype': np.int32, 'is_const': True, 'nonnegative': True,
+                  'root': self.root, 'local': self.local}
+
+        names = ["%s_%stkn" % (self.parent.name, s) for s in ('l', 'r')]
+        sides = [LEFT, RIGHT]
+        return SubDimensionThickness(*[Thickness(name=n, side=s, value=t, **kwargs)
+                                       for n, s, t in zip(names, sides, thickness)])
+
+    @cached_property
+    def _interval(self):
+        if self.thickness.right.value is None:  # Left SubDimension
+            left = self.parent.symbolic_min
+            right = self.parent.symbolic_min + self.ltkn - 1
+        elif self.thickness.left.value is None:  # Right SubDimension
+            left = self.parent.symbolic_max - self.rtkn + 1
+            right = self.parent.symbolic_max
+        else:  # Middle SubDimension
+            return super()._interval
+
+        return sympy.Interval(left, right)
 
     @property
     def local(self):
@@ -697,11 +776,11 @@ def local(self):
 
     @property
     def is_left(self):
-        return self.thickness.right[1] is None
+        return self.thickness.right.value is None
 
     @property
     def is_right(self):
-        return self.thickness.left[1] is None
+        return self.thickness.left.value is None
 
     @property
     def is_middle(self):
@@ -716,89 +795,17 @@ def bound_symbols(self):
     def _maybe_distributed(self):
         return not self.local
 
-    @cached_property
-    def _offset_left(self):
-        # The left extreme of the SubDimension can be related to either the
-        # min or max of the parent dimension
-        try:
-            symbolic_thickness = self.symbolic_min - self.parent.symbolic_min
-            val = symbolic_thickness.subs(self._thickness_map)
-            return SubDimensionOffset(
-                int(val),
-                self.parent.symbolic_min,
-                symbolic_thickness
-            )
-        except TypeError:
-            symbolic_thickness = self.symbolic_min - self.parent.symbolic_max
-            val = symbolic_thickness.subs(self._thickness_map)
-            return SubDimensionOffset(
-                int(val),
-                self.parent.symbolic_max,
-                symbolic_thickness
-            )
-
-    @cached_property
-    def _offset_right(self):
-        # The right extreme of the SubDimension can be related to either the
-        # min or max of the parent dimension
-        try:
-            symbolic_thickness = self.symbolic_max - self.parent.symbolic_min
-            val = symbolic_thickness.subs(self._thickness_map)
-            return SubDimensionOffset(
-                int(val),
-                self.parent.symbolic_min,
-                symbolic_thickness
-            )
-        except TypeError:
-            symbolic_thickness = self.symbolic_max - self.parent.symbolic_max
-            val = symbolic_thickness.subs(self._thickness_map)
-            return SubDimensionOffset(
-                int(val),
-                self.parent.symbolic_max,
-                symbolic_thickness
-            )
-
     @property
     def _arg_names(self):
-        return tuple(k.name for k, _ in self.thickness) + self.parent._arg_names
+        return tuple(k.name for k in self.thickness) + self.parent._arg_names
 
     def _arg_defaults(self, grid=None, **kwargs):
         return {}
 
     def _arg_values(self, interval, grid=None, **kwargs):
-        # Allow override of thickness values to disable BCs
-        # However, arguments from the user are considered global
-        # So overriding the thickness to a nonzero value should not cause
-        # boundaries to exist between ranks where they did not before
-        r_ltkn, r_rtkn = (
-            kwargs.get(k.name, v) for k, v in self.thickness
-        )
-
-        if grid is not None and grid.is_distributed(self.root):
-            # Get local thickness
-            if self.local:
-                # dimension is of type ``left``/right`` - compute the 'offset'
-                # and then add 1 to get the appropriate thickness
-                if r_ltkn is not None:
-                    ltkn = grid.distributor.glb_to_loc(self.root, r_ltkn-1, LEFT)
-                    ltkn = ltkn+1 if ltkn is not None else 0
-                else:
-                    ltkn = 0
-
-                if r_rtkn is not None:
-                    rtkn = grid.distributor.glb_to_loc(self.root, r_rtkn-1, RIGHT)
-                    rtkn = rtkn+1 if rtkn is not None else 0
-                else:
-                    rtkn = 0
-            else:
-                # dimension is of type ``middle``
-                ltkn = grid.distributor.glb_to_loc(self.root, r_ltkn, LEFT) or 0
-                rtkn = grid.distributor.glb_to_loc(self.root, r_rtkn, RIGHT) or 0
-        else:
-            ltkn = r_ltkn or 0
-            rtkn = r_rtkn or 0
-
-        return {i.name: v for i, v in zip(self._thickness_map, (ltkn, rtkn))}
+        # SubDimension thicknesses at runtime are calculated by the thicknesses
+        # themselves
+        return {}
 
 
 class MultiSubDimension(AbstractSubDimension):
@@ -809,51 +816,16 @@ class MultiSubDimension(AbstractSubDimension):
 
     is_MultiSub = True
 
-    __rargs__ = (DerivedDimension.__rargs__ + ('thickness',))
     __rkwargs__ = ('functions', 'bounds_indices', 'implicit_dimension')
 
     def __init_finalize__(self, name, parent, thickness, functions=None,
                           bounds_indices=None, implicit_dimension=None):
-        # Canonicalize thickness
-        if thickness is None:
-            # Using dummy left/right is the only thing we can do for such
-            # an abstract MultiSubDimension
-            thickness = ((None, None), (None, None))
-
-            left = sympy.S.NegativeInfinity
-            right = sympy.S.Infinity
-        elif isinstance(thickness, tuple):
-            if all(isinstance(i, Symbol) for i in thickness):
-                ltkn, rtkn = thickness
-                thickness = ((ltkn, None), (rtkn, None))
-            elif all(isinstance(i, tuple) and len(i) == 2 for i in thickness):
-                (ltkn, _), (rtkn, _) = thickness
 
-            try:
-                left = parent.symbolic_min + ltkn
-                right = parent.symbolic_max - rtkn
-            except TypeError:
-                # May end up here after a reconstruction
-                left = sympy.S.NegativeInfinity
-                right = sympy.S.Infinity
-        else:
-            raise ValueError("MultiSubDimension expects a tuple of thicknesses")
-
-        super().__init_finalize__(name, parent, left, right, thickness)
+        super().__init_finalize__(name, parent, thickness)
         self.functions = functions
         self.bounds_indices = bounds_indices
         self.implicit_dimension = implicit_dimension
 
-    def __hash__(self):
-        # There is no possibility for two MultiSubDimensions to ever hash the
-        # same, since a MultiSubDimension carries a reference to a MultiSubDomain,
-        # which is unique
-        return id(self)
-
-    @classmethod
-    def _symbolic_thickness(cls, name):
-        return super()._symbolic_thickness(name, stype=Symbol)
-
     @cached_property
     def bound_symbols(self):
         return self.parent.bound_symbols
diff --git a/examples/cfd/01_convection_revisited.ipynb b/examples/cfd/01_convection_revisited.ipynb
index 0d5d405618..259a47744d 100644
--- a/examples/cfd/01_convection_revisited.ipynb
+++ b/examples/cfd/01_convection_revisited.ipynb
@@ -426,9 +426,9 @@
       "  for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))\n",
       "  {\n",
       "    START(section0)\n",
-      "    for (int x = x_ltkn0 + x_m; x <= x_M - x_rtkn0; x += 1)\n",
+      "    for (int x = x_m + x_ltkn0; x <= x_M - x_rtkn0; x += 1)\n",
       "    {\n",
-      "      for (int y = y_ltkn0 + y_m; y <= y_M - y_rtkn0; y += 1)\n",
+      "      for (int y = y_m + y_ltkn0; y <= y_M - y_rtkn0; y += 1)\n",
       "      {\n",
       "        u[t1][x + 1][y + 1] = dt*(-(-u[t0][x][y + 1]/h_x + u[t0][x + 1][y + 1]/h_x) - (-u[t0][x + 1][y]/h_y + u[t0][x + 1][y + 1]/h_y) + u[t0][x + 1][y + 1]/dt);\n",
       "      }\n",
diff --git a/tests/test_caching.py b/tests/test_caching.py
index 8dca69fa60..93200d3d73 100644
--- a/tests/test_caching.py
+++ b/tests/test_caching.py
@@ -663,8 +663,8 @@ def test_sparse_function(self, operate_on_empty_cache):
         i = u.inject(expr=u, field=u)
 
         # created: rux, ruy (radius dimensions) and spacings
-        # posx, posy, px, py, u_coords (as indexified),
-        ncreated = 2+1+2+2+2+1
+        # posx, posy, px, py, u_coords (as indexified), x_m, x_M, y_m, y_M
+        ncreated = 2+1+2+2+2+1+4
         # Note that injection is now lazy so no new symbols should be created
         assert len(_SymbolCache) == cur_cache_size
         i.evaluate
@@ -684,14 +684,21 @@ def test_sparse_function(self, operate_on_empty_cache):
         # in the first clear_cache they were still referenced by their "parent" objects
         # (e.g., ru* by ConditionalDimensions, through `condition`)
 
-        assert len(_SymbolCache) == init_cache_size + 8
+        assert len(_SymbolCache) == init_cache_size + 12
         clear_cache()
         # Now we should be back to the original state except for
-        # pos* that belong to the abstract class
-        assert len(_SymbolCache) == init_cache_size + 2
+        # pos* that belong to the abstract class and the dimension
+        # bounds (x_m, x_M, y_m, y_M)
+        assert len(_SymbolCache) == init_cache_size + 6
         clear_cache()
-        # Now we should be back to the original state
-        assert len(_SymbolCache) == init_cache_size
+        # Now we should be back to the original state plus the dimension bounds
+        # (x_m, x_M, y_m, y_M)
+        assert len(_SymbolCache) == init_cache_size + 4
+        # Delete the grid and check that all symbols are subsequently garbage collected
+        del grid
+        for n in (10, 3, 0):
+            clear_cache()
+            assert len(_SymbolCache) == n
 
     def test_after_indexification(self):
         """
@@ -842,10 +849,10 @@ def test_solve(self, operate_on_empty_cache):
         # created by the finite difference (u.dt, u.dx2). We would have had
         # three extra references to u(t + dt), u(x - h_x) and u(x + h_x).
         # But this is not the case anymore!
-        assert len(_SymbolCache) == 12
-        clear_cache()
         assert len(_SymbolCache) == 8
         clear_cache()
+        assert len(_SymbolCache) == 4
+        clear_cache()
         assert len(_SymbolCache) == 2
         clear_cache()
         assert len(_SymbolCache) == 0
diff --git a/tests/test_dimension.py b/tests/test_dimension.py
index 6e63159cdb..84df27ecef 100644
--- a/tests/test_dimension.py
+++ b/tests/test_dimension.py
@@ -323,15 +323,15 @@ def test_symbolic_size(self):
         xleft = SubDimension.left(name='xleft', parent=x, thickness=thickness)
         assert xleft.is_left
         assert not xleft.is_middle
-        assert xleft.symbolic_size == xleft.thickness.left[0]
+        assert xleft.symbolic_size == xleft.thickness.left
 
         xi = SubDimension.middle(name='xi', parent=x,
                                  thickness_left=thickness, thickness_right=thickness)
         assert xi.symbolic_size == (x.symbolic_max - x.symbolic_min -
-                                    xi.thickness.left[0] - xi.thickness.right[0] + 1)
+                                    xi.thickness.left - xi.thickness.right + 1)
 
         xright = SubDimension.right(name='xright', parent=x, thickness=thickness)
-        assert xright.symbolic_size == xright.thickness.right[0]
+        assert xright.symbolic_size == xright.thickness.right
 
     @pytest.mark.parametrize('opt', opts_tiling)
     def test_bcs(self, opt):
@@ -787,6 +787,18 @@ def test_expandingbox_like(self, opt):
         assert np.all(u.data[:, :, :2] == 0.)
         assert np.all(u.data[:, :, -2:] == 0.)
 
+    def test_standalone_thickness(self):
+        x = Dimension('x')
+        ix = SubDimension.left('ix', x, 2)
+        f = Function(name="f", dimensions=(ix,), shape=(5,), dtype=np.int32)
+
+        eqns = Eq(f[ix.symbolic_max], 1)
+
+        op = Operator(eqns)
+        assert 'x_ltkn0' in str(op.ccode)
+        op(x_m=0)
+        assert np.all(f.data == np.array([0, 1, 0, 0, 0]))
+
 
 class TestConditionalDimension:
 
diff --git a/tests/test_dle.py b/tests/test_dle.py
index 45583d4346..7f1ae73186 100644
--- a/tests/test_dle.py
+++ b/tests/test_dle.py
@@ -146,7 +146,7 @@ def test_cache_blocking_structure_subdims():
     # zi is rebuilt with name z, so check symbolic max and min are preserved
     # Also check the zi was rebuilt
     assert not tree[4].dim.is_Block and tree[4].dim is not zi and\
-        str(tree[4].dim.symbolic_min) == 'z_ltkn0 + z_m' and\
+        str(tree[4].dim.symbolic_min) == 'z_m + z_ltkn0' and\
         str(tree[4].dim.symbolic_max) == 'z_M - z_rtkn0' and\
         tree[4].dim.parent is z
 
@@ -1374,7 +1374,7 @@ def test_nested_cache_blocking_structure_subdims(self, blocklevels):
 
         if blocklevels == 1:
             assert not tree[4].dim.is_Block and tree[4].dim is not zi and\
-                str(tree[4].dim.symbolic_min) == 'z_ltkn0 + z_m' and\
+                str(tree[4].dim.symbolic_min) == 'z_m + z_ltkn0' and\
                 str(tree[4].dim.symbolic_max) == 'z_M - z_rtkn0' and\
                 tree[4].dim.parent is z
         elif blocklevels == 2:
@@ -1385,7 +1385,7 @@ def test_nested_cache_blocking_structure_subdims(self, blocklevels):
             assert tree[5].dim.is_Block and tree[5].dim.parent is tree[3].dim and\
                 tree[5].dim.root is y
             assert not tree[6].dim.is_Block and tree[6].dim is not zi and\
-                str(tree[6].dim.symbolic_min) == 'z_ltkn0 + z_m' and\
+                str(tree[6].dim.symbolic_min) == 'z_m + z_ltkn0' and\
                 str(tree[6].dim.symbolic_max) == 'z_M - z_rtkn0' and\
                 tree[6].dim.parent is z
 
diff --git a/tests/test_subdomains.py b/tests/test_subdomains.py
index d95fec84bb..784dc1c3ad 100644
--- a/tests/test_subdomains.py
+++ b/tests/test_subdomains.py
@@ -197,7 +197,7 @@ def define(self, dimensions):
         check = np.zeros(grid.shape)
 
         mM_map = {x.symbolic_min: 0, x.symbolic_max: grid.shape[0]-1}
-        t_map = {k: v for k, v in xd0._thickness_map.items() if v is not None}
+        t_map = {tkn: tkn.value for tkn in xd0.thickness if tkn.value is not None}
         start = int(xd0.symbolic_min.subs({**mM_map, **t_map}))
         stop = int(xd0.symbolic_max.subs({**mM_map, **t_map})+1)
 
diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py
index 353fdc934c..496e325387 100644
--- a/tests/test_symbolics.py
+++ b/tests/test_symbolics.py
@@ -130,15 +130,15 @@ def test_subdimension():
 
     di = SubDimension.middle(name='di', parent=d, thickness_left=4, thickness_right=4)
     assert di.free_symbols == {di}
-    assert di.bound_symbols == {d.symbolic_min, d.symbolic_max} | set(di._thickness_map)
+    assert di.bound_symbols == {d.symbolic_min, d.symbolic_max} | set(di.thickness)
 
     dl = SubDimension.left(name='dl', parent=d, thickness=4)
     assert dl.free_symbols == {dl}
-    assert dl.bound_symbols == {d.symbolic_min, dl.thickness.left[0]}
+    assert dl.bound_symbols == {d.symbolic_min, dl.thickness.left}
 
     dr = SubDimension.right(name='dr', parent=d, thickness=4)
     assert dr.free_symbols == {dr}
-    assert dr.bound_symbols == {d.symbolic_max, dr.thickness.right[0]}
+    assert dr.bound_symbols == {d.symbolic_max, dr.thickness.right}
 
 
 def test_timefunction():