Repository: devitocodes/devito
Issue: Thickness overriding when using a different grid is broken

Explanation. Currently:

* `Grid` creates `Dimensions`; each `Dimension` is provided with a `Constant` spacing. Therefore, the physical value is attached to `Dimension` via the spacing.
* Two grids created via `Grid(shape=(...))` are defined over different `Dimension`s .
* However, this means that `grid.is_distributed(self.root)` in `SubDimension._arg_values(...)` will return false when `grid` is not the Grid used to construct the `Operator`. This breaks the thickness overriding.
* One could then say "OK, I'm then going to create the second grid via `Grid(shape=(...), dimensions=grid1.dimensions)`", but the issue here is that then both grids will **share the same state** (in particular, the Dimension spacing).

Reproducer:

```
from devito import *
import numpy as np

grid = Grid(shape=(10, 10))
x, y = grid.dimensions

xi = SubDimension.middle(name='xi', parent=x, thickness_left=3, thickness_right=3)
yi = SubDimension.middle(name='yi', parent=y, thickness_left=3, thickness_right=3)

u = TimeFunction(name='u', grid=grid, space_order=2, time_order=0)

eqn = Eq(u.forward, u + 1).subs({x: xi, y: yi})

op = Operator(eqn)

grid2 = Grid(shape=(10, 10), dimensions=(x, y)) 

u2 = TimeFunction(name='u', grid=grid2, space_order=2, time_order=0)

op.apply(time_M=0, u=u2)

assert np.all(u2.data[0, 3:-3, 3:-3] == 1.) 
```

Changed files:
tests/test_derivatives.py, devito/operator/operator.py, devito/types/sparse.py, devito/types/dense.py, tests/test_caching.py, tests/test_pickle.py, tests/test_dimension.py, devito/types/dimension.py, tests/test_interpolation.py, devito/types/basic.py, tests/test_mpi.py, devito/types/grid.py, tests/test_autotuner.py

diff --git a/devito/operator/operator.py b/devito/operator/operator.py
index 570702b01..a14340e0c 100644
--- a/devito/operator/operator.py
+++ b/devito/operator/operator.py
@@ -456,6 +456,7 @@ class Operator(Callable):
             raise ValueError("Multiple Grids found")
         try:
             grid = grids.pop()
+            args.update(grid._arg_values(**kwargs))
         except KeyError:
             grid = None
 
@@ -517,8 +518,16 @@ class Operator(Callable):
     @cached_property
     def _known_arguments(self):
         """The arguments that can be passed to ``apply`` when running the Operator."""
-        ret = set.union(*[set(i._arg_names) for i in self.input + self.dimensions])
-        return tuple(sorted(ret))
+        ret = set()
+        for i in self.input:
+            ret.update(i._arg_names)
+            try:
+                ret.update(i.grid._arg_names)
+            except AttributeError:
+                pass
+        for d in self.dimensions:
+            ret.update(d._arg_names)
+        return frozenset(ret)
 
     def _autotune(self, args, setup):
         """Auto-tuning to improve runtime performance."""
diff --git a/devito/types/basic.py b/devito/types/basic.py
index 5ac0d061b..205860362 100644
--- a/devito/types/basic.py
+++ b/devito/types/basic.py
@@ -466,6 +466,15 @@ class Scalar(Symbol, ArgProvider):
     def __dtype_setup__(cls, **kwargs):
         return kwargs.get('dtype', np.float32)
 
+    def _arg_defaults(self):
+        return {}
+
+    def _arg_values(self, **kwargs):
+        if self.name in kwargs:
+            return {self.name: kwargs.pop(self.name)}
+        else:
+            return {}
+
 
 class AbstractTensor(sympy.ImmutableDenseMatrix, Basic, Pickable, Evaluable):
 
diff --git a/devito/types/dense.py b/devito/types/dense.py
index 0d3d22536..ba3095f17 100644
--- a/devito/types/dense.py
+++ b/devito/types/dense.py
@@ -810,10 +810,6 @@ class DiscreteFunction(AbstractFunction, ArgProvider, Differentiable):
         for i, s in zip(key.dimensions, self.shape):
             args.update(i._arg_defaults(_min=0, size=s))
 
-        # Add value overrides associated with the Grid
-        if self.grid is not None:
-            args.update(self.grid._arg_defaults())
-
         return args
 
     def _arg_values(self, **kwargs):
@@ -840,9 +836,6 @@ class DiscreteFunction(AbstractFunction, ArgProvider, Differentiable):
                 for i, s in zip(self.dimensions, new.shape):
                     size = s - sum(self._size_nodomain[i])
                     values.update(i._arg_defaults(size=size))
-                # Add value overrides associated with the Grid
-                if self.grid is not None:
-                    values.update(self.grid._arg_defaults())
         else:
             values = self._arg_defaults(alias=self).reduce_all()
 
diff --git a/devito/types/dimension.py b/devito/types/dimension.py
index c4e096cbd..d2d2d599d 100644
--- a/devito/types/dimension.py
+++ b/devito/types/dimension.py
@@ -6,7 +6,7 @@ from cached_property import cached_property
 
 from devito.data import LEFT, RIGHT
 from devito.exceptions import InvalidArgument
-from devito.logger import debug, warning
+from devito.logger import debug
 from devito.tools import Pickable, dtype_to_cstr, is_integer, memoized_meth
 from devito.types.args import ArgProvider
 from devito.types.basic import Symbol, DataSymbol, Scalar
@@ -121,6 +121,15 @@ class Dimension(ArgProvider):
         else:
             return BasicDimension.__new__(cls, *args, **kwargs)
 
+    @classmethod
+    def class_key(cls):
+        """
+        Overrides sympy.Symbol.class_key such that Dimensions always
+        preceed other symbols when printed (e.g. x + h_x, not h_x + x).
+        """
+        a, b, c = super().class_key()
+        return a, b - 1, c
+
     @classmethod
     def __dtype_setup__(cls, **kwargs):
         # Unlike other Symbols, Dimensions can only be integers
@@ -261,20 +270,7 @@ class Dimension(ArgProvider):
             except (AttributeError, TypeError):
                 pass
 
-        args = {self.min_name: loc_minv, self.max_name: loc_maxv}
-
-        # Maybe override spacing
-        if grid is not None:
-            try:
-                spacing_map = {k.name: v for k, v in grid.spacing_map.items()}
-                args[self.spacing.name] = spacing_map[self.spacing.name]
-            except KeyError:
-                pass
-            except AttributeError:
-                # See issue #1524
-                warning("Unable to override spacing")
-
-        return args
+        return {self.min_name: loc_minv, self.max_name: loc_maxv}
 
     def _arg_check(self, args, size, interval):
         """
@@ -1204,7 +1200,7 @@ class CustomDimension(BasicDimension):
         if self.is_Derived:
             return self.parent.spacing
         else:
-            return self.spacing
+            return self._spacing
 
     @cached_property
     def _defines(self):
diff --git a/devito/types/grid.py b/devito/types/grid.py
index 9ff0853c6..d97d2f7f5 100644
--- a/devito/types/grid.py
+++ b/devito/types/grid.py
@@ -1,14 +1,15 @@
 from collections import namedtuple
+from math import floor
 
 import numpy as np
 from sympy import prod
-from math import floor
+from cached_property import cached_property
 
 from devito.data import LEFT, RIGHT
 from devito.mpi import Distributor, MPI
 from devito.tools import ReducerMap, as_tuple, memoized_meth
 from devito.types.args import ArgProvider
-from devito.types.constant import Constant
+from devito.types.basic import Scalar
 from devito.types.dense import Function
 from devito.types.dimension import (Dimension, SpaceDimension, TimeDimension,
                                     SteppingDimension, SubDimension)
@@ -111,11 +112,19 @@ class Grid(ArgProvider):
             # Create the spatial dimensions and constant spacing symbols
             assert(self.dim <= 3)
             dim_names = self._default_dimensions[:self.dim]
-            dim_spacing = tuple(self._const(name='h_%s' % n, value=v, dtype=self.dtype)
-                                for n, v in zip(dim_names, self.spacing))
+            dim_spacing = tuple(Scalar(name='h_%s' % n, dtype=self.dtype, is_const=True)
+                                for n in dim_names)
             self._dimensions = tuple(SpaceDimension(name=n, spacing=s)
                                      for n, s in zip(dim_names, dim_spacing))
         else:
+            # Sanity check
+            for d in dimensions:
+                if not d.is_Space:
+                    raise ValueError("Cannot create Grid with Dimension `%s` "
+                                     "since it's not a SpaceDimension" % d)
+                if d.is_Derived and not d.is_Conditional:
+                    raise ValueError("Cannot create Grid with derived Dimension `%s` "
+                                     "of type `%s`" % (d, type(d)))
             self._dimensions = dimensions
 
         self._distributor = Distributor(self.shape, self.dimensions, comm, topology)
@@ -127,21 +136,23 @@ class Grid(ArgProvider):
                                      distributor=self._distributor, counter=counter)
         self._subdomains = subdomains
 
-        origin = as_tuple(origin or tuple(0. for _ in self.shape))
-        self._origin = tuple(self._const(name='o_%s' % d.name, value=v, dtype=self.dtype)
-                             for d, v in zip(self.dimensions, origin))
+        self._origin = as_tuple(origin or tuple(0. for _ in self.shape))
+        self._origin_symbols = tuple(Scalar(name='o_%s' % d.name, dtype=dtype,
+                                            is_const=True)
+                                     for d in self.dimensions)
 
         # Sanity check
         assert (self.dim == len(self.origin) == len(self.extent) == len(self.spacing))
 
         # Store or create default symbols for time and stepping dimensions
         if time_dimension is None:
-            spacing = self._const(name='dt', dtype=self.dtype)
-            self._time_dim = self._make_time_dim(spacing)
-            self._stepping_dim = self._make_stepping_dim(self.time_dim, name='t')
+            spacing = Scalar(name='dt', dtype=dtype, is_const=True)
+            self._time_dim = TimeDimension(name='time', spacing=spacing)
+            self._stepping_dim = SteppingDimension(name='t', parent=self.time_dim)
         elif isinstance(time_dimension, TimeDimension):
             self._time_dim = time_dimension
-            self._stepping_dim = self._make_stepping_dim(self.time_dim)
+            self._stepping_dim = SteppingDimension(name='%s_s' % self.time_dim.name,
+                                                   parent=self.time_dim)
         else:
             raise ValueError("`time_dimension` must be None or of type TimeDimension")
 
@@ -165,10 +176,22 @@ class Grid(ArgProvider):
         """Physical coordinates of the domain origin."""
         return self._origin
 
+    @property
+    def origin_symbols(self):
+        """Symbols representing the grid origin in each SpaceDimension."""
+        return self._origin_symbols
+
     @property
     def origin_map(self):
-        """Map between origin symbols and their values"""
-        return {o: o.data for o in self.origin}
+        """Map between origin symbols and their values."""
+        return dict(zip(self.origin_symbols, self.origin))
+
+    @property
+    def origin_offset(self):
+        """Offset of the local (per-process) origin from the domain origin."""
+        grid_origin = [min(i) for i in self.distributor.glb_numb]
+        assert len(grid_origin) == len(self.spacing)
+        return tuple(i*h for i, h in zip(grid_origin, self.spacing))
 
     @property
     def dimensions(self):
@@ -211,22 +234,29 @@ class Grid(ArgProvider):
         spacing = (np.array(self.extent) / (np.array(self.shape) - 1)).astype(self.dtype)
         return as_tuple(spacing)
 
-    @property
+    @cached_property
     def spacing_symbols(self):
-        """Symbols representing the grid spacing in each SpaceDimension"""
-        return as_tuple(d.spacing for d in self.dimensions)
+        """Symbols representing the grid spacing in each SpaceDimension."""
+        return as_tuple(d.root.spacing for d in self.dimensions)
 
-    @property
+    @cached_property
     def spacing_map(self):
         """Map between spacing symbols and their values for each SpaceDimension."""
-        return dict(zip(self.spacing_symbols, self.spacing))
+        mapper = {}
+        for d, s in zip(self.dimensions, self.spacing):
+            if d.is_Conditional:
+                # Special case subsampling: `Grid.dimensions` -> (xb, yb, zb)`
+                # where `xb, yb, zb` are ConditionalDimensions whose parents
+                # are SpaceDimensions
+                mapper[d.root.spacing] = s/self.dtype(d.factor)
+            elif d.is_Space:
+                # Typical case: `Grid.dimensions` -> (x, y, z)` where `x, y, z` are
+                # the SpaceDimensions
+                mapper[d.spacing] = s
+            else:
+                assert False
 
-    @property
-    def origin_offset(self):
-        """Offset of the local (per-process) origin from the domain origin."""
-        grid_origin = [min(i) for i in self.distributor.glb_numb]
-        assert len(grid_origin) == len(self.spacing)
-        return tuple(i*h for i, h in zip(grid_origin, self.spacing))
+        return mapper
 
     @property
     def shape(self):
@@ -258,29 +288,35 @@ class Grid(ArgProvider):
         """True if ``dim`` is a distributed Dimension, False otherwise."""
         return any(dim is d for d in self.distributor.dimensions)
 
-    @property
-    def _const(self):
-        """The type to be used to create constant symbols."""
-        return Constant
-
-    def _make_stepping_dim(self, time_dim, name=None):
-        """Create a SteppingDimension for this Grid."""
-        if name is None:
-            name = '%s_s' % time_dim.name
-        return SteppingDimension(name=name, parent=time_dim)
-
-    def _make_time_dim(self, spacing):
-        """Create a TimeDimension for this Grid."""
-        return TimeDimension(name='time', spacing=spacing)
+    @cached_property
+    def _arg_names(self):
+        ret = []
+        ret.append(self.time_dim.spacing.name)
+        ret.extend([i.name for i in self.origin_map])
+        for i in self.spacing_map:
+            try:
+                ret.append(i.name)
+            except AttributeError:
+                # E.g., {n*h_x: v} (the case of ConditionalDimension)
+                ret.extend([a.name for a in i.free_symbols])
+        return tuple(ret)
 
     @memoized_meth
     def _arg_defaults(self):
         """A map of default argument values defined by this Grid."""
         args = ReducerMap()
 
+        # Dimensions size
         for k, v in self.dimension_map.items():
             args.update(k._arg_defaults(_min=0, size=v.loc))
 
+        # Dimensions spacing
+        args.update({k.name: v for k, v in self.spacing_map.items()})
+
+        # Grid origin
+        args.update({k.name: v for k, v in self.origin_map.items()})
+
+        # MPI-related objects
         if self.distributor.is_parallel:
             distributor = self.distributor
             args[distributor._obj_comm.name] = distributor._obj_comm.value
@@ -288,6 +324,14 @@ class Grid(ArgProvider):
 
         return args
 
+    def _arg_values(self, **kwargs):
+        values = dict(self._arg_defaults())
+
+        # Override spacing and origin if necessary
+        values.update({i: kwargs[i] for i in self._arg_names if i in kwargs})
+
+        return values
+
     def __getstate__(self):
         state = self.__dict__.copy()
         # A Distributor wraps an MPI communicator, which can't and shouldn't be pickled
diff --git a/devito/types/sparse.py b/devito/types/sparse.py
index 20d102d76..a94772c71 100644
--- a/devito/types/sparse.py
+++ b/devito/types/sparse.py
@@ -293,9 +293,6 @@ class AbstractSparseFunction(DiscreteFunction):
             for i, s in zip(mapper[k].indices, v.shape):
                 args.update(i._arg_defaults(_min=0, size=s))
 
-        # Add MPI-related data structures
-        args.update(self.grid._arg_defaults())
-
         return args
 
     def _eval_at(self, func):
@@ -317,8 +314,6 @@ class AbstractSparseFunction(DiscreteFunction):
                     for i, s in zip(k.indices, v.shape):
                         size = s - sum(k._size_nodomain[i])
                         values.update(i._arg_defaults(size=size))
-                # Add value overrides associated with the Grid
-                values.update(self.grid._arg_defaults())
         else:
             values = self._arg_defaults(alias=self).reduce_all()
 
@@ -567,8 +562,9 @@ class SparseFunction(AbstractSparseFunction):
         """
         symbols = [Scalar(name='pos%s' % d, dtype=self.dtype)
                    for d in self.grid.dimensions]
-        return OrderedDict([(c - o, p) for p, c, o in
-                            zip(symbols, self._coordinate_symbols, self.grid.origin)])
+        return OrderedDict([(c - o, p) for p, c, o in zip(symbols,
+                                                          self._coordinate_symbols,
+                                                          self.grid.origin_symbols)])
 
     @cached_property
     def _point_increments(self):
@@ -587,14 +583,14 @@ class SparseFunction(AbstractSparseFunction):
         """Symbol for each grid index according to the coordinates."""
         return tuple([INT(FLOOR((c - o) / i.spacing))
                       for c, o, i in zip(self._coordinate_symbols,
-                                         self.grid.origin,
+                                         self.grid.origin_symbols,
                                          self.grid.dimensions[:self.grid.dim])])
 
     def _coordinate_bases(self, field_offset):
         """Symbol for the base coordinates of the reference grid point."""
         return tuple([cast_mapper[self.dtype](c - o - idx * i.spacing)
                       for c, o, idx, i, of in zip(self._coordinate_symbols,
-                                                  self.grid.origin,
+                                                  self.grid.origin_symbols,
                                                   self._coordinate_indices,
                                                   self.grid.dimensions[:self.grid.dim],
                                                   field_offset)])
@@ -624,8 +620,8 @@ class SparseFunction(AbstractSparseFunction):
             raise ValueError("No coordinates attached to this SparseFunction")
         ret = []
         for coords in self.coordinates.data._local:
-            ret.append(tuple(int(np.floor(c - o.data)/i.spacing.data) for c, o, i in
-                             zip(coords, self.grid.origin, self.grid.dimensions)))
+            ret.append(tuple(int(np.floor(c - o)/s) for c, o, s in
+                             zip(coords, self.grid.origin, self.grid.spacing)))
         return tuple(ret)
 
     def guard(self, expr=None, offset=0):
diff --git a/tests/test_autotuner.py b/tests/test_autotuner.py
index c5b94272e..d510dab24 100644
--- a/tests/test_autotuner.py
+++ b/tests/test_autotuner.py
@@ -127,7 +127,7 @@ def test_tti_aggressive():
     from test_dse import TestTTI
     wave_solver = TestTTI().tti_operator(opt='advanced')
     op = wave_solver.op_fwd(kernel='centered')
-    op.apply(time=0, autotune='aggressive')
+    op.apply(time=0, autotune='aggressive', dt=0.1)
     assert op._state['autotuning'][0]['runs'] == 30
 
 
diff --git a/tests/test_caching.py b/tests/test_caching.py
index 778a5f041..00a0aabc2 100644
--- a/tests/test_caching.py
+++ b/tests/test_caching.py
@@ -358,8 +358,9 @@ class TestCaching(object):
 
     def test_grid_objs(self):
         """
-        Test that two different Grids use different Symbols/Dimensions. This is
-        because objects such as spacing and origin are Constants carrying a value.
+        Test that two different Grids use the same Symbols/Dimensions if possible
+        (i.e., if already in cache). This is because objects such as spacing and origin
+        are Scalars, which carry no value.
         """
         grid0 = Grid(shape=(4, 4))
         x0, y0 = grid0.dimensions
@@ -369,12 +370,12 @@ class TestCaching(object):
         x1, y1 = grid1.dimensions
         ox1, oy1 = grid1.origin
 
-        assert x0 is not x1
-        assert y0 is not y1
-        assert x0.spacing is not x1.spacing
-        assert y0.spacing is not y1.spacing
-        assert ox0 is not ox1
-        assert oy0 is not oy1
+        assert x0 is x1
+        assert y0 is y1
+        assert x0.spacing is x1.spacing
+        assert y0.spacing is y1.spacing
+        assert ox0 is ox1
+        assert oy0 is oy1
 
     def test_symbol_aliasing(self):
         """Test to assert that our aliasing cache isn't defeated by sympys
diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py
index a0a3ebcbe..e252fafdc 100644
--- a/tests/test_derivatives.py
+++ b/tests/test_derivatives.py
@@ -384,7 +384,7 @@ class TestFD(object):
         # Creates subsampled spatial dimensions and according grid
         dims = tuple([ConditionalDimension(d.name+'sub', parent=d, factor=2)
                       for d in u.grid.dimensions])
-        grid2 = Grid((6, 6), dimensions=dims)
+        grid2 = Grid((6, 6), dimensions=dims, extent=(10, 10))
         u2 = TimeFunction(name='u2', grid=grid2, save=nt, space_order=1)
         for i in range(nt):
             for j in range(u2.data_with_halo.shape[2]):
diff --git a/tests/test_dimension.py b/tests/test_dimension.py
index 76bb76153..0f734ccfd 100644
--- a/tests/test_dimension.py
+++ b/tests/test_dimension.py
@@ -589,7 +589,7 @@ class TestConditionalDimension(object):
         u2[x, y] = u[2*x, 2*y]
         """
         nt = 19
-        grid = Grid(shape=(12, 12))
+        grid = Grid(shape=(11, 11))
         time = grid.time_dim
 
         u = TimeFunction(name='u', grid=grid, save=nt)
@@ -606,7 +606,7 @@ class TestConditionalDimension(object):
         op = Operator(eqns)
         op.apply(time_M=nt-2)
         # Verify that u2[x,y]= u[2*x, 2*y]
-        assert np.allclose(u.data[:-1, 0:-1:2, 0:-1:2], u2.data[:-1, :, :])
+        assert np.allclose(u.data[:-1, 0::2, 0::2], u2.data[:-1, :, :])
 
     def test_time_subsampling_fd(self):
         nt = 19
@@ -640,7 +640,7 @@ class TestConditionalDimension(object):
         """
         Test that the FD shortcuts are handled correctly with ConditionalDimensions
         """
-        grid = Grid(shape=(21, 21))
+        grid = Grid(shape=(11, 11))
         time = grid.time_dim
         # Creates subsampled spatial dimensions and accordine grid
         dims = tuple([ConditionalDimension(d.name+'sub', parent=d, factor=2)
@@ -654,13 +654,13 @@ class TestConditionalDimension(object):
         # Verify that u2 contains subsampled fd values
         assert np.all(u2.data[0, :, :] == 2.)
         assert np.all(u2.data[1, 0, 0] == 0.)
-        assert np.all(u2.data[1, -1, -1] == -40.)
-        assert np.all(u2.data[1, 0, -1] == -20.)
-        assert np.all(u2.data[1, -1, 0] == -20.)
+        assert np.all(u2.data[1, -1, -1] == -20.)
+        assert np.all(u2.data[1, 0, -1] == -10.)
+        assert np.all(u2.data[1, -1, 0] == -10.)
         assert np.all(u2.data[1, 1:-1, 0] == 0.)
         assert np.all(u2.data[1, 0, 1:-1] == 0.)
-        assert np.all(u2.data[1, 1:-1, -1] == -20.)
-        assert np.all(u2.data[1, -1, 1:-1] == -20.)
+        assert np.all(u2.data[1, 1:-1, -1] == -10.)
+        assert np.all(u2.data[1, -1, 1:-1] == -10.)
         assert np.all(u2.data[1, 1:4, 1:4] == 0.)
 
     # This test generates an openmp loop form which makes older gccs upset
diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py
index 9c166b846..6c96eca02 100644
--- a/tests/test_interpolation.py
+++ b/tests/test_interpolation.py
@@ -468,7 +468,7 @@ def test_position(shape):
                                               num=nrec)
     rec2.coordinates.data[:, 1:] = src.coordinates.data[0, 1:]
 
-    ox_g, oy_g, oz_g = tuple(o.dtype(o.data+100.) for o in model.grid.origin)
+    ox_g, oy_g, oz_g = tuple(o + 100. for o in model.grid.origin)
 
     rec1, u1, _ = solver.forward(save=False, src=src, rec=rec2,
                                  o_x=ox_g, o_y=oy_g, o_z=oz_g)
diff --git a/tests/test_mpi.py b/tests/test_mpi.py
index f75a704a8..0b0ba42cc 100644
--- a/tests/test_mpi.py
+++ b/tests/test_mpi.py
@@ -2067,6 +2067,32 @@ class TestOperatorAdvanced(object):
         assert np.isclose(norm(u), 23.70654, atol=1e-5, rtol=0)
         assert np.isclose(norm(v), 21.14994, atol=1e-5, rtol=0)
 
+    @pytest.mark.parallel(mode=2)
+    def test_overriding_from_different_grid(self):
+        """
+        MFE for issue #1629.
+        """
+        grid = Grid(shape=(10, 10))
+        x, y = grid.dimensions
+        xi = SubDimension.middle(name='xi', parent=x, thickness_left=3, thickness_right=3)
+        yi = SubDimension.middle(name='yi', parent=y, thickness_left=3, thickness_right=3)
+        u = TimeFunction(name='u', grid=grid, space_order=2, time_order=0)
+
+        eqn = Eq(u.forward, u + 1).subs({x: xi, y: yi})
+        op = Operator(eqn)
+
+        grid2 = Grid(shape=(10, 10), dimensions=(x, y))
+        u2 = TimeFunction(name='u', grid=grid2, space_order=2, time_order=0)
+
+        op.apply(time_M=0, u=u2)
+        assert np.all(u2.data[0, 3:-3, 3:-3] == 1.)
+
+        grid3 = Grid(shape=(10, 10))
+        u3 = TimeFunction(name='u', grid=grid3, space_order=2, time_order=0)
+
+        op.apply(time_M=0, u=u3)
+        assert np.all(u3.data[0, 3:-3, 3:-3] == 1.)
+
 
 def gen_serial_norms(shape, so):
     """
diff --git a/tests/test_pickle.py b/tests/test_pickle.py
index 487ea640a..4dc61b520 100644
--- a/tests/test_pickle.py
+++ b/tests/test_pickle.py
@@ -625,10 +625,10 @@ def test_full_model():
                       np.linalg.norm(new_time_range.time_values))
 
     # Test Class Constant pickling
-    pkl_origin = pickle.dumps(model.grid.origin)
+    pkl_origin = pickle.dumps(model.grid.origin_symbols)
     new_origin = pickle.loads(pkl_origin)
 
-    for a, b in zip(model.grid.origin, new_origin):
+    for a, b in zip(model.grid.origin_symbols, new_origin):
         assert a.compare(b) == 0
 
     # Test Class TimeDimension pickling