diff --git a/devito/__init__.py b/devito/__init__.py
index 01279350a2..f435a2ef0f 100644
--- a/devito/__init__.py
+++ b/devito/__init__.py
@@ -11,6 +11,7 @@
 from devito.logger import error, warning, info  # noqa
 from devito.parameters import (configuration, init_configuration,  # noqa
                                env_vars_mapper)
+from devito.tools import *  # noqa
 
 
 def clear_cache():
diff --git a/devito/dle/blocking_utils.py b/devito/dle/blocking_utils.py
index bca3d77bb0..77ff736855 100644
--- a/devito/dle/blocking_utils.py
+++ b/devito/dle/blocking_utils.py
@@ -40,9 +40,13 @@ def fold_blockable_tree(node, exclude_innermost=False):
             pairwise_folds = list(zip(*reversed(trees)))
             if any(not is_foldable(j) for j in pairwise_folds):
                 continue
-            # Perform folding
+            # Maybe heuristically exclude innermost Iteration
             if exclude_innermost is True:
                 pairwise_folds = pairwise_folds[:-1]
+            # Perhaps there's nothing to fold
+            if len(pairwise_folds) == 1:
+                continue
+            # Perform folding
             for j in pairwise_folds:
                 root, remainder = j[0], j[1:]
                 folds = [(tuple(y-x for x, y in zip(i.offsets, root.offsets)), i.nodes)
diff --git a/devito/dse/backends/advanced.py b/devito/dse/backends/advanced.py
index 89da9c7d82..637d6f461f 100644
--- a/devito/dse/backends/advanced.py
+++ b/devito/dse/backends/advanced.py
@@ -9,7 +9,7 @@
 from devito.dse.clusterizer import clusterize
 from devito.dse.inspection import estimate_cost
 from devito.dse.manipulation import (common_subexprs_elimination, collect_nested,
-                                     xreplace_constrained)
+                                     xreplace_constrained, compact_temporaries)
 from devito.dse.queries import iq_timeinvariant
 from devito.interfaces import Indexed, ScalarFunction, TensorFunction
 
@@ -34,14 +34,18 @@ def _extract_time_invariants(self, cluster, template, with_cse=True,
         rule = iq_timeinvariant(cluster.trace)
         costmodel = costmodel or (lambda e: estimate_cost(e) > 0)
         processed, found = xreplace_constrained(cluster.exprs, make, rule, costmodel)
-        leaves = [i for i in processed if i not in found]
 
-        # Search for common sub-expressions amongst them (and only them)
         if with_cse:
+            leaves = [i for i in processed if i not in found]
+
+            # Search for common sub-expressions amongst them (and only them)
             make = lambda i: ScalarFunction(name=template(i + len(found))).indexify()
             found = common_subexprs_elimination(found, make)
 
-        return cluster.reschedule(found + leaves)
+            # Some temporaries may be droppable at this point
+            processed = compact_temporaries(found + leaves)
+
+        return cluster.reschedule(processed)
 
     @dse_pass
     def _factorize(self, cluster, *args, **kwargs):
diff --git a/devito/dse/clusterizer.py b/devito/dse/clusterizer.py
index 6773e17b49..e593050efa 100644
--- a/devito/dse/clusterizer.py
+++ b/devito/dse/clusterizer.py
@@ -105,7 +105,8 @@ def optimize(clusters):
         mapper = {}
         temporaries = []
         for k, v in c1.trace.items():
-            if v.is_tensor and not any(v.function in c2.unknown for c2 in clusters):
+            if v.function.is_TensorFunction and\
+                    not any(v.function in c2.unknown for c2 in clusters):
                 for i in c1.tensors[v.function]:
                     # LHS scalarization
                     scalarized = ScalarFunction(name='s%d' % len(mapper)).indexify()
diff --git a/devito/dse/graph.py b/devito/dse/graph.py
index 95e0587dac..db1d6e4858 100644
--- a/devito/dse/graph.py
+++ b/devito/dse/graph.py
@@ -92,10 +92,6 @@ def is_tensor(self):
     def is_scalar(self):
         return not self.is_tensor
 
-    @property
-    def is_dead(self):
-        return self.is_scalar and self.is_terminal and len(self.reads) == 1
-
     def construct(self, rule):
         """
         Create a new temporary starting from ``self`` replacing symbols in
diff --git a/devito/dse/manipulation.py b/devito/dse/manipulation.py
index 483a9621b3..b640a3c29c 100644
--- a/devito/dse/manipulation.py
+++ b/devito/dse/manipulation.py
@@ -10,7 +10,7 @@
 from devito.dse.extended_sympy import Add, Eq, Mul
 from devito.dse.inspection import count, estimate_cost, retrieve_indexed
 from devito.dse.graph import temporaries_graph
-from devito.dse.queries import q_indexed, q_op
+from devito.dse.queries import q_indexed, q_op, q_leaf
 from devito.interfaces import Indexed, TensorFunction
 from devito.tools import as_tuple
 
@@ -105,7 +105,7 @@ def run(expr):
     return run(expr)[0]
 
 
-def xreplace_constrained(exprs, make, rule, costmodel=lambda e: True, repeat=False):
+def xreplace_constrained(exprs, make, rule=None, costmodel=lambda e: True, repeat=False):
     """
     Unlike ``xreplace``, which replaces all objects specified in a mapper,
     this function replaces all objects satisfying two criteria: ::
@@ -122,28 +122,38 @@ def xreplace_constrained(exprs, make, rule, costmodel=lambda e: True, repeat=Fal
     ``a + b`` satisfies the cost model.
 
     :param exprs: The target SymPy expression, or a collection of SymPy expressions.
-    :param make: A function to construct symbols used for replacement.
-                 The function takes as input an integer ID; ID is computed internally
-                 and used as a unique identifier for the constructed symbols.
-    :param rule: The matching rule (a lambda function).
+    :param make: Either a mapper M: K -> V, indicating how to replace an expression
+                 in K with a symbol in V, or a function, used to construct new, unique
+                 symbols. Such a function should take as input a parameter, used to
+                 enumerate the new symbols.
+    :param rule: The matching rule (a lambda function). May be left unspecified if
+                 ``make`` is a mapper.
     :param costmodel: The cost model (a lambda function, optional).
     :param repeat: Repeatedly apply ``xreplace`` until no more replacements are
                    possible (optional, defaults to False).
     """
-
     found = OrderedDict()
     rebuilt = []
 
-    def replace(expr):
-        temporary = found.get(expr)
-        if temporary:
-            return temporary
-        else:
-            temporary = make(replace.c)
-            found[expr] = temporary
-            replace.c += 1
-            return temporary
-    replace.c = 0  # Unique identifier for new temporaries
+    # Define /replace()/ based on the user-provided /make/
+    if isinstance(make, dict):
+        rule = rule if rule is not None else (lambda i: i in make)
+        replace = lambda i: make[i]
+    else:
+        assert callable(make) and callable(rule)
+
+        def replace(expr):
+            if isinstance(make, dict):
+                return make[expr]
+            temporary = found.get(expr)
+            if temporary:
+                return temporary
+            else:
+                temporary = make(replace.c)
+                found[expr] = temporary
+                replace.c += 1
+                return temporary
+        replace.c = 0  # Unique identifier for new temporaries
 
     def run(expr):
         if expr.is_Atom or q_indexed(expr):
@@ -254,9 +264,6 @@ def common_subexprs_elimination(exprs, make, mode='default'):
     mapper = {i.lhs: j.lhs for i, j in zip(mapped, reversed(mapped))}
     processed = [e.xreplace(mapper) for e in processed]
 
-    # Some temporaries may be droppable at this point
-    processed = compact_temporaries(processed)
-
     return processed
 
 
@@ -266,14 +273,16 @@ def compact_temporaries(exprs):
     """
     g = temporaries_graph(exprs)
 
-    mapper = {list(v.reads)[0]: k for k, v in g.items() if v.is_dead}
+    mapper = {k: v.rhs for k, v in g.items()
+              if v.is_scalar and (q_leaf(v.rhs) or v.rhs.is_Function)}
 
     processed = []
     for k, v in g.items():
-        if k in mapper:
-            processed.append(Eq(mapper[k], v.rhs))
-        elif not v.is_dead:
-            processed.append(v.xreplace(mapper))
+        if k not in mapper:
+            # The temporary /v/ is retained, and substitutions may be applied
+            handle, _ = xreplace_constrained(v, mapper, repeat=True)
+            assert len(handle) == 1
+            processed.extend(handle)
 
     return processed
 
diff --git a/devito/tools.py b/devito/tools.py
index c126b43765..d6898c58be 100644
--- a/devito/tools.py
+++ b/devito/tools.py
@@ -1,6 +1,7 @@
 import os
 import ctypes
-from collections import Callable, Iterable, OrderedDict
+from collections import Callable, Iterable, OrderedDict, Hashable
+from functools import partial
 try:
     from itertools import izip_longest as zip_longest
 except ImportError:
@@ -9,6 +10,8 @@
 
 import numpy as np
 
+__all__ = ['memoized']
+
 
 def as_tuple(item, type=None, length=None):
     """
@@ -210,3 +213,39 @@ def __enter__(self):
 
     def __exit__(self, etype, value, traceback):
         os.chdir(self.savedPath)
+
+
+class memoized(object):
+    """
+    Decorator. Caches a function's return value each time it is called.
+    If called later with the same arguments, the cached value is returned
+    (not reevaluated).
+
+    Adapted from: ::
+
+        https://wiki.python.org/moin/PythonDecoratorLibrary#Memoize
+    """
+
+    def __init__(self, func):
+        self.func = func
+        self.cache = {}
+
+    def __call__(self, *args):
+        if not isinstance(args, Hashable):
+            # Uncacheable, a list, for instance.
+            # Better to not cache than blow up.
+            return self.func(*args)
+        if args in self.cache:
+            return self.cache[args]
+        else:
+            value = self.func(*args)
+            self.cache[args] = value
+            return value
+
+    def __repr__(self):
+        """Return the function's docstring."""
+        return self.func.__doc__
+
+    def __get__(self, obj, objtype):
+        """Support instance methods."""
+        return partial(self.__call__, obj)
diff --git a/examples/seismic/acoustic/acoustic_example.py b/examples/seismic/acoustic/acoustic_example.py
index 29ae2d6cb7..332f0d0a7a 100644
--- a/examples/seismic/acoustic/acoustic_example.py
+++ b/examples/seismic/acoustic/acoustic_example.py
@@ -21,10 +21,10 @@ def smooth10(vel, shape):
     return out
 
 
-def acoustic_setup(dimensions=(50, 50, 50), spacing=(15.0, 15.0, 15.0),
+def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0),
                    tn=500., time_order=2, space_order=4, nbpml=10, **kwargs):
-    nrec = dimensions[0]
-    model = demo_model('layers', ratio=3, shape=dimensions,
+    nrec = shape[0]
+    model = demo_model('layers-isotropic', ratio=3, shape=shape,
                        spacing=spacing, nbpml=nbpml)
 
     # Derive timestepping from model spacing
@@ -50,11 +50,11 @@ def acoustic_setup(dimensions=(50, 50, 50), spacing=(15.0, 15.0, 15.0),
     return solver
 
 
-def run(dimensions=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
+def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
         time_order=2, space_order=4, nbpml=40, full_run=False,
         autotune=False, **kwargs):
 
-    solver = acoustic_setup(dimensions=dimensions, spacing=spacing,
+    solver = acoustic_setup(dimensions=shape, spacing=spacing,
                             nbpml=nbpml, tn=tn, space_order=space_order,
                             time_order=time_order, **kwargs)
 
@@ -89,22 +89,25 @@ def run(dimensions=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
                         type=int, help="Space order of the simulation")
     parser.add_argument("--nbpml", default=40,
                         type=int, help="Number of PML layers around the domain")
-    parser.add_argument("-dse", default='advanced',
-                        type=str, help="DSE backend choice")
-    parser.add_argument("-dle", default='advanced',
-                        type=str, help="DLE backend choice")
+    parser.add_argument("-dse", "-dse", default="advanced",
+                        choices=["noop", "basic", "advanced",
+                                 "speculative", "aggressive"],
+                        help="Devito symbolic engine (DSE) mode")
+    parser.add_argument("-dle", default="advanced",
+                        choices=["noop", "advanced", "speculative"],
+                        help="Devito loop engine (DSE) mode")
     args = parser.parse_args()
 
     # 3D preset parameters
     if args.dim2:
-        dimensions = (150, 150)
+        shape = (150, 150)
         spacing = (15.0, 15.0)
         tn = 750.0
     else:
-        dimensions = (50, 50, 50)
+        shape = (50, 50, 50)
         spacing = (20.0, 20.0, 20.0)
         tn = 250.0
 
-    run(dimensions=dimensions, spacing=spacing, nbpml=args.nbpml, tn=tn,
+    run(shape=shape, spacing=spacing, nbpml=args.nbpml, tn=tn,
         space_order=args.space_order, time_order=args.time_order,
         autotune=args.autotune, dse=args.dse, dle=args.dle, full_run=args.full)
diff --git a/examples/seismic/benchmark.py b/examples/seismic/benchmark.py
index 2a1eefd890..c1c2af1fb7 100644
--- a/examples/seismic/benchmark.py
+++ b/examples/seismic/benchmark.py
@@ -30,7 +30,7 @@
                         choices=["acoustic", "tti"], help="Problem")
 
     simulation = parser.add_argument_group("Simulation")
-    simulation.add_argument("-d", "--dimensions", nargs=3, default=[50, 50, 50],
+    simulation.add_argument("-d", "--shape", nargs=3, default=[50, 50, 50],
                             type=int, help="Number of grid points along each axis",
                             metavar=("dim1", "dim2", "dim3"))
     simulation.add_argument("-s", "--spacing", nargs=3, default=[20.0, 20.0, 20.0],
@@ -91,7 +91,7 @@
     del parameters["max_flops"]
     del parameters["point_runtime"]
 
-    parameters["dimensions"] = tuple(parameters["dimensions"])
+    parameters["shape"] = tuple(parameters["shape"])
     parameters["spacing"] = tuple(parameters["spacing"])
 
     if args.execmode == "run":
diff --git a/examples/seismic/model.py b/examples/seismic/model.py
index e8d7962130..fd04234b17 100644
--- a/examples/seismic/model.py
+++ b/examples/seismic/model.py
@@ -19,7 +19,7 @@ def demo_model(preset, **kwargs):
                     filepath. Requires the ``opesci/data`` repository
                     to be available on your machine.
     """
-    if preset.lower() in ['constant']:
+    if preset.lower() in ['constant-isotropic']:
         # A constant single-layer model in a 2D or 3D domain
         # with velocity 1.5km/s.
         shape = kwargs.pop('shape', (101, 101))
@@ -31,7 +31,29 @@ def demo_model(preset, **kwargs):
         return Model(vp=vp, origin=origin, shape=shape,
                      spacing=spacing, nbpml=nbpml, **kwargs)
 
-    elif preset.lower() in ['layers', 'twolayer', '2layer']:
+    elif preset.lower() in ['constant-tti']:
+        # A constant single-layer model in a 2D or 3D domain
+        # with velocity 1.5km/s.
+        shape = kwargs.pop('shape', (101, 101))
+        spacing = kwargs.pop('spacing', tuple([10. for _ in shape]))
+        origin = kwargs.pop('origin', tuple([0. for _ in shape]))
+        nbpml = kwargs.pop('nbpml', 10)
+        v = np.empty(shape, dtype=np.float32)
+        v[:] = 1.5
+        epsilon = .3*np.ones(shape)
+        delta = .2*np.ones(shape)
+        theta = .7*np.ones(shape)
+        phi = None
+        if len(shape) > 2:
+            phi = .35*np.ones(shape)
+
+        return Model(vp=v, origin=origin, shape=shape,
+                     spacing=spacing, nbpml=nbpml,
+                     epsilon=epsilon, delta=delta, theta=theta, phi=phi,
+                     **kwargs)
+
+    elif preset.lower() in ['layers-isotropic', 'twolayer-isotropic',
+                            '2layer-isotropic']:
         # A two-layer model in a 2D or 3D domain with two different
         # velocities split across the height dimension:
         # By default, the top part of the domain has 1.5 km/s,
@@ -52,7 +74,37 @@ def demo_model(preset, **kwargs):
         return Model(vp=v, origin=origin, shape=shape,
                      spacing=spacing, nbpml=nbpml, **kwargs)
 
-    elif preset.lower() in ['circle']:
+    elif preset.lower() in ['layers-tti', 'twolayer-tti', '2layer-tti']:
+        # A two-layer model in a 2D or 3D domain with two different
+        # velocities split across the height dimension:
+        # By default, the top part of the domain has 1.5 km/s,
+        # and the bottom part of the domain has 2.5 km/s.\
+        shape = kwargs.pop('shape', (101, 101))
+        spacing = kwargs.pop('spacing', tuple([10. for _ in shape]))
+        origin = kwargs.pop('origin', tuple([0. for _ in shape]))
+        nbpml = kwargs.pop('nbpml', 10)
+        ratio = kwargs.pop('ratio', 2)
+        vp_top = kwargs.pop('vp_top', 1.5)
+        vp_bottom = kwargs.pop('vp_bottom', 2.5)
+
+        # Define a velocity profile in km/s
+        v = np.empty(shape, dtype=np.float32)
+        v[:] = vp_top  # Top velocity (background)
+        v[..., int(shape[-1] / ratio):] = vp_bottom  # Bottom velocity
+
+        epsilon = .3*(v - 1.5)
+        delta = .2*(v - 1.5)
+        theta = .5*(v - 1.5)
+        phi = None
+        if len(shape) > 2:
+            phi = .1*(v - 1.5)
+
+        return Model(vp=v, origin=origin, shape=shape,
+                     spacing=spacing, nbpml=nbpml,
+                     epsilon=epsilon, delta=delta, theta=theta, phi=phi,
+                     **kwargs)
+
+    elif preset.lower() in ['circle-isotropic']:
         # A simple circle in a 2D domain with a background velocity.
         # By default, the circle velocity is 2.5 km/s,
         # and the background veloity is 3.0 km/s.
@@ -77,7 +129,7 @@ def demo_model(preset, **kwargs):
         return Model(vp=v, origin=origin, shape=shape,
                      spacing=spacing, nbpml=nbpml)
 
-    elif preset.lower() in ['marmousi', 'marmousi2d']:
+    elif preset.lower() in ['marmousi-isotropic', 'marmousi2d-isotropic']:
         shape = (1601, 401)
         spacing = (7.5, 7.5)
         origin = (0., 0.)
@@ -98,6 +150,89 @@ def demo_model(preset, **kwargs):
         return Model(vp=v, origin=origin, shape=v.shape,
                      spacing=spacing, nbpml=20)
 
+    elif preset.lower() in ['marmousi-tti2d', 'marmousi2d-tti']:
+
+        shape_full = (201, 201, 70)
+        shape = (201, 70)
+        spacing = (10., 10.)
+        origin = (0., 0.)
+        nbpml = kwargs.pop('nbpml', 20)
+
+        # Read 2D Marmousi model from opesc/data repo
+        data_path = kwargs.pop('data_path', None)
+        if data_path is None:
+            error("Path to opesci/data not found! Please specify with "
+                  "'data_path=<path/to/opesci/data>'")
+            raise ValueError("Path to model data unspecified")
+        path = os.path.join(data_path, 'marmousi3D/vp_marmousi_bi')
+
+        # velocity
+        vp = 1e-3 * np.fromfile(os.path.join(data_path, 'marmousi3D/MarmousiVP.raw'),
+                                dtype='float32', sep="")
+        vp = vp.reshape(shape_full)
+        vp = vp[101, :, :]
+        # Epsilon, in % in file, resale between 0 and 1
+        epsilon = np.fromfile(os.path.join(data_path, 'marmousi3D/MarmousiEps.raw'),
+                              dtype='float32', sep="") * 1e-2
+        epsilon = epsilon.reshape(shape_full)
+        epsilon = epsilon[101, :, :]
+        # Delta, in % in file, resale between 0 and 1
+        delta = np.fromfile(os.path.join(data_path, 'marmousi3D/MarmousiDelta.raw'),
+                            dtype='float32', sep="") * 1e-2
+        delta = delta.reshape(shape_full)
+        delta = delta[101, :, :]
+        # Theta, in degrees in file, resale in radian
+        theta = np.fromfile(os.path.join(data_path, 'marmousi3D/MarmousiTilt.raw'),
+                            dtype='float32', sep="")
+        theta = np.float32(np.pi / 180 * theta.reshape(shape_full))
+        theta = theta[101, :, :]
+
+        return Model(vp=vp, origin=origin, shape=shape,
+                     spacing=spacing, nbpml=nbpml,
+                     epsilon=epsilon, delta=delta, theta=theta,
+                     **kwargs)
+
+    elif preset.lower() in ['marmousi-tti3d', 'marmousi3d-tti']:
+
+        shape = (201, 201, 70)
+        spacing = (10., 10., 10.)
+        origin = (0., 0., 0.)
+        nbpml = kwargs.pop('nbpml', 20)
+
+        # Read 2D Marmousi model from opesc/data repo
+        data_path = kwargs.pop('data_path', None)
+        if data_path is None:
+            error("Path to opesci/data not found! Please specify with "
+                  "'data_path=<path/to/opesci/data>'")
+            raise ValueError("Path to model data unspecified")
+        path = os.path.join(data_path, 'marmousi3D/vp_marmousi_bi')
+
+        # Velcoity
+        vp = 1e-3 * np.fromfile(os.path.join(data_path, 'marmousi3D/MarmousiVP.raw'),
+                                dtype='float32', sep="")
+        vp = vp.reshape(shape)
+        # Epsilon, in % in file, resale between 0 and 1
+        epsilon = np.fromfile(os.path.join(data_path, 'marmousi3D/MarmousiEps.raw'),
+                              dtype='float32', sep="") * 1e-2
+        epsilon = epsilon.reshape(shape)
+        # Delta, in % in file, resale between 0 and 1
+        delta = np.fromfile(os.path.join(data_path, 'marmousi3D/MarmousiDelta.raw'),
+                            dtype='float32', sep="") * 1e-2
+        delta = delta.reshape(shape)
+        # Theta, in degrees in file, resale in radian
+        theta = np.fromfile(os.path.join(data_path, 'marmousi3D/MarmousiTilt.raw'),
+                            dtype='float32', sep="")
+        theta = np.float32(np.pi / 180 * theta.reshape(shape))
+        # Phi, in degrees in file, resale in radian
+        phi = np.fromfile(os.path.join(data_path, 'marmousi3D/Azimuth.raw'),
+                          dtype='float32', sep="")
+        phi = np.float32(np.pi / 180 * phi.reshape(shape))
+
+        return Model(vp=vp, origin=origin, shape=shape,
+                     spacing=spacing, nbpml=nbpml,
+                     epsilon=epsilon, delta=delta, theta=theta, phi=phi,
+                     **kwargs)
+
     else:
         error('Unknown model preset name %s' % preset)
 
diff --git a/examples/seismic/tti/operators.py b/examples/seismic/tti/operators.py
index de784a3802..f2ed91cbb3 100644
--- a/examples/seismic/tti/operators.py
+++ b/examples/seismic/tti/operators.py
@@ -6,18 +6,362 @@
 from devito.dimension import x, y, z, t, time
 
 
-def ForwardOperator(model, source, receiver, time_order=2, space_order=4,
-                    save=False, **kwargs):
+def Gxx_shifted(field, costheta, sintheta, cosphi, sinphi, space_order):
+    """
+    3D rotated second order derivative in the direction x as an average of
+    two non-centered rotated second order derivative in the direction x
+    :param field: symbolic data whose derivative we are computing
+    :param costheta: cosine of the tilt angle
+    :param sintheta:  sine of the tilt angle
+    :param cosphi: cosine of the azymuth angle
+    :param sinphi: sine of the azymuth angle
+    :param space_order: discretization order
+    :return: rotated second order derivative wrt x
+    """
+    Gx1 = (costheta * cosphi * field.dx + costheta * sinphi * field.dyr -
+           sintheta * field.dzr)
+    Gxx1 = (first_derivative(Gx1 * costheta * cosphi,
+                             dim=x, side=centered, order=space_order,
+                             matvec=transpose) +
+            first_derivative(Gx1 * costheta * sinphi,
+                             dim=y, side=right, order=space_order,
+                             matvec=transpose) -
+            first_derivative(Gx1 * sintheta,
+                             dim=z, side=right, order=space_order,
+                             matvec=transpose))
+    Gx2 = (costheta * cosphi * field.dxr + costheta * sinphi * field.dy -
+           sintheta * field.dz)
+    Gxx2 = (first_derivative(Gx2 * costheta * cosphi,
+                             dim=x, side=right, order=space_order,
+                             matvec=transpose) +
+            first_derivative(Gx2 * costheta * sinphi,
+                             dim=y, side=centered, order=space_order,
+                             matvec=transpose) -
+            first_derivative(Gx2 * sintheta,
+                             dim=z, side=centered, order=space_order,
+                             matvec=transpose))
+    return -.5 * (Gxx1 + Gxx2)
+
+
+def Gxx_shifted_2d(field, costheta, sintheta, space_order):
+    """
+    2D rotated second order derivative in the direction x as an average of
+    two non-centered rotated second order derivative in the direction x
+    :param field: symbolic data whose derivative we are computing
+    :param costheta: cosine of the tilt angle
+    :param sintheta:  sine of the tilt angle
+    :param space_order: discretization order
+    :return: rotated second order derivative wrt x
+    """
+    Gx1 = (costheta * field.dxr - sintheta * field.dy)
+    Gxx1 = (first_derivative(Gx1 * costheta, dim=x,
+                             side=right, order=space_order,
+                             matvec=transpose) -
+            first_derivative(Gx1 * sintheta, dim=y,
+                             side=centered, order=space_order,
+                             matvec=transpose))
+    Gx2p = (costheta * field.dx - sintheta * field.dyr)
+    Gxx2 = (first_derivative(Gx2p * costheta, dim=x,
+                             side=centered, order=space_order,
+                             matvec=transpose) -
+            first_derivative(Gx2p * sintheta, dim=y,
+                             side=right, order=space_order,
+                             matvec=transpose))
+
+    return -.5 * (Gxx1 + Gxx2)
+
+
+def Gyy_shifted(field, cosphi, sinphi, space_order):
+    """
+    3D rotated second order derivative in the direction y as an average of
+    two non-centered rotated second order derivative in the direction y
+    :param field: symbolic data whose derivative we are computing
+    :param cosphi: cosine of the azymuth angle
+    :param sinphi: sine of the azymuth angle
+    :param space_order: discretization order
+    :return: rotated second order derivative wrt y
+    """
+    Gyp = (sinphi * field.dx - cosphi * field.dyr)
+    Gyy = (first_derivative(Gyp * sinphi,
+                            dim=x, side=centered, order=space_order,
+                            matvec=transpose) -
+           first_derivative(Gyp * cosphi,
+                            dim=y, side=right, order=space_order,
+                            matvec=transpose))
+    Gyp2 = (sinphi * field.dxr - cosphi * field.dy)
+    Gyy2 = (first_derivative(Gyp2 * sinphi,
+                             dim=x, side=right, order=space_order,
+                             matvec=transpose) -
+            first_derivative(Gyp2 * cosphi,
+                             dim=y, side=centered, order=space_order,
+                             matvec=transpose))
+    return -.5 * (Gyy + Gyy2)
+
+
+def Gzz_shited(field, costheta, sintheta, cosphi, sinphi, space_order):
+    """
+    3D rotated second order derivative in the direction z as an average of
+    two non-centered rotated second order derivative in the direction z
+    :param field: symbolic data whose derivative we are computing
+    :param costheta: cosine of the tilt angle
+    :param sintheta:  sine of the tilt angle
+    :param cosphi: cosine of the azymuth angle
+    :param sinphi: sine of the azymuth angle
+    :param space_order: discretization order
+    :return: rotated second order derivative wrt z
+    """
+    Gzr = (sintheta * cosphi * field.dx + sintheta * sinphi * field.dyr +
+           costheta * field.dzr)
+    Gzz = (first_derivative(Gzr * sintheta * cosphi,
+                            dim=x, side=centered, order=space_order,
+                            matvec=transpose) +
+           first_derivative(Gzr * sintheta * sinphi,
+                            dim=y, side=right, order=space_order,
+                            matvec=transpose) +
+           first_derivative(Gzr * costheta,
+                            dim=z, side=right, order=space_order,
+                            matvec=transpose))
+    Gzr2 = (sintheta * cosphi * field.dxr + sintheta * sinphi * field.dy +
+            costheta * field.dz)
+    Gzz2 = (first_derivative(Gzr2 * sintheta * cosphi,
+                             dim=x, side=right, order=space_order,
+                             matvec=transpose) +
+            first_derivative(Gzr2 * sintheta * sinphi,
+                             dim=y, side=centered, order=space_order,
+                             matvec=transpose) +
+            first_derivative(Gzr2 * costheta,
+                             dim=z, side=centered, order=space_order,
+                             matvec=transpose))
+    return -.5 * (Gzz + Gzz2)
+
+
+def Gzz_shited_2d(field, costheta, sintheta, space_order):
+    """
+    2D rotated second order derivative in the direction z as an average of
+    two non-centered rotated second order derivative in the direction z
+    :param field: symbolic data whose derivative we are computing
+    :param costheta: cosine of the tilt
+    :param sintheta:  sine of the tilt
+    :param space_order: discretization order
+    :return: rotated second order derivative wrt z
+    """
+    Gz1r = (sintheta * field.dxr + costheta * field.dy)
+    Gzz1 = (first_derivative(Gz1r * sintheta, dim=x,
+                             side=right, order=space_order,
+                             matvec=transpose) +
+            first_derivative(Gz1r * costheta, dim=y,
+                             side=centered, order=space_order,
+                             matvec=transpose))
+    Gz2r = (sintheta * field.dx + costheta * field.dyr)
+    Gzz2 = (first_derivative(Gz2r * sintheta, dim=x,
+                             side=centered, order=space_order,
+                             matvec=transpose) +
+            first_derivative(Gz2r * costheta, dim=y,
+                             side=right, order=space_order,
+                             matvec=transpose))
+
+    return -.5 * (Gzz1 + Gzz2)
+
+
+def Gzz_centered(field, costheta, sintheta, cosphi, sinphi, space_order):
+    """
+    3D rotated second order derivative in the direction z
+    :param field: symbolic data whose derivative we are computing
+    :param costheta: cosine of the tilt angle
+    :param sintheta:  sine of the tilt angle
+    :param cosphi: cosine of the azymuth angle
+    :param sinphi: sine of the azymuth angle
+    :param space_order: discretization order
+    :return: rotated second order derivative wrt z
+    """
+    order1 = space_order / 2
+    Gz = -(sintheta * cosphi * first_derivative(field, dim=x,
+                                                side=centered, order=order1) +
+           sintheta * sinphi * first_derivative(field, dim=y,
+                                                side=centered, order=order1) +
+           costheta * first_derivative(field, dim=z,
+                                       side=centered, order=order1))
+    Gzz = (first_derivative(Gz * sintheta * cosphi,
+                            dim=x, side=centered, order=order1,
+                            matvec=transpose) +
+           first_derivative(Gz * sintheta * sinphi,
+                            dim=y, side=centered, order=order1,
+                            matvec=transpose) +
+           first_derivative(Gz * costheta,
+                            dim=z, side=centered, order=order1,
+                            matvec=transpose))
+    return Gzz
+
+
+def Gzz_centered_2d(field, costheta, sintheta, space_order):
+    """
+    2D rotated second order derivative in the direction z
+    :param field: symbolic data whose derivative we are computing
+    :param costheta: cosine of the tilt angle
+    :param sintheta:  sine of the tilt angle
+    :param space_order: discretization order
+    :return: rotated second order derivative wrt z
+    """
+    order1 = space_order / 2
+    Gz = -(sintheta * first_derivative(field, dim=x, side=centered, order=order1) +
+           costheta * first_derivative(field, dim=y, side=centered, order=order1))
+    Gzz = (first_derivative(Gz * sintheta, dim=x,
+                            side=centered, order=order1,
+                            matvec=transpose) +
+           first_derivative(Gz * costheta, dim=y,
+                            side=centered, order=order1,
+                            matvec=transpose))
+    return Gzz
+
+
+# Centered case produces directly Gxx + Gyy
+def Gxxyy_centered(field, costheta, sintheta, cosphi, sinphi, space_order):
+    """
+    Sum of the 3D rotated second order derivative in the direction x and y.
+    As the Laplacian is rotation invariant, it is computed as the conventional
+    Laplacian minus the second order rotated second order derivative in the direction z
+    Gxx + Gyy = field.laplace - Gzz
+    :param field: symbolic data whose derivative we are computing
+    :param costheta: cosine of the tilt angle
+    :param sintheta:  sine of the tilt angle
+    :param cosphi: cosine of the azymuth angle
+    :param sinphi: sine of the azymuth angle
+    :param space_order: discretization order
+    :return: Sum of the 3D rotated second order derivative in the direction x and y
+    """
+    Gzz = Gzz_centered(field, costheta, sintheta, cosphi, sinphi, space_order)
+    return field.laplace - Gzz
+
+
+def Gxx_centered_2d(field, costheta, sintheta, space_order):
+    """
+    2D rotated second order derivative in the direction x.
+    As the Laplacian is rotation invariant, it is computed as the conventional
+    Laplacian minus the second order rotated second order derivative in the direction z
+    Gxx = field.laplace - Gzz
+    :param field: symbolic data whose derivative we are computing
+    :param costheta: cosine of the tilt angle
+    :param sintheta:  sine of the tilt angle
+    :param cosphi: cosine of the azymuth angle
+    :param sinphi: sine of the azymuth angle
+    :param space_order: discretization order
+    :return: Sum of the 3D rotated second order derivative in the direction x
+    """
+    return field.laplace - Gzz_centered_2d(field, costheta, sintheta, space_order)
+
+
+def kernel_shited_2d(u, v, costheta, sintheta, cosphi, sinphi, space_order):
+    """
+    TTI finite difference kernel. The equation we solve is:
+
+    u.dt2 = (1+2 *epsilon) (Gxx(u)) + sqrt(1+ 2*delta) Gzz(v)
+    v.dt2 = sqrt(1+ 2*delta) (Gxx(u)) +  Gzz(v)
+
+    where epsilon and delta are the thomsen parameters. This function computes
+    H0 = Gxx(u) + Gyy(u)
+    Hz = Gzz(v)
+
+    :param u: first TTI field
+    :param v: second TTI field
+    :param costheta: cosine of the tilt angle
+    :param sintheta:  sine of the tilt angle
+    :param cosphi: cosine of the azymuth angle, has to be 0 in 2D
+    :param sinphi: sine of the azymuth angle, has to be 0 in 2D
+    :param space_order: discretization order
+    :return: u and v component of the rotated Laplacian in 2D
     """
-    Constructor method for the forward modelling operator in an acoustic media
+    Gxx = Gxx_shifted_2d(u, costheta, sintheta, space_order)
+    Gzz = Gzz_shited_2d(v, costheta, sintheta, space_order)
+    return Gxx, Gzz
+
 
-    :param model: :class:`Model` object containing the physical parameters
-    :param src: None ot IShot() (not currently supported properly)
-    :param data: IShot() object containing the acquisition geometry and field data
-    :param: time_order: Time discretization order
-    :param: spc_order: Space discretization order
-    :param: u_ini : wavefield at the three first time step for non-zero initial condition
+def kernel_shited_3d(u, v, costheta, sintheta, cosphi, sinphi, space_order):
     """
+    TTI finite difference kernel. The equation we solve is:
+
+    u.dt2 = (1+2 *epsilon) (Gxx(u)+Gyy(u)) + sqrt(1+ 2*delta) Gzz(v)
+    v.dt2 = sqrt(1+ 2*delta) (Gxx(u)+Gyy(u)) +  Gzz(v)
+
+    where epsilon and delta are the thomsen parameters. This function computes
+    H0 = Gxx(u) + Gyy(u)
+    Hz = Gzz(v)
+
+    :param u: first TTI field
+    :param v: second TTI field
+    :param costheta: cosine of the tilt angle
+    :param sintheta:  sine of the tilt angle
+    :param cosphi: cosine of the azymuth angle
+    :param sinphi: sine of the azymuth angle
+    :param space_order: discretization order
+    :return: u and v component of the rotated Laplacian in 3D
+    """
+    Gxx = Gxx_shifted(u, costheta, sintheta, cosphi, sinphi, space_order)
+    Gyy = Gyy_shifted(u, cosphi, sinphi, space_order)
+    Gzz = Gzz_shited(v, costheta, sintheta, cosphi, sinphi, space_order)
+    return Gxx + Gyy, Gzz
+
+
+def kernel_centered_2d(u, v, costheta, sintheta, cosphi, sinphi, space_order):
+    """
+    TTI finite difference kernel. The equation we solve is:
+
+    u.dt2 = (1+2 *epsilon) (Gxx(u)) + sqrt(1+ 2*delta) Gzz(v)
+    v.dt2 = sqrt(1+ 2*delta) (Gxx(u)) +  Gzz(v)
+
+    where epsilon and delta are the thomsen parameters. This function computes
+    H0 = Gxx(u) + Gyy(u)
+    Hz = Gzz(v)
+
+    :param u: first TTI field
+    :param v: second TTI field
+    :param costheta: cosine of the tilt angle
+    :param sintheta:  sine of the tilt angle
+    :param cosphi: cosine of the azymuth angle, has to be 0 in 2D
+    :param sinphi: sine of the azymuth angle, has to be 0 in 2D
+    :param space_order: discretization order
+    :return: u and v component of the rotated Laplacian in 2D
+    """
+    Gxx = Gxx_centered_2d(u, costheta, sintheta, space_order)
+    Gzz = Gzz_centered_2d(v, costheta, sintheta, space_order)
+    return Gxx, Gzz
+
+
+def kernel_centered_3d(u, v, costheta, sintheta, cosphi, sinphi, space_order):
+    """
+    TTI finite difference kernel. The equation we solve is:
+
+    u.dt2 = (1+2 *epsilon) (Gxx(u)+Gyy(u)) + sqrt(1+ 2*delta) Gzz(v)
+    v.dt2 = sqrt(1+ 2*delta) (Gxx(u)+Gyy(u)) +  Gzz(v)
+
+    where epsilon and delta are the thomsen parameters. This function computes
+    H0 = Gxx(u) + Gyy(u)
+    Hz = Gzz(v)
+
+    :param u: first TTI field
+    :param v: second TTI field
+    :param costheta: cosine of the tilt angle
+    :param sintheta:  sine of the tilt angle
+    :param cosphi: cosine of the azymuth angle
+    :param sinphi: sine of the azymuth angle
+    :param space_order: discretization order
+    :return: u and v component of the rotated Laplacian in 2D
+    """
+    Gxx = Gxxyy_centered(u, costheta, sintheta, cosphi, sinphi, space_order)
+    Gzz = Gzz_centered(v, costheta, sintheta, cosphi, sinphi, space_order)
+    return Gxx, Gzz
+
+
+def ForwardOperator(model, source, receiver, time_order=2, space_order=4,
+                    save=False, kernel='centered', **kwargs):
+    """
+       Constructor method for the forward modelling operator in an acoustic media
+
+       :param model: :class:`Model` object containing the physical parameters
+       :param src: None ot IShot() (not currently supported properly)
+       :param data: IShot() object containing the acquisition geometry and field data
+       :param: time_order: Time discretization order
+       :param: spc_order: Space discretization order
+       """
     dt = model.critical_dt
 
     m, damp, epsilon, delta, theta, phi = (model.m, model.damp, model.epsilon,
@@ -35,124 +379,42 @@ def ForwardOperator(model, source, receiver, time_order=2, space_order=4,
     rec = Receiver(name='rec', ntime=receiver.nt, ndim=receiver.ndim,
                    npoint=receiver.npoint)
 
+    # Tilt and azymuth setup
     ang0 = cos(theta)
     ang1 = sin(theta)
+    ang2 = 0
+    ang3 = 0
     if len(model.shape) == 3:
         ang2 = cos(phi)
         ang3 = sin(phi)
 
-        # Derive stencil from symbolic equation
-        Gyp = (ang3 * u.dx - ang2 * u.dyr)
-        Gyy = (first_derivative(Gyp * ang3,
-                                dim=x, side=centered, order=space_order,
-                                matvec=transpose) -
-               first_derivative(Gyp * ang2,
-                                dim=y, side=right, order=space_order,
-                                matvec=transpose))
-        Gyp2 = (ang3 * u.dxr - ang2 * u.dy)
-        Gyy2 = (first_derivative(Gyp2 * ang3,
-                                 dim=x, side=right, order=space_order,
-                                 matvec=transpose) -
-                first_derivative(Gyp2 * ang2,
-                                 dim=y, side=centered, order=space_order,
-                                 matvec=transpose))
-
-        Gxp = (ang0 * ang2 * u.dx + ang0 * ang3 * u.dyr - ang1 * u.dzr)
-        Gzr = (ang1 * ang2 * v.dx + ang1 * ang3 * v.dyr + ang0 * v.dzr)
-        Gxx = (first_derivative(Gxp * ang0 * ang2,
-                                dim=x, side=centered, order=space_order,
-                                matvec=transpose) +
-               first_derivative(Gxp * ang0 * ang3,
-                                dim=y, side=right, order=space_order,
-                                matvec=transpose) -
-               first_derivative(Gxp * ang1,
-                                dim=z, side=right, order=space_order,
-                                matvec=transpose))
-        Gzz = (first_derivative(Gzr * ang1 * ang2,
-                                dim=x, side=centered, order=space_order,
-                                matvec=transpose) +
-               first_derivative(Gzr * ang1 * ang3,
-                                dim=y, side=right, order=space_order,
-                                matvec=transpose) +
-               first_derivative(Gzr * ang0,
-                                dim=z, side=right, order=space_order,
-                                matvec=transpose))
-        Gxp2 = (ang0 * ang2 * u.dxr + ang0 * ang3 * u.dy - ang1 * u.dz)
-        Gzr2 = (ang1 * ang2 * v.dxr + ang1 * ang3 * v.dy + ang0 * v.dz)
-        Gxx2 = (first_derivative(Gxp2 * ang0 * ang2,
-                                 dim=x, side=right, order=space_order,
-                                 matvec=transpose) +
-                first_derivative(Gxp2 * ang0 * ang3,
-                                 dim=y, side=centered, order=space_order,
-                                 matvec=transpose) -
-                first_derivative(Gxp2 * ang1,
-                                 dim=z, side=centered, order=space_order,
-                                 matvec=transpose))
-        Gzz2 = (first_derivative(Gzr2 * ang1 * ang2,
-                                 dim=x, side=right, order=space_order,
-                                 matvec=transpose) +
-                first_derivative(Gzr2 * ang1 * ang3,
-                                 dim=y, side=centered, order=space_order,
-                                 matvec=transpose) +
-                first_derivative(Gzr2 * ang0,
-                                 dim=z, side=centered, order=space_order,
-                                 matvec=transpose))
-        Hp = -(.5*Gxx + .5*Gxx2 + .5 * Gyy + .5*Gyy2)
-        Hzr = -(.5*Gzz + .5 * Gzz2)
-
-    else:
-        Gx1p = (ang0 * u.dxr - ang1 * u.dy)
-        Gz1r = (ang1 * v.dxr + ang0 * v.dy)
-        Gxx1 = (first_derivative(Gx1p * ang0, dim=x,
-                                 side=right, order=space_order,
-                                 matvec=transpose) -
-                first_derivative(Gx1p * ang1, dim=y,
-                                 side=centered, order=space_order,
-                                 matvec=transpose))
-        Gzz1 = (first_derivative(Gz1r * ang1, dim=x,
-                                 side=right, order=space_order,
-                                 matvec=transpose) +
-                first_derivative(Gz1r * ang0, dim=y,
-                                 side=centered, order=space_order,
-                                 matvec=transpose))
-        Gx2p = (ang0 * u.dx - ang1 * u.dyr)
-        Gz2r = (ang1 * v.dx + ang0 * v.dyr)
-        Gxx2 = (first_derivative(Gx2p * ang0, dim=x,
-                                 side=centered, order=space_order,
-                                 matvec=transpose) -
-                first_derivative(Gx2p * ang1, dim=y,
-                                 side=right, order=space_order,
-                                 matvec=transpose))
-        Gzz2 = (first_derivative(Gz2r * ang1, dim=x,
-                                 side=centered, order=space_order,
-                                 matvec=transpose) +
-                first_derivative(Gz2r * ang0, dim=y,
-                                 side=right, order=space_order,
-                                 matvec=transpose))
-
-        Hp = -(.5 * Gxx1 + .5 * Gxx2)
-        Hzr = -(.5 * Gzz1 + .5 * Gzz2)
-
+    FD_kernel = kernels[(kernel, len(model.shape))]
+    H0, Hz = FD_kernel(u, v, ang0, ang1, ang2, ang3, space_order)
     s = t.spacing
-
+    # Stencils
     stencilp = 1.0 / (2.0 * m + s * damp) * \
         (4.0 * m * u + (s * damp - 2.0 * m) *
-         u.backward + 2.0 * s**2 * (epsilon * Hp + delta * Hzr))
+         u.backward + 2.0 * s ** 2 * (epsilon * H0 + delta * Hz))
     stencilr = 1.0 / (2.0 * m + s * damp) * \
         (4.0 * m * v + (s * damp - 2.0 * m) *
-         v.backward + 2.0 * s**2 * (delta * Hp + Hzr))
-
-    # Add substitutions for spacing (temporal and spatial)
+         v.backward + 2.0 * s ** 2 * (delta * H0 + Hz))
     first_stencil = Eq(u.forward, stencilp)
     second_stencil = Eq(v.forward, stencilr)
     stencils = [first_stencil, second_stencil]
 
+    # Source and receivers
     stencils += src.inject(field=u.forward, expr=src * dt * dt / m,
                            offset=model.nbpml)
     stencils += src.inject(field=v.forward, expr=src * dt * dt / m,
                            offset=model.nbpml)
     stencils += rec.interpolate(expr=u + v, offset=model.nbpml)
+    # Add substitutions for spacing (temporal and spatial)
     subs = dict([(t.spacing, dt)] + [(time.spacing, dt)] +
                 [(i.spacing, model.get_spacing()[j]) for i, j
                  in zip(u.indices[1:], range(len(model.shape)))])
+    # Operator
     return Operator(stencils, subs=subs, name='ForwardTTI', **kwargs)
+
+
+kernels = {('shifted', 3): kernel_shited_3d, ('shifted', 2): kernel_shited_2d,
+           ('centered', 3): kernel_centered_3d, ('centered', 2): kernel_centered_2d}
diff --git a/examples/seismic/tti/tti_example.py b/examples/seismic/tti/tti_example.py
index 0bec1ff154..fe5230a825 100644
--- a/examples/seismic/tti/tti_example.py
+++ b/examples/seismic/tti/tti_example.py
@@ -2,23 +2,16 @@
 from argparse import ArgumentParser
 
 from devito.logger import warning
-from examples.seismic import demo_model, GaborSource, Receiver
+from examples.seismic import demo_model, Receiver, RickerSource
 from examples.seismic.tti import AnisotropicWaveSolver
 
 
-def tti_setup(dimensions=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=250.0,
+def tti_setup(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=250.0,
               time_order=2, space_order=4, nbpml=10, **kwargs):
 
     nrec = 101
-
     # Two layer model for true velocity
-    model = demo_model('layers', ratio=3, nbpml=nbpml,
-                       shape=dimensions, spacing=spacing,
-                       epsilon=.4*np.ones(dimensions),
-                       delta=-.1*np.ones(dimensions),
-                       theta=-np.pi/7*np.ones(dimensions),
-                       phi=np.pi/5*np.ones(dimensions))
-
+    model = demo_model('layers-tti', shape=shape, spacing=spacing, nbpml=nbpml)
     # Derive timestepping from model spacing
     dt = model.critical_dt
     t0 = 0.0
@@ -26,7 +19,7 @@ def tti_setup(dimensions=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=250.0,
     time = np.linspace(t0, tn, nt)
 
     # Define source geometry (center of domain, just below surface)
-    src = GaborSource(name='src', ndim=model.dim, f0=0.01, time=time)
+    src = RickerSource(name='src', ndim=model.dim, f0=0.015, time=time)
     src.coordinates.data[0, :] = np.array(model.domain_size) * .5
     src.coordinates.data[0, -1] = model.origin[-1] + 2 * spacing[-1]
 
@@ -40,15 +33,16 @@ def tti_setup(dimensions=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=250.0,
                                  space_order=space_order, **kwargs)
 
 
-def run(dimensions=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=250.0,
-        autotune=False, time_order=2, space_order=4, nbpml=10, **kwargs):
+def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=250.0,
+        autotune=False, time_order=2, space_order=4, nbpml=10,
+        kernel='centered', **kwargs):
 
-    solver = tti_setup(dimensions, spacing, tn, time_order, space_order, nbpml, **kwargs)
+    solver = tti_setup(shape, spacing, tn, time_order, space_order, nbpml, **kwargs)
 
     if space_order % 4 != 0:
         warning('WARNING: TTI requires a space_order that is a multiple of 4!')
 
-    rec, u, v, summary = solver.forward(autotune=autotune)
+    rec, u, v, summary = solver.forward(autotune=autotune, kernel=kernel)
 
     return summary.gflopss, summary.oi, summary.timings, [rec, u, v]
 
@@ -64,24 +58,30 @@ def run(dimensions=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=250.0,
                         type=int, help="Time order of the simulation")
     parser.add_argument("-so", "--space_order", default=4,
                         type=int, help="Space order of the simulation")
-    parser.add_argument("--nbpml", default=10,
+    parser.add_argument("--nbpml", default=40,
                         type=int, help="Number of PML layers around the domain")
-    parser.add_argument("-dse", default='advanced',
-                        type=str, help="DSE backend choice")
-    parser.add_argument("-dle", default='advanced',
-                        type=str, help="DLE backend choice")
+    parser.add_argument("-k", dest="kernel", default='centered',
+                        choices=['centered', 'shifted'],
+                        help="Choice of finite-difference kernel")
+    parser.add_argument("-dse", "-dse", default="advanced",
+                        choices=["noop", "basic", "advanced",
+                                 "speculative", "aggressive"],
+                        help="Devito symbolic engine (DSE) mode")
+    parser.add_argument("-dle", default="advanced",
+                        choices=["noop", "advanced", "speculative"],
+                        help="Devito loop engine (DSE) mode")
     args = parser.parse_args()
 
     # 3D preset parameters
     if args.dim2:
-        dimensions = (150, 150)
-        spacing = (15.0, 15.0)
+        shape = (150, 150)
+        spacing = (10.0, 10.0)
         tn = 750.0
     else:
-        dimensions = (50, 50, 50)
-        spacing = (20.0, 20.0, 20.0)
+        shape = (50, 50, 50)
+        spacing = (10.0, 10.0, 10.0)
         tn = 250.0
 
-    run(dimensions=dimensions, spacing=spacing, nbpml=args.nbpml, tn=tn,
+    run(shape=shape, spacing=spacing, nbpml=args.nbpml, tn=tn,
         space_order=args.space_order, time_order=args.time_order,
-        autotune=args.autotune, dse=args.dse, dle=args.dle)
+        autotune=args.autotune, dse=args.dse, dle=args.dle, kernel=args.kernel)
diff --git a/examples/seismic/tti/wavesolver.py b/examples/seismic/tti/wavesolver.py
index 081497a286..ea807ac805 100644
--- a/examples/seismic/tti/wavesolver.py
+++ b/examples/seismic/tti/wavesolver.py
@@ -1,7 +1,5 @@
 # coding: utf-8
-from cached_property import cached_property
-
-from devito import TimeData
+from devito import TimeData, memoized
 from examples.seismic.tti.operators import ForwardOperator
 from examples.seismic import Receiver
 
@@ -27,33 +25,23 @@ def __init__(self, model, source, receiver,
         self.receiver = receiver
 
         self.time_order = time_order
-        self.space_order = space_order/2
-
-        # Time step can be \sqrt{3}=1.73 bigger with 4th order
+        self.space_order = space_order
         self.dt = self.model.critical_dt
-        if self.time_order == 4:
-            self.dt *= 1.73
 
         # Cache compiler options
         self._kwargs = kwargs
 
-    @cached_property
-    def op_fwd(self):
+    @memoized
+    def op_fwd(self, kernel='shifted', save=False):
         """Cached operator for forward runs with buffered wavefield"""
-        return ForwardOperator(self.model, save=False, source=self.source,
-                               receiver=self.receiver, time_order=self.time_order,
-                               space_order=self.space_order, **self._kwargs)
-
-    @cached_property
-    def op_fwd_save(self):
-        """Cached operator for forward runs with unrolled wavefield"""
-        return ForwardOperator(self.model, save=True, source=self.source,
+        return ForwardOperator(self.model, save=save, source=self.source,
                                receiver=self.receiver, time_order=self.time_order,
-                               space_order=self.space_order, **self._kwargs)
+                               space_order=self.space_order,
+                               kernel=kernel, **self._kwargs)
 
     def forward(self, src=None, rec=None, u=None, v=None, m=None,
                 epsilon=None, delta=None, theta=None, phi=None,
-                save=False, **kwargs):
+                save=False, kernel='centered', **kwargs):
         """
         Forward modelling function that creates the necessary
         data objects for running a forward modelling operator.
@@ -63,9 +51,15 @@ def forward(self, src=None, rec=None, u=None, v=None, m=None,
         :param u: (Optional) Symbol to store the computed wavefield
         :param m: (Optional) Symbol for the time-constant square slowness
         :param save: Option to store the entire (unrolled) wavefield
+        :param kernel: type of discretization, centered or shifted
 
         :returns: Receiver, wavefield and performance summary
         """
+
+        # Space order needs to be halved in the shifted case to have an
+        # overall space_order discretization
+        self.space_order = self.space_order / 2 if kernel == 'shifted' \
+            else self.space_order
         # Source term is read-only, so re-use the default
         if src is None:
             src = self.source
@@ -84,7 +78,7 @@ def forward(self, src=None, rec=None, u=None, v=None, m=None,
         # Create the forward wavefield if not provided
         if v is None:
             v = TimeData(name='v', shape=self.model.shape_domain,
-                         save=save, time_dim=self.source.nt,
+                         save=save, time_dim=self.source.nt if save else None,
                          time_order=self.time_order,
                          space_order=self.space_order,
                          dtype=self.model.dtype)
@@ -101,10 +95,7 @@ def forward(self, src=None, rec=None, u=None, v=None, m=None,
             phi = phi or self.model.phi
 
         # Execute operator and return wavefield and receiver data
-        if save:
-            op = self.op_fwd_save
-        else:
-            op = self.op_fwd
+        op = self.op_fwd(kernel, save)
 
         if len(m.shape) == 2:
             summary = op.apply(src=src, rec=rec, u=u, v=v, m=m, epsilon=epsilon,
diff --git a/examples/seismic/tutorials/test_02_rtm.ipynb b/examples/seismic/tutorials/test_02_rtm.ipynb
index 86ebfe2923..d2a46da873 100644
--- a/examples/seismic/tutorials/test_02_rtm.ipynb
+++ b/examples/seismic/tutorials/test_02_rtm.ipynb
@@ -73,13 +73,13 @@
     "from examples.seismic import demo_model\n",
     "\n",
     "# Enable model presets here:\n",
-    "preset = 'twolayer'  # A simple but cheap model (recommended)\n",
+    "preset = 'twolayer-isotropic'  # A simple but cheap model (recommended)\n",
     "# preset = 'marmousi2d'  # A larger more realistic model\n",
     "\n",
     "# Standard preset with a simple two-layer model\n",
-    "if preset == 'twolayer':\n",
+    "if preset == 'twolayer-isotropic':\n",
     "    def create_model():\n",
-    "        return demo_model('twolayer', origin=(0., 0.), shape=(101, 101),\n",
+    "        return demo_model('twolayer-isotropic', origin=(0., 0.), shape=(101, 101),\n",
     "                          spacing=(10., 10.), nbpml=20)\n",
     "    filter_sigma = (1, 1)\n",
     "    nshots = 21\n",
@@ -90,9 +90,9 @@
     "\n",
     "\n",
     "# A more computationally demanding preset based on the 2D Marmousi model\n",
-    "if preset == 'marmousi2d':\n",
+    "if preset == 'marmousi2d-isotropic':\n",
     "    def create_model():\n",
-    "        return demo_model('marmousi2d', data_path='../../../../opesci-data/')\n",
+    "        return demo_model('marmousi2d-isotropic', data_path='../../../../opesci-data/')\n",
     "    filter_sigma = (6, 6)\n",
     "    nshots = 301  # Need good covergae in shots, one every two grid points\n",
     "    nreceivers = 601  # One recevier every grid point\n",
@@ -549,7 +549,7 @@
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 2.0
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
@@ -560,5 +560,5 @@
   }
  },
  "nbformat": 4,
- "nbformat_minor": 1
-}
+ "nbformat_minor": 0
+}
\ No newline at end of file
diff --git a/examples/seismic/tutorials/test_03_fwi.ipynb b/examples/seismic/tutorials/test_03_fwi.ipynb
index cea90346d5..ee349ae85b 100644
--- a/examples/seismic/tutorials/test_03_fwi.ipynb
+++ b/examples/seismic/tutorials/test_03_fwi.ipynb
@@ -125,10 +125,10 @@
     "spacing = (10., 10.)  # Grid spacing in m. The domain size is now 1km by 1km\n",
     "origin = (0., 0.)  # Need origin to define relative source and receiver locations\n",
     "\n",
-    "model = demo_model('circle', vp=3.0, vp_background=2.5,\n",
+    "model = demo_model('circle-isotropic', vp=3.0, vp_background=2.5,\n",
     "                    origin=origin, shape=shape, spacing=spacing, nbpml=40)\n",
     "\n",
-    "model0 = demo_model('circle', vp=2.5, vp_background=2.5,\n",
+    "model0 = demo_model('circle-isotropic', vp=2.5, vp_background=2.5,\n",
     "                     origin=origin, shape=shape, spacing=spacing, nbpml=40)\n",
     "\n",
     "plot_velocity(model)\n",
@@ -634,7 +634,7 @@
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 2.0
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
@@ -645,5 +645,5 @@
   }
  },
  "nbformat": 4,
- "nbformat_minor": 1
-}
+ "nbformat_minor": 0
+}
\ No newline at end of file
diff --git a/tests/test_adjointA.py b/tests/test_adjointA.py
index a82864f8db..8b98c38d45 100644
--- a/tests/test_adjointA.py
+++ b/tests/test_adjointA.py
@@ -9,12 +9,12 @@
 
 
 presets = {
-    'constant': {'preset': 'constant'},
-    'layers': {'preset': 'layers', 'ratio': 3},
+    'constant': {'preset': 'constant-isotropic'},
+    'layers': {'preset': 'layers-isotropic', 'ratio': 3},
 }
 
 
-@pytest.mark.parametrize('mkey, dimensions, time_order, space_order, nbpml, fix_dim', [
+@pytest.mark.parametrize('mkey, shape, time_order, space_order, nbpml, fix_dim', [
     # 2D tests with varying time and space orders
     ('layers', (60, 70), 2, 4, 10, False), ('layers', (60, 70), 2, 8, 10, False),
     ('layers', (60, 70), 2, 12, 10, False), ('layers', (60, 70), 4, 4, 10, False),
@@ -28,14 +28,14 @@
     # Constant model in 2D and 3D
     ('constant', (60, 70), 2, 8, 14, False), ('constant', (60, 70, 80), 2, 8, 14, False),
 ])
-def test_acoustic(mkey, dimensions, time_order, space_order, nbpml, fix_dim):
+def test_acoustic(mkey, shape, time_order, space_order, nbpml, fix_dim):
     t0 = 0.0  # Start time
     tn = 500.  # Final time
     nrec = 130  # Number of receivers
 
     # Create model from preset
-    model = demo_model(spacing=[15. for _ in dimensions],
-                       shape=dimensions, nbpml=nbpml, **(presets[mkey]))
+    model = demo_model(spacing=[15. for _ in shape],
+                       shape=shape, nbpml=nbpml, **(presets[mkey]))
 
     # Derive timestepping from model spacing
     dt = model.critical_dt * (1.73 if time_order == 4 else 1.0)
diff --git a/tests/test_adjointJ.py b/tests/test_adjointJ.py
index 72b8e2891e..f252cba5b6 100644
--- a/tests/test_adjointJ.py
+++ b/tests/test_adjointJ.py
@@ -8,17 +8,17 @@
 
 
 @pytest.mark.parametrize('space_order', [4, 8, 12])
-@pytest.mark.parametrize('dimensions', [(60, 70), (40, 50, 30)])
-def test_acousticJ(dimensions, space_order):
+@pytest.mark.parametrize('shape', [(60, 70), (40, 50, 30)])
+def test_acousticJ(shape, space_order):
     t0 = 0.0  # Start time
     tn = 500.  # Final time
-    nrec = dimensions[0]  # Number of receivers
+    nrec = shape[0]  # Number of receivers
     nbpml = 10 + space_order / 2
-    spacing = [15. for _ in dimensions]
+    spacing = [15. for _ in shape]
 
     # Create two-layer "true" model from preset with a fault 1/3 way down
-    model = demo_model('layers', ratio=3, vp_top=1.5, vp_bottom=2.5,
-                       spacing=spacing, shape=dimensions, nbpml=nbpml)
+    model = demo_model('layers-isotropic', ratio=3, vp_top=1.5, vp_bottom=2.5,
+                       spacing=spacing, shape=shape, nbpml=nbpml)
 
     # Derive timestepping from model spacing
     dt = model.critical_dt
@@ -40,8 +40,8 @@ def test_acousticJ(dimensions, space_order):
                                 time_order=2, space_order=space_order)
 
     # Create initial model (m0) with a constant velocity throughout
-    model0 = demo_model('layers', ratio=3, vp_top=1.5, vp_bottom=1.5,
-                        spacing=spacing, shape=dimensions, nbpml=nbpml)
+    model0 = demo_model('layers-isotropic', ratio=3, vp_top=1.5, vp_bottom=1.5,
+                        spacing=spacing, shape=shape, nbpml=nbpml)
 
     # Compute the full wavefield u0
     _, u0, _ = solver.forward(save=True, m=model0.m)
diff --git a/tests/test_dse.py b/tests/test_dse.py
index b60c048cee..69a7e7a232 100644
--- a/tests/test_dse.py
+++ b/tests/test_dse.py
@@ -20,7 +20,7 @@
 # Acoustic
 
 def run_acoustic_forward(dse=None):
-    dimensions = (50, 50, 50)
+    shape = (50, 50, 50)
     spacing = (10., 10., 10.)
     nbpml = 10
     nrec = 101
@@ -28,8 +28,8 @@ def run_acoustic_forward(dse=None):
     tn = 250.0
 
     # Create two-layer model from preset
-    model = demo_model(preset='layers', vp_top=3., vp_bottom=4.5,
-                       spacing=spacing, shape=dimensions, nbpml=nbpml)
+    model = demo_model(preset='layers-isotropic', vp_top=3., vp_bottom=4.5,
+                       spacing=spacing, shape=shape, nbpml=nbpml)
 
     # Derive timestepping from model spacing
     dt = model.critical_dt
@@ -67,17 +67,14 @@ def tti_operator(dse=False):
     t0 = 0.0
     tn = 250.
     nbpml = 10
-    dimensions = (50, 50, 50)
+    shape = (50, 50, 50)
     spacing = (20., 20., 20.)
 
     # Two layer model for true velocity
-    model = demo_model('layers', ratio=3, nbpml=nbpml,
-                       shape=dimensions, spacing=spacing,
-                       epsilon=.4*np.ones(dimensions),
-                       delta=-.1*np.ones(dimensions),
-                       theta=-np.pi/7*np.ones(dimensions),
-                       phi=np.pi/5*np.ones(dimensions))
+    model = demo_model('layers-tti', ratio=3, nbpml=nbpml,
+                       shape=shape, spacing=spacing)
 
+    # Derive timestepping from model spacing
     # Derive timestepping from model spacing
     dt = model.critical_dt
     nt = int(1 + (tn-t0) / dt)  # Number of timesteps
@@ -107,10 +104,10 @@ def tti_nodse():
 def test_tti_clusters_to_graph():
     solver = tti_operator()
 
-    nodes = FindNodes(Expression).visit(solver.op_fwd.elemental_functions +
-                                        (solver.op_fwd,))
+    nodes = FindNodes(Expression).visit(solver.op_fwd('centered').elemental_functions +
+                                        (solver.op_fwd('centered'),))
     expressions = [n.expr for n in nodes]
-    stencils = solver.op_fwd._retrieve_stencils(expressions)
+    stencils = solver.op_fwd('centered')._retrieve_stencils(expressions)
     clusters = clusterize(expressions, stencils)
     assert len(clusters) == 3
 
diff --git a/tests/test_gradient.py b/tests/test_gradient.py
index 04859c2ad5..be4066ed3c 100644
--- a/tests/test_gradient.py
+++ b/tests/test_gradient.py
@@ -9,8 +9,8 @@
 
 @pytest.mark.parametrize('space_order', [4])
 @pytest.mark.parametrize('time_order', [2])
-@pytest.mark.parametrize('dimensions', [(70, 80)])
-def test_gradientFWI(dimensions, time_order, space_order):
+@pytest.mark.parametrize('shape', [(70, 80)])
+def test_gradientFWI(shape, time_order, space_order):
     """
     This test ensure that the FWI gradient computed with devito
     satisfies the Taylor expansion property:
@@ -31,8 +31,8 @@ def test_gradientFWI(dimensions, time_order, space_order):
     :param space_order: order of the spacial discretization scheme
     :return: assertion that the Taylor properties are satisfied
     """
-    spacing = tuple(15. for _ in dimensions)
-    wave = setup(dimensions=dimensions, spacing=spacing,
+    spacing = tuple(15. for _ in shape)
+    wave = setup(shape=shape, spacing=spacing,
                  time_order=time_order, space_order=space_order,
                  nbpml=10+space_order/2)
     m0 = smooth10(wave.model.m.data, wave.model.shape_domain)
@@ -76,8 +76,8 @@ def test_gradientFWI(dimensions, time_order, space_order):
 
 @pytest.mark.parametrize('space_order', [4])
 @pytest.mark.parametrize('time_order', [2])
-@pytest.mark.parametrize('dimensions', [(70, 80)])
-def test_gradientJ(dimensions, time_order, space_order):
+@pytest.mark.parametrize('shape', [(70, 80)])
+def test_gradientJ(shape, time_order, space_order):
     """
     This test ensure that the Jacobian computed with devito
     satisfies the Taylor expansion property:
@@ -92,8 +92,8 @@ def test_gradientJ(dimensions, time_order, space_order):
     :param space_order: order of the spacial discretization scheme
     :return: assertion that the Taylor properties are satisfied
     """
-    spacing = tuple(15. for _ in dimensions)
-    wave = setup(dimensions=dimensions, spacing=spacing,
+    spacing = tuple(15. for _ in shape)
+    wave = setup(shape=shape, spacing=spacing,
                  time_order=time_order, space_order=space_order,
                  tn=1000., nbpml=10+space_order/2)
     m0 = smooth10(wave.model.m.data, wave.model.shape_domain)
@@ -128,4 +128,4 @@ def test_gradientJ(dimensions, time_order, space_order):
 
 
 if __name__ == "__main__":
-    test_gradientJ(dimensions=(60, 70), time_order=2, space_order=4)
+    test_gradientJ(shape=(60, 70), time_order=2, space_order=4)
diff --git a/tests/test_tti.py b/tests/test_tti.py
index 3db68cb445..afe285bc46 100644
--- a/tests/test_tti.py
+++ b/tests/test_tti.py
@@ -9,33 +9,33 @@
 from examples.seismic.tti import AnisotropicWaveSolver
 
 
-@pytest.mark.parametrize('dimensions', [(120, 140), (120, 140, 150)])
+@pytest.mark.parametrize('shape', [(120, 140), (120, 140, 150)])
 @pytest.mark.parametrize('space_order', [4, 8])
-def test_tti(dimensions, space_order):
+def test_tti(shape, space_order):
     nbpml = 10
-    ndim = len(dimensions)
-    origin = [0. for _ in dimensions]
-    spacing = [10. for _ in dimensions]
+    ndim = len(shape)
+    origin = [0. for _ in shape]
+    spacing = [10. for _ in shape]
 
     # Source location
     location = np.zeros((1, ndim), dtype=np.float32)
-    location[0, :-1] = [origin[i] + dimensions[i] * spacing[i] * .5
+    location[0, :-1] = [origin[i] + shape[i] * spacing[i] * .5
                         for i in range(ndim-1)]
     location[0, -1] = origin[-1] + 2 * spacing[-1]
     # Receivers locations
-    receiver_coords = np.zeros((dimensions[0], ndim), dtype=np.float32)
+    receiver_coords = np.zeros((shape[0], ndim), dtype=np.float32)
     receiver_coords[:, 0] = np.linspace(0, origin[0] +
-                                        (dimensions[0]-1) * spacing[0],
-                                        num=dimensions[0])
+                                        (shape[0]-1) * spacing[0],
+                                        num=shape[0])
     receiver_coords[:, 1:] = location[0, 1:]
 
     # Two layer model for true velocity
-    model = demo_model('layers', ratio=3, shape=dimensions,
+    model = demo_model('layers-isotropic', ratio=3, shape=shape,
                        spacing=spacing, nbpml=nbpml,
-                       epsilon=np.zeros(dimensions),
-                       delta=np.zeros(dimensions),
-                       theta=np.zeros(dimensions),
-                       phi=np.zeros(dimensions))
+                       epsilon=np.zeros(shape),
+                       delta=np.zeros(shape),
+                       theta=np.zeros(shape),
+                       phi=np.zeros(shape))
 
     # Define seismic data and parameters
     f0 = .010
@@ -61,7 +61,7 @@ def ricker_source(t, f0):
                                   time_order=2, space_order=space_order)
     rec, u1, _ = acoustic.forward(save=False)
 
-    tn = 50.0
+    tn = 100.0
     nt = int(1 + (tn - t0) / dt)
     # Source geometry
     time_series = np.zeros((nt, 1))
@@ -93,7 +93,7 @@ def ricker_source(t, f0):
                       .5 * u_tti.data.reshape(-1) - .5 * v_tti.data.reshape(-1))
     res /= linalg.norm(u.data.reshape(-1))
     log("Difference between acoustic and TTI with all coefficients to 0 %f" % res)
-    assert np.isclose(res, 0.0, atol=1e-1)
+    assert np.isclose(res, 0.0, atol=1e-4)
 
 
 if __name__ == "__main__":