diff --git a/devito/cgen_wrapper.py b/devito/cgen_wrapper.py
index 07ddd6d5a9..39d1d8984c 100644
--- a/devito/cgen_wrapper.py
+++ b/devito/cgen_wrapper.py
@@ -1,5 +1,4 @@
 import ctypes
-
 from cgen import *
 
 
diff --git a/devito/compiler.py b/devito/compiler.py
old mode 100644
new mode 100755
index 888d03ef47..f6419ef371
--- a/devito/compiler.py
+++ b/devito/compiler.py
@@ -2,9 +2,11 @@
 from tempfile import gettempdir
 
 import numpy.ctypeslib as npct
+
 from cgen import Pragma
 from codepy.jit import extension_file_from_string
 from codepy.toolchain import GCCToolchain
+from devito.logger import log
 
 __all__ = ['get_tmp_dir', 'get_compiler_from_env',
            'jit_compile', 'load', 'jit_compile_and_load',
@@ -90,7 +92,7 @@ def __init__(self, *args, **kwargs):
         self.ldflags = ['-shared']
 
         if self.openmp:
-            print "WARNING: Disabling OpenMP because clang does not support it."
+            log("WARNING: Disabling OpenMP because clang does not support it.")
             self.openmp = False
 
 
@@ -130,8 +132,7 @@ def __init__(self, *args, **kwargs):
         if self.openmp:
             self.ldflags += ['-qopenmp']
         else:
-            print "WARNING: Running on Intel MIC without OpenMP is highly discouraged"
-
+            log("WARNING: Running on Intel MIC without OpenMP is highly discouraged")
         self._mic = __import__('pymic')
 
 
@@ -147,7 +148,7 @@ def __init__(self, *args, **kwargs):
         if self.openmp:
             self.ldflags += ['-qopenmp']
         else:
-            print "WARNING: Running on Intel KNL without OpenMP is highly discouraged"
+            log("WARNING: Running on Intel KNL without OpenMP is highly discouraged")
 
 # Registry dict for deriving Compiler classes according to
 # environment variable DEVITO_ARCH. Developers should add
@@ -206,7 +207,7 @@ def jit_compile(ccode, basename, compiler=GNUCompiler):
     """
     src_file = "%s.cpp" % basename
     lib_file = "%s.so" % basename
-    print "%s: Compiling %s" % (compiler, src_file)
+    log("%s: Compiling %s" % (compiler, src_file))
     extension_file_from_string(toolchain=compiler, ext_file=lib_file,
                                source_string=ccode, source_name=src_file)
 
diff --git a/devito/iteration.py b/devito/iteration.py
index cbf0481a27..581148c913 100644
--- a/devito/iteration.py
+++ b/devito/iteration.py
@@ -1,7 +1,6 @@
 from collections import Iterable
 
 import cgen
-
 from devito.codeprinter import ccode
 
 __all__ = ['Iteration']
diff --git a/devito/operator.py b/devito/operator.py
index 5295a5e293..f210e9064e 100644
--- a/devito/operator.py
+++ b/devito/operator.py
@@ -208,7 +208,7 @@ def apply(self, debug=False):
         args = [param.data for param in self.signature]
         self.propagator.run(args)
 
-        return tuple([param.data for param in self.output_params])
+        return tuple([param for param in self.output_params])
 
     def apply_python(self):
         """Uses Python to apply the operator
diff --git a/devito/propagator.py b/devito/propagator.py
old mode 100644
new mode 100755
diff --git a/devito/tools.py b/devito/tools.py
index 174bf9badb..80eba2e0df 100644
--- a/devito/tools.py
+++ b/devito/tools.py
@@ -1,5 +1,4 @@
 import ctypes
-
 import numpy as np
 from sympy import symbols
 
diff --git a/examples/AcousticWave2D.py b/examples/AcousticWave2D.py
index 08e25430ef..0ac6523287 100644
--- a/examples/AcousticWave2D.py
+++ b/examples/AcousticWave2D.py
@@ -434,17 +434,3 @@ def Born(self, dm):
             U1, U2, U3 = U2, U3, U1
 
         return rec
-
-    def run(self):
-        nt = self.data.nsamples
-
-        print('Starting forward')
-        rec, u = self.Forward(nt)
-
-        res = rec - np.transpose(self.data.traces)
-        f = 0.5*linalg.norm(res)**2
-
-        print('Residual is ', f, 'starting gradient')
-        g = self.Gradient(nt, res, u)
-
-        return f, g[40:-40, 40:-40]
diff --git a/examples/Acoustic_codegen.py b/examples/Acoustic_codegen.py
old mode 100644
new mode 100755
index 60fcb722c2..651c8252a6
--- a/examples/Acoustic_codegen.py
+++ b/examples/Acoustic_codegen.py
@@ -11,41 +11,27 @@ class Acoustic_cg:
     """ Class to setup the problem for the Acoustic Wave
         Note: s_order must always be greater than t_order
     """
-    def __init__(self, model, data, dm_initializer=None, source=None, nbpml=40, t_order=2, s_order=2):
+    def __init__(self, model, data, source=None, nbpml=40, t_order=2, s_order=2):
         self.model = model
         self.t_order = t_order
         self.s_order = s_order
         self.data = data
-        self.dtype = np.float32
+        self.dtype = np.float64
         self.dt = model.get_critical_dt()
-        self.h = model.get_spacing()
-        self.nbpml = nbpml
-        dimensions = self.model.get_shape()
-        pad_list = []
-
-        for dim_index in range(len(dimensions)):
-            pad_list.append((nbpml, nbpml))
-
-        self.model.vp = np.pad(self.model.vp, tuple(pad_list), 'edge')
-        self.data.reinterpolate(self.dt)
-        self.nrec, self.nt = self.data.traces.shape
+        self.model.nbpml = nbpml
         self.model.set_origin(nbpml)
-        self.dm_initializer = dm_initializer
-
         if source is not None:
             self.source = source.read()
             self.source.reinterpolate(self.dt)
             source_time = self.source.traces[0, :]
-
             while len(source_time) < self.data.nsamples:
                 source_time = np.append(source_time, [0.0])
-
             self.data.set_source(source_time, self.dt, self.data.source_coords)
 
         def damp_boundary(damp):
-            h = self.h
+            h = self.model.get_spacing()
             dampcoeff = 1.5 * np.log(1.0 / 0.001) / (40 * h)
-            nbpml = self.nbpml
+            nbpml = self.model.nbpml
             num_dim = len(damp.shape)
 
             for i in range(nbpml):
@@ -64,60 +50,42 @@ def damp_boundary(damp):
                     damp[:, -(i + 1), :] += val
                     damp[:, :, i] += val
                     damp[:, :, -(i + 1)] += val
-        self.m = DenseData(name="m", shape=self.model.vp.shape, dtype=self.dtype)
-        self.m.data[:] = self.model.vp**(-2)
-
-        self.damp = DenseData(name="damp", shape=self.model.vp.shape, dtype=self.dtype)
+        self.damp = DenseData(name="damp", shape=self.model.get_shape_comp(),
+                              dtype=self.dtype)
 
         # Initialize damp by calling the function that can precompute damping
         damp_boundary(self.damp.data)
-
-        self.src = SourceLike(name="src", npoint=1, nt=self.nt, dt=self.dt, h=self.h,
-                              coordinates=np.array(self.data.source_coords, dtype=self.dtype)[np.newaxis, :],
-                              ndim=len(dimensions), dtype=self.dtype, nbpml=nbpml)
-        self.rec = SourceLike(name="rec", npoint=self.nrec, nt=self.nt, dt=self.dt, h=self.h,
-                              coordinates=self.data.receiver_coords, ndim=len(dimensions), dtype=self.dtype,
-                              nbpml=nbpml)
-        self.src.data[:] = self.data.get_source()[:, np.newaxis]
-
-        self.u = TimeData(name="u", shape=self.m.shape, time_dim=self.src.nt, time_order=t_order,
-                          save=True, dtype=self.m.dtype)
-        self.srca = SourceLike(name="srca", npoint=1, nt=self.nt, dt=self.dt, h=self.h,
-                               coordinates=np.array(self.data.source_coords, dtype=self.dtype)[np.newaxis, :],
-                               ndim=len(dimensions), dtype=self.dtype, nbpml=nbpml)
-        if dm_initializer is not None:
-            self.dm = DenseData(name="dm", shape=self.model.vp.shape, dtype=self.dtype)
-            self.dm.data[:] = np.pad(dm_initializer, tuple(pad_list), 'edge')
-
-    def Forward(self):
-        fw = ForwardOperator(
-            self.m, self.src, self.damp, self.rec, self.u, time_order=self.t_order, spc_order=self.s_order)
-        fw.apply()
-
-        return (self.rec.data, self.u.data)
+        srccoord = np.array(self.data.source_coords, dtype=self.dtype)[np.newaxis, :]
+        self.src = SourceLike(name="src", npoint=1, nt=data.traces.shape[1],
+                              dt=self.dt, h=self.model.get_spacing(),
+                              coordinates=srccoord, ndim=len(self.damp.shape),
+                              dtype=self.dtype, nbpml=nbpml)
+        self.src.data[:] = data.get_source()[:, np.newaxis]
+
+    def Forward(self, save=False):
+        fw = ForwardOperator(self.model, self.src, self.damp, self.data,
+                             time_order=self.t_order, spc_order=self.s_order,
+                             save=save)
+        u, rec = fw.apply()
+        return rec.data, u
 
     def Adjoint(self, rec):
-        self.rec.data[:] = rec
-        adj = AdjointOperator(self.m, self.rec, self.damp, self.srca, time_order=self.t_order, spc_order=self.s_order)
+        adj = AdjointOperator(self.model, self.damp, self.data, rec,
+                              time_order=self.t_order, spc_order=self.s_order)
         v = adj.apply()[0]
-
-        return (self.srca.data, v)
+        return v.data
 
     def Gradient(self, rec, u):
-        self.u.data[:] = u
-        self.rec.data[:] = rec
-        grad_op = GradientOperator(self.u, self.m, self.rec, self.damp, time_order=self.t_order, spc_order=self.s_order)
-        dt = self.dt
+        grad_op = GradientOperator(self.model, self.damp, self.data, rec, u,
+                                   time_order=self.t_order, spc_order=self.s_order)
         grad = grad_op.apply()[0]
+        return grad.data
 
-        return (dt**-2)*grad
-
-    def Born(self):
-        born_op = BornOperator(
-            self.dm, self.m, self.src, self.damp, self.rec, time_order=self.t_order, spc_order=self.s_order)
-        born_op.apply()
-
-        return self.rec.data
+    def Born(self, dm):
+        born_op = BornOperator(self.model, self.src, self.damp, self.data, dm,
+                               time_order=self.t_order, spc_order=self.s_order)
+        rec = born_op.apply()[0]
+        return rec.data
 
     def run(self):
         print('Starting forward')
diff --git a/examples/TTI_codegen.py b/examples/TTI_codegen.py
index 3afa327715..3fe31af810 100755
--- a/examples/TTI_codegen.py
+++ b/examples/TTI_codegen.py
@@ -3,7 +3,6 @@
 
 import numpy as np
 
-from devito.interfaces import DenseData
 from examples.tti_operators import *
 
 
@@ -11,28 +10,17 @@ class TTI_cg:
     """ Class to setup the problem for the Acoustic Wave
         Note: s_order must always be greater than t_order
     """
-    def __init__(self, model, data, dm_initializer=None, source=None, t_order=2, s_order=2, nbpml=40):
+    def __init__(self, model, data, source=None, t_order=2, s_order=2, nbpml=40, save=False):
         self.model = model
         self.t_order = t_order
         self.s_order = s_order
         self.data = data
-        self.dtype = np.float64
-        self.dt = self.model.get_critical_dt()
-        self.h = model.get_spacing()
-        self.nbpml = nbpml
-        dimensions = self.model.get_shape()
-        pad_list = []
-        for dim_index in range(len(dimensions)):
-            pad_list.append((nbpml, nbpml))
-        self.model.vp = np.pad(self.model.vp, tuple(pad_list), 'edge')
-        self.model.epsilon = np.pad(self.model.epsilon, tuple(pad_list), 'edge')
-        self.model.delta = np.pad(self.model.delta, tuple(pad_list), 'edge')
-        self.model.theta = np.pad(self.model.theta, tuple(pad_list), 'edge')
+        self.dtype = np.float32
+        self.dt = model.get_critical_dt()
+        self.model.nbpml = nbpml
         self.model.set_origin(nbpml)
         self.data.reinterpolate(self.dt)
-        self.nrec, self.nt = self.data.traces.shape
-        self.model.set_origin(nbpml)
-        self.dm_initializer = dm_initializer
+
         if source is not None:
             self.source = source.read()
             self.source.reinterpolate(self.dt)
@@ -42,9 +30,9 @@ def __init__(self, model, data, dm_initializer=None, source=None, t_order=2, s_o
             self.data.set_source(source_time, self.dt, self.data.source_coords)
 
         def damp_boundary(damp):
-            h = self.h
+            h = self.model.get_spacing()
             dampcoeff = 1.5 * np.log(1.0 / 0.001) / (40 * h)
-            nbpml = self.nbpml
+            nbpml = self.model.nbpml
             num_dim = len(damp.shape)
             for i in range(nbpml):
                 pos = np.abs((nbpml-i)/float(nbpml))
@@ -62,70 +50,26 @@ def damp_boundary(damp):
                     damp[:, :, i] += val
                     damp[:, :, -(i + 1)] += val
 
-        self.m = DenseData(name="m", shape=self.model.vp.shape, dtype=self.dtype)
-        self.m.data[:] = self.model.vp**(-2)
-        self.a = DenseData(name="a", shape=self.model.vp.shape, dtype=self.dtype)
-        self.a.data[:] = 1.0 + 2.0 * self.model.epsilon
-        self.b = DenseData(name="b", shape=self.model.vp.shape, dtype=self.dtype)
-        self.b.data[:] = np.sqrt(1.0 + 2.0 * self.model.delta)
-        self.th = DenseData(name="th", shape=self.model.vp.shape, dtype=self.dtype)
-        self.th.data[:] = self.model.theta
-        if len(dimensions) == 2:
-            self.ph = None
-        else:
-            self.model.phi = np.pad(self.model.phi, tuple(pad_list), 'edge')
-            self.ph = DenseData(name="ph", shape=self.model.vp.shape, dtype=self.dtype)
-            self.ph.data[:] = self.model.phi
-        self.damp = DenseData(name="damp", shape=self.model.vp.shape, dtype=self.dtype)
+        self.damp = DenseData(name="damp", shape=self.model.get_shape_comp(),
+                              dtype=self.dtype)
         # Initialize damp by calling the function that can precompute damping
         damp_boundary(self.damp.data)
-        self.src = SourceLikeTTI(name="src", npoint=1, nt=self.nt, dt=self.dt, h=self.h,
-                                 coordinates=np.array(self.data.source_coords, dtype=self.dtype)[np.newaxis, :],
-                                 ndim=len(dimensions), dtype=self.dtype, nbpml=nbpml)
-        self.rec = SourceLikeTTI(name="rec", npoint=self.nrec, nt=self.nt, dt=self.dt, h=self.h,
-                                 coordinates=self.data.receiver_coords, ndim=len(dimensions), dtype=self.dtype,
-                                 nbpml=nbpml)
-        self.src.data[:] = self.data.get_source()[:, np.newaxis]
-        self.u = TimeData(name="u", shape=self.m.shape, time_dim=self.src.nt, time_order=t_order,
-                          save=False, dtype=self.m.dtype)
-        self.v = TimeData(name="v", shape=self.m.shape, time_dim=self.src.nt, time_order=t_order,
-                          save=False, dtype=self.m.dtype)
-        self.srca = SourceLikeTTI(name="srca", npoint=1, nt=self.nt, dt=self.dt, h=self.h,
-                                  coordinates=np.array(self.data.source_coords, dtype=self.dtype)[np.newaxis, :],
-                                  ndim=len(dimensions), dtype=self.dtype, nbpml=nbpml)
-        if dm_initializer is not None:
-            self.dm = DenseData(name="dm", shape=self.model.vp.shape, dtype=self.dtype)
-            self.dm.data[:] = np.pad(dm_initializer, tuple(pad_list), 'edge')
+        srccoord = np.array(self.data.source_coords, dtype=self.dtype)[np.newaxis, :]
+        self.src = SourceLike(name="src", npoint=1, nt=data.traces.shape[1],
+                              dt=self.dt, h=self.model.get_spacing(),
+                              coordinates=srccoord, ndim=len(self.damp.shape),
+                              dtype=self.dtype, nbpml=nbpml)
+        self.src.data[:] = data.get_source()[:, np.newaxis]
 
-    def Forward(self):
-        fw = ForwardOperator(self.m, self.src, self.damp, self.rec, self.u, self.v, self.a,
-                             self.b, self.th, self.ph, time_order=self.t_order, spc_order=self.s_order)
-        fw.apply()
-        return (self.rec.data, self.u.data, self.v.data)
+    def Forward(self, save=False):
+        fw = ForwardOperator(self.model, self.src, self.damp, self.data,
+                             time_order=self.t_order, spc_order=self.s_order,
+                             save=save)
+        u, v, rec = fw.apply()
+        return (rec.data, u.data, v.data)
 
     def Adjoint(self, rec):
-        adj = AdjointOperator(self.m, self.rec, self.damp, self.srca, self.u, self.v, self.a,
-                              self.b, self.th, self.ph, time_order=self.t_order, spc_order=self.s_order)
-        adj.apply()
-        return (self.srca.data, self.u.data, self.v.data)
-
-    def Gradient(self, rec, u):
-        grad_op = GradientOperator(self.u, self.m, self.rec, self.damp, time_order=self.t_order, spc_order=self.s_order)
-        dt = self.dt
-        grad = grad_op.apply()[0]
-        return (dt**-2)*grad
-
-    def Born(self):
-        born_op = BornOperator(
-            self.dm, self.m, self.src, self.damp, self.rec, time_order=self.t_order, spc_order=self.s_order)
-        born_op.apply()
-        return self.rec.data
-
-    def run(self):
-        print('Starting forward')
-        rec, u = self.Forward()
-        res = rec - np.transpose(self.data.traces)
-        f = 0.5*np.linalg.norm(res)**2
-        print('Residual is ', f, 'starting gradient')
-        g = self.Gradient(res, u)
-        return f, g[self.nbpml:-self.nbpml, self.nbpml:-self.nbpml]
+        adj = AdjointOperator(self.model, self.damp, self.data, rec,
+                              time_order=self.t_order, spc_order=self.s_order)
+        srca = adj.apply()[0]
+        return srca.data
diff --git a/examples/acoustic_example.py b/examples/acoustic_example.py
index 2c9485d06f..e7f2f3d5eb 100644
--- a/examples/acoustic_example.py
+++ b/examples/acoustic_example.py
@@ -36,11 +36,7 @@ def smooth10(vel, shape):
 # Smooth velocity
 initial_vp = smooth10(true_vp, dimensions)
 
-dm_orig = true_vp**-2 - initial_vp**-2
-
-
-def create_dm(dm):
-    np.copyto(dm, dm_orig)
+dm = true_vp**-2 - initial_vp**-2
 
 dv = -true_vp + initial_vp
 
@@ -55,7 +51,6 @@ def create_dm(dm):
 tn = 250.0
 nt = int(1+(tn-t0)/dt)
 h = model.get_spacing()
-data.reinterpolate(dt)
 
 
 # Set up the source as Ricker wavelet for f0
@@ -76,15 +71,11 @@ def source(t, f0):
 receiver_coords[:, 2] = location[2]
 data.set_receiver_pos(receiver_coords)
 data.set_shape(nt, 101)
-
-Acoustic = Acoustic_cg(model, data, dm_orig, None, nbpml=10, t_order=2, s_order=2)
-print("Preparing Forward")
-print("Applying")
-(rec, u) = Acoustic.Forward()
-
+Acoustic = Acoustic_cg(model, data)
+(rec, u) = Acoustic.Forward(save=True)
 print("Preparing adjoint")
 print("Applying")
-(srca, v) = Acoustic.Adjoint(rec)
+srca = Acoustic.Adjoint(rec)
 
 print("Preparing Gradient")
 print("Applying")
@@ -92,4 +83,4 @@ def source(t, f0):
 
 print("Preparing Born")
 print("Applying")
-LinRec = Acoustic.Born()
+LinRec = Acoustic.Born(dm)
diff --git a/examples/containers.py b/examples/containers.py
old mode 100644
new mode 100755
index 86b09bccf0..823594b8d0
--- a/examples/containers.py
+++ b/examples/containers.py
@@ -37,17 +37,25 @@ def get_spacing(self):
 
     def create_model(self, origin, spacing, vp, epsilon=None, delta=None, theta=None, phi=None):
         self.vp = vp
-        self.epsilon = epsilon
-        self.delta = delta
-        self.theta = theta
-        self.phi = phi
         self.spacing = spacing
+        self.dimensions = vp.shape
         if epsilon is not None:
+            self.epsilon = 1 + 2 * epsilon
             self.scale = np.sqrt(1 + 2 * np.max(self.epsilon))
         else:
             self.scale = 1
+        if delta is not None:
+            self.delta = np.sqrt(1 + 2 * delta)
+        if phi is not None:
+            self.theta = theta
+        if theta is not None:
+            self.phi = phi
+
         self.origin = origin
 
+    def set_vp(self, vp):
+        self.vp = vp
+
     def set_origin(self, shift):
         norig = len(self.origin)
         aux = []
@@ -60,6 +68,22 @@ def set_origin(self, shift):
     def get_origin(self):
         return self.origin
 
+    def padm(self):
+        return self.pad(self.vp**(-2))
+
+    def pad(self, m):
+        pad_list = []
+        for dim_index in range(len(self.vp.shape)):
+            pad_list.append((self.nbpml, self.nbpml))
+        return np.pad(m, pad_list, 'edge')
+
+    def get_shape_comp(self):
+        dim = self.dimensions
+        if len(dim) == 3:
+            return (dim[0] + 2 * self.nbpml, dim[1] + 2 * self.nbpml, dim[2] + 2 * self.nbpml)
+        else:
+            return (dim[0] + 2 * self.nbpml, dim[1] + 2 * self.nbpml)
+
 
 class ISource:
     def get_source(self):
diff --git a/examples/fwi_operators.py b/examples/fwi_operators.py
old mode 100644
new mode 100755
index 7e9972c067..c7b096cc43
--- a/examples/fwi_operators.py
+++ b/examples/fwi_operators.py
@@ -1,3 +1,4 @@
+import numpy as np
 from sympy import Eq, Function, Matrix, solve, symbols
 from sympy.abc import p, t
 
@@ -79,7 +80,7 @@ def sym_coord_bases(self):
                       for x, idx in zip(self.sym_coordinates,
                                         self.sym_coord_indices)])
 
-    def point2grid(self, u, m):
+    def point2grid(self, u, m, t=t):
         """Generates an expression for generic point-to-grid interpolation"""
         dt = self.dt
         subs = dict(zip(self.rs, self.sym_coord_bases))
@@ -91,7 +92,7 @@ def point2grid(self, u, m):
                 for idx, b in zip(index_matrix, self.bs)]
         return eqns
 
-    def grid2point(self, u):
+    def grid2point(self, u, t=t):
         """Generates an expression for generic grid-to-point interpolation"""
         subs = dict(zip(self.rs, self.sym_coord_bases))
         index_matrix = [tuple([idx + ii + self.nbpml for ii, idx
@@ -105,35 +106,47 @@ def read(self, u):
         interp_expr = Eq(self.indexed[t, p], self.grid2point(u))
         return [Iteration(interp_expr, variable=p, limits=self.shape[1])]
 
-    def add(self, m, u):
+    def read2(self, u, v):
+        """Iteration loop over points performing grid-to-point interpolation."""
+        interp_expr = Eq(self.indexed[t, p], self.grid2point(u) + self.grid2point(v))
+        return [Iteration(interp_expr, variable=p, limits=self.shape[1])]
+
+    def add(self, m, u, t=t):
         """Iteration loop over points performing point-to-grid interpolation."""
-        return [Iteration(self.point2grid(u, m), variable=p, limits=self.shape[1])]
+        return [Iteration(self.point2grid(u, m, t), variable=p, limits=self.shape[1])]
 
 
 class ForwardOperator(Operator):
-    def __init__(self, m, src, damp, rec, u, time_order=2, spc_order=6, **kwargs):
-        assert(m.shape == damp.shape)
-
-        u.pad_time = True
-        # Set time and space orders
-        u.time_order = time_order
-        u.space_order = spc_order
-
+    def __init__(self, model, src, damp, data, time_order=2, spc_order=6, save=False, **kwargs):
+        nrec, nt = data.traces.shape
+        dt = model.get_critical_dt()
+        u = TimeData(name="u", shape=model.get_shape_comp(), time_dim=nt,
+                     time_order=time_order, space_order=spc_order, save=save,
+                     dtype=damp.dtype)
+        m = DenseData(name="m", shape=model.get_shape_comp(), dtype=damp.dtype)
+        m.data[:] = model.padm()
+        u.pad_time = save
+        rec = SourceLike(name="rec", npoint=nrec, nt=nt, dt=dt, h=model.get_spacing(),
+                         coordinates=data.receiver_coords, ndim=len(damp.shape),
+                         dtype=damp.dtype, nbpml=model.nbpml)
         # Derive stencil from symbolic equation
         eqn = m * u.dt2 - u.laplace + damp * u.dt
         stencil = solve(eqn, u.forward)[0]
-
         # Add substitutions for spacing (temporal and spatial)
         s, h = symbols('s h')
-        subs = {s: src.dt, h: src.h}
-
-        super(ForwardOperator, self).__init__(src.nt, m.shape, stencils=Eq(u.forward, stencil),
-                                              substitutions=subs, spc_border=spc_order/2,
-                                              time_order=time_order, forward=True, dtype=m.dtype,
+        subs = {s: dt, h: model.get_spacing()}
+        super(ForwardOperator, self).__init__(nt, m.shape,
+                                              stencils=Eq(u.forward, stencil),
+                                              substitutions=subs,
+                                              spc_border=spc_order/2,
+                                              time_order=time_order,
+                                              forward=True,
+                                              dtype=m.dtype,
                                               **kwargs)
 
         # Insert source and receiver terms post-hoc
         self.input_params += [src, src.coordinates, rec, rec.coordinates]
+        self.output_params += [rec]
         self.propagator.time_loop_stencils_a = src.add(m, u) + rec.read(u)
         self.propagator.add_devito_param(src)
         self.propagator.add_devito_param(src.coordinates)
@@ -142,27 +155,42 @@ def __init__(self, m, src, damp, rec, u, time_order=2, spc_order=6, **kwargs):
 
 
 class AdjointOperator(Operator):
-    def __init__(self, m, rec, damp, srca, time_order=2, spc_order=6, **kwargs):
-        assert(m.shape == damp.shape)
-
-        # Create v with given time and space orders
-        v = TimeData(name="v", shape=m.shape, dtype=m.dtype, time_dim=rec.nt,
-                     time_order=time_order, space_order=spc_order, save=False)
+    def __init__(self, model, damp, data, recin, time_order=2, spc_order=6, **kwargs):
+        nrec, nt = data.traces.shape
+        dt = model.get_critical_dt()
+        v = TimeData(name="v", shape=model.get_shape_comp(), time_dim=nt,
+                     time_order=time_order, space_order=spc_order,
+                     save=False, dtype=damp.dtype)
+        m = DenseData(name="m", shape=model.get_shape_comp(), dtype=damp.dtype)
+        m.data[:] = model.padm()
+        v.pad_time = False
+        srca = SourceLike(name="srca", npoint=1, nt=nt, dt=dt, h=model.get_spacing(),
+                          coordinates=np.array(data.source_coords, dtype=damp.dtype)[np.newaxis, :],
+                          ndim=len(damp.shape), dtype=damp.dtype, nbpml=model.nbpml)
+        rec = SourceLike(name="rec", npoint=nrec, nt=nt, dt=dt, h=model.get_spacing(),
+                         coordinates=data.receiver_coords, ndim=len(damp.shape),
+                         dtype=damp.dtype, nbpml=model.nbpml)
+        rec.data[:] = recin[:]
         # Derive stencil from symbolic equation
         eqn = m * v.dt2 - v.laplace - damp * v.dt
         stencil = solve(eqn, v.backward)[0]
 
         # Add substitutions for spacing (temporal and spatial)
         s, h = symbols('s h')
-        subs = {s: rec.dt, h: rec.h}
-        super(AdjointOperator, self).__init__(rec.nt, m.shape, stencils=Eq(v.backward, stencil),
-                                              substitutions=subs, spc_border=spc_order/2,
-                                              time_order=time_order, forward=False, dtype=m.dtype,
+        subs = {s: model.get_critical_dt(), h: model.get_spacing()}
+        super(AdjointOperator, self).__init__(nt, m.shape,
+                                              stencils=Eq(v.backward, stencil),
+                                              substitutions=subs,
+                                              spc_border=spc_order/2,
+                                              time_order=time_order,
+                                              forward=False,
+                                              dtype=m.dtype,
                                               **kwargs)
 
         # Insert source and receiver terms post-hoc
         self.input_params += [srca, srca.coordinates, rec, rec.coordinates]
         self.propagator.time_loop_stencils_a = rec.add(m, v) + srca.read(v)
+        self.output_params = [srca]
         self.propagator.add_devito_param(srca)
         self.propagator.add_devito_param(srca.coordinates)
         self.propagator.add_devito_param(rec)
@@ -170,11 +198,19 @@ def __init__(self, m, rec, damp, srca, time_order=2, spc_order=6, **kwargs):
 
 
 class GradientOperator(Operator):
-    def __init__(self, u, m, rec, damp, time_order=2, spc_order=6, **kwargs):
-        assert(m.shape == damp.shape)
-
-        v = TimeData(name="v", shape=m.shape, dtype=m.dtype, time_dim=rec.nt,
-                     time_order=time_order, space_order=spc_order, save=False, )
+    def __init__(self, model, damp, data, recin, u, time_order=2, spc_order=6, **kwargs):
+        nrec, nt = data.traces.shape
+        dt = model.get_critical_dt()
+        v = TimeData(name="v", shape=model.get_shape_comp(), time_dim=nt,
+                     time_order=time_order, space_order=spc_order,
+                     save=False, dtype=damp.dtype)
+        m = DenseData(name="m", shape=model.get_shape_comp(), dtype=damp.dtype)
+        m.data[:] = model.padm()
+        v.pad_time = False
+        rec = SourceLike(name="rec", npoint=nrec, nt=nt, dt=dt, h=model.get_spacing(),
+                         coordinates=data.receiver_coords, ndim=len(damp.shape),
+                         dtype=damp.dtype, nbpml=model.nbpml)
+        rec.data[:] = recin
         grad = DenseData(name="grad", shape=m.shape, dtype=m.dtype)
 
         # Derive stencil from symbolic equation
@@ -183,35 +219,46 @@ def __init__(self, u, m, rec, damp, time_order=2, spc_order=6, **kwargs):
 
         # Add substitutions for spacing (temporal and spatial)
         s, h = symbols('s h')
-        subs = {s: rec.dt, h: rec.h}
-
-        # Add Gradient-specific updates and resets
-        gradient_update = Eq(grad, grad - (v - 2 * v.forward + v.forward.forward) * u)
-        # reset_v = Eq(v.indexed[tuple((t + 2,) + space_dim)], 0)
-        stencils = [Eq(v.backward, stencil), gradient_update]
-        super(GradientOperator, self).__init__(rec.nt, m.shape, stencils=stencils,
-                                               substitutions=[subs, subs, {}], spc_border=spc_order/2,
-                                               time_order=time_order, forward=False, dtype=m.dtype,
+        subs = {s: model.get_critical_dt(), h: model.get_spacing()}
+        # Add Gradient-specific updates. The dt2 is currently hacky as it has to match the cyclic indices
+        gradient_update = Eq(grad, grad - s**-2*(v + v.forward - 2 * v.forward.forward) * u.forward)
+        stencils = [gradient_update, Eq(v.backward, stencil)]
+        super(GradientOperator, self).__init__(rec.nt - 1, m.shape,
+                                               stencils=stencils,
+                                               substitutions=[subs, subs, {}],
+                                               spc_border=spc_order/2,
+                                               time_order=time_order,
+                                               forward=False,
+                                               dtype=m.dtype,
+                                               input_params=[m, v, damp, u],
                                                **kwargs)
-
         # Insert receiver term post-hoc
-        self.input_params += [u, rec, rec.coordinates]
+        self.input_params += [grad, rec, rec.coordinates]
         self.output_params = [grad]
-        self.propagator.time_loop_stencils_a = rec.add(m, v)
-        self.propagator.add_devito_param(u)
+        self.propagator.time_loop_stencils_b = rec.add(m, v, t + 1)
         self.propagator.add_devito_param(rec)
         self.propagator.add_devito_param(rec.coordinates)
-        self.propagator.add_devito_param(grad)
 
 
 class BornOperator(Operator):
-    def __init__(self, dm, m, src, damp, rec, time_order=2, spc_order=6, **kwargs):
-        assert(m.shape == damp.shape)
-
-        u = TimeData(name="u", shape=m.shape, dtype=m.dtype, time_dim=src.nt,
-                     time_order=time_order, space_order=spc_order, save=False)
-        U = TimeData(name="U", shape=m.shape, dtype=m.dtype, time_dim=src.nt,
-                     time_order=time_order, space_order=spc_order, save=False)
+    def __init__(self, model, src, damp, data, dmin, time_order=2, spc_order=6, **kwargs):
+        nrec, nt = data.traces.shape
+        dt = model.get_critical_dt()
+        u = TimeData(name="u", shape=model.get_shape_comp(), time_dim=nt,
+                     time_order=time_order, space_order=spc_order,
+                     save=False, dtype=damp.dtype)
+        U = TimeData(name="U", shape=model.get_shape_comp(), time_dim=nt,
+                     time_order=time_order, space_order=spc_order,
+                     save=False, dtype=damp.dtype)
+        m = DenseData(name="m", shape=model.get_shape_comp(), dtype=damp.dtype)
+        m.data[:] = model.padm()
+
+        dm = DenseData(name="dm", shape=model.get_shape_comp(), dtype=damp.dtype)
+        dm.data[:] = model.pad(dmin)
+
+        rec = SourceLike(name="rec", npoint=nrec, nt=nt, dt=dt, h=model.get_spacing(),
+                         coordinates=data.receiver_coords, ndim=len(damp.shape),
+                         dtype=damp.dtype, nbpml=model.nbpml)
 
         # Derive stencils from symbolic equation
         first_eqn = m * u.dt2 - u.laplace - damp * u.dt
@@ -224,24 +271,23 @@ def __init__(self, dm, m, src, damp, rec, time_order=2, spc_order=6, **kwargs):
         subs = {s: src.dt, h: src.h}
 
         # Add Born-specific updates and resets
-        total_dim = tuple(u.indices(m.shape))
-        space_dim = tuple(m.indices(m.shape))
-        src2 = -(src.dt**-2) * (u.indexed[total_dim]-2*u.indexed[tuple((t - 1,) + space_dim)]
-                                + u.indexed[tuple((t - 2,) + space_dim)]) * dm.indexed[space_dim]
-        insert_second_source = Eq(U.indexed[total_dim], U.indexed[total_dim] + (src.dt * src.dt)
-                                  / m.indexed[space_dim]*src2)
-        reset_u = Eq(u.indexed[tuple((t - 2,) + space_dim)], 0)
+        src2 = -(src.dt**-2) * (- 2 * u + u.forward + u.backward) * dm
+        insert_second_source = Eq(U, U + (src.dt * src.dt) / m*src2)
         stencils = [Eq(u.forward, first_stencil), Eq(U.forward, second_stencil),
-                    insert_second_source, reset_u]
-        super(BornOperator, self).__init__(src.nt, m.shape, stencils=stencils,
-                                           substitutions=[subs, subs, {}, {}], spc_border=spc_order/2,
-                                           time_order=time_order, forward=True, dtype=m.dtype,
+                    insert_second_source]
+        super(BornOperator, self).__init__(src.nt, m.shape,
+                                           stencils=stencils,
+                                           substitutions=[subs, subs, {}, {}],
+                                           spc_border=spc_order/2,
+                                           time_order=time_order,
+                                           forward=True,
+                                           dtype=m.dtype,
                                            **kwargs)
 
         # Insert source and receiver terms post-hoc
-        self.input_params += [dm, src, src.coordinates, rec, rec.coordinates]
-        self.output_params += [U]
-        self.propagator.time_loop_stencils_b = src.add(m, u)
+        self.input_params += [dm, src, src.coordinates, rec, rec.coordinates, U]
+        self.output_params = [rec]
+        self.propagator.time_loop_stencils_b = src.add(m, u, t - 1)
         self.propagator.time_loop_stencils_a = rec.read(U)
         self.propagator.add_devito_param(dm)
         self.propagator.add_devito_param(src)
diff --git a/examples/tti_example.py b/examples/tti_example.py
old mode 100644
new mode 100755
index f2778c5431..1e72536b28
--- a/examples/tti_example.py
+++ b/examples/tti_example.py
@@ -50,7 +50,7 @@ def source(t, f0):
 data.set_receiver_pos(receiver_coords)
 data.set_shape(nt, 101)
 
-TTI = TTI_cg(model, data, None, None, t_order=2, s_order=2, nbpml=10)
+TTI = TTI_cg(model, data, None, t_order=2, s_order=2, nbpml=10)
 (rec, u, v) = TTI.Forward()
 
 # recf = open('RecTTI','w')
diff --git a/examples/tti_example2D.py b/examples/tti_example2D.py
old mode 100644
new mode 100755
index c21df26b65..365522d5b7
--- a/examples/tti_example2D.py
+++ b/examples/tti_example2D.py
@@ -44,7 +44,7 @@ def source(t, f0):
 data.set_receiver_pos(receiver_coords)
 data.set_shape(nt, 301)
 
-TTI = TTI_cg(model, data, None, None, t_order=2, s_order=spc_order, nbpml=10)
+TTI = TTI_cg(model, data, None, t_order=2, s_order=spc_order, nbpml=10)
 (rec, u, v) = TTI.Forward()
 
 # recw = open('RecTTI', 'w')
diff --git a/examples/tti_operators.py b/examples/tti_operators.py
old mode 100644
new mode 100755
index 651b9171f0..a63b155e57
--- a/examples/tti_operators.py
+++ b/examples/tti_operators.py
@@ -2,118 +2,64 @@
 from sympy.abc import *
 
 from devito.finite_difference import first_derivative
-from devito.interfaces import DenseData, PointData, TimeData
-from devito.iteration import Iteration
+from devito.interfaces import DenseData, TimeData
 from devito.operator import *
+from examples.fwi_operators import SourceLike
 
 
-class SourceLikeTTI(PointData):
-    """Defines the behaviour of sources and receivers.
-    """
-    def __init__(self, *args, **kwargs):
-        self.dt = kwargs.get('dt')
-        self.h = kwargs.get('h')
-        self.ndim = kwargs.get('ndim')
-        self.nbpml = kwargs.get('nbpml')
-        PointData.__init__(self, *args, **kwargs)
-        x1, y1, z1, x2, y2, z2 = symbols('x1, y1, z1, x2, y2, z2')
-
-        if self.ndim == 2:
-            A = Matrix([[1, x1, z1, x1*z1],
-                        [1, x1, z2, x1*z2],
-                        [1, x2, z1, x2*z1],
-                        [1, x2, z2, x2*z2]])
-            self.increments = (0, 0), (0, 1), (1, 0), (1, 1)
-            self.rs = symbols('rx, rz')
-            rx, rz = self.rs
-            p = Matrix([[1],
-                        [rx],
-                        [rz],
-                        [rx*rz]])
+class ForwardOperator(Operator):
+    def __init__(self, model, src, damp, data, time_order=2, spc_order=4, save=False, **kwargs):
+        nrec, nt = data.traces.shape
+        dt = model.get_critical_dt()
+        u = TimeData(name="u", shape=model.get_shape_comp(),
+                     time_dim=nt, time_order=time_order,
+                     space_order=spc_order,
+                     save=save, dtype=damp.dtype)
+        v = TimeData(name="v", shape=model.get_shape_comp(),
+                     time_dim=nt, time_order=time_order,
+                     space_order=spc_order,
+                     save=save, dtype=damp.dtype)
+        m = DenseData(name="m", shape=model.get_shape_comp(),
+                      dtype=damp.dtype)
+        m.data[:] = model.padm()
+
+        if model.epsilon is not None:
+            epsilon = DenseData(name="epsilon", shape=model.get_shape_comp(),
+                                dtype=damp.dtype)
+            epsilon.data[:] = model.pad(model.epsilon)
         else:
-            A = Matrix([[1, x1, y1, z1, x1*y1, x1*z1, y1*z1, x1*y1*z1],
-                        [1, x1, y2, z1, x1*y2, x1*z1, y2*z1, x1*y2*z1],
-                        [1, x2, y1, z1, x2*y1, x2*z1, y2*z1, x2*y1*z1],
-                        [1, x1, y1, z2, x1*y1, x1*z2, y1*z2, x1*y1*z2],
-                        [1, x2, y2, z1, x2*y2, x2*z1, y2*z1, x2*y2*z1],
-                        [1, x1, y2, z2, x1*y2, x1*z2, y2*z2, x1*y2*z2],
-                        [1, x2, y1, z2, x2*y1, x2*z2, y1*z2, x2*y1*z2],
-                        [1, x2, y2, z2, x2*y2, x2*z2, y2*z2, x2*y2*z2]])
-            self.increments = (0, 0, 0), (0, 1, 0), (1, 0, 0), (0, 0, 1), (1, 1, 0), (0, 1, 1), (1, 0, 1), (1, 1, 1)
-            self.rs = symbols('rx, ry, rz')
-            rx, ry, rz = self.rs
-            p = Matrix([[1],
-                        [rx],
-                        [ry],
-                        [rz],
-                        [rx*ry],
-                        [rx*rz],
-                        [ry*rz],
-                        [rx*ry*rz]])
-
-        # Map to reference cell
-        reference_cell = [(x1, 0),
-                          (y1, 0),
-                          (z1, 0),
-                          (x2, self.h),
-                          (y2, self.h),
-                          (z2, self.h)]
-
-        A = A.subs(reference_cell)
-        self.bs = A.inv().T.dot(p)
-
-    @property
-    def sym_coordinates(self):
-        """Symbol representing the coordinate values in each dimension"""
-        return tuple([self.coordinates.indexed[p, i]
-                      for i in range(self.ndim)])
-
-    @property
-    def sym_coord_indices(self):
-        """Symbol for each grid index according to the coordinates"""
-        return tuple([Function('INT')(Function('floor')(x / self.h))
-                      for x in self.sym_coordinates])
-
-    @property
-    def sym_coord_bases(self):
-        """Symbol for the base coordinates of the reference grid point"""
-        return tuple([Function('FLOAT')(x - idx * self.h)
-                      for x, idx in zip(self.sym_coordinates,
-                                        self.sym_coord_indices)])
-
-    def point2grid(self, u, m):
-        """Generates an expression for generic point-to-grid interpolation"""
-        dt = self.dt
-        subs = dict(zip(self.rs, self.sym_coord_bases))
-        index_matrix = [tuple([idx + ii + self.nbpml for ii, idx
-                               in zip(inc, self.sym_coord_indices)])
-                        for inc in self.increments]
-        eqns = [Eq(u.indexed[(t, ) + idx], u.indexed[(t, ) + idx]
-                   + self.indexed[t, p] * dt * dt / m.indexed[idx] * b.subs(subs))
-                for idx, b in zip(index_matrix, self.bs)]
-        return eqns
-
-    def grid2point(self, u):
-        """Generates an expression for generic grid-to-point interpolation"""
-        subs = dict(zip(self.rs, self.sym_coord_bases))
-        index_matrix = [tuple([idx + ii + self.nbpml for ii, idx
-                               in zip(inc, self.sym_coord_indices)])
-                        for inc in self.increments]
-        return sum([b.subs(subs) * u.indexed[(t, ) + idx]
-                    for idx, b in zip(index_matrix, self.bs)])
-
-    def read(self, u, v):
-        """Iteration loop over points performing grid-to-point interpolation."""
-        interp_expr = Eq(self.indexed[t, p], self.grid2point(u) + self.grid2point(v))
-        return [Iteration(interp_expr, variable=p, limits=self.shape[1])]
-
-    def add(self, m, u):
-        """Iteration loop over points performing point-to-grid interpolation."""
-        return [Iteration(self.point2grid(u, m), variable=p, limits=self.shape[1])]
+            epsilon = 1.0
 
+        if model.delta is not None:
+            delta = DenseData(name="delta", shape=model.get_shape_comp(),
+                              dtype=damp.dtype)
+            delta.data[:] = model.pad(model.delta)
+        else:
+            delta = 1.0
+        if model.theta is not None:
+            theta = DenseData(name="theta", shape=model.get_shape_comp(),
+                              dtype=damp.dtype)
+            theta.data[:] = model.pad(model.theta)
+        else:
+            theta = 0
+
+        if len(model.get_shape_comp()) == 3:
+            if model.phi is not None:
+                phi = DenseData(name="phi", shape=model.get_shape_comp(),
+                                dtype=damp.dtype)
+                phi.data[:] = model.pad(model.phi)
+            else:
+                phi = 0
+
+        u.pad_time = save
+        v.pad_time = save
+        rec = SourceLike(name="rec", npoint=nrec, nt=nt, dt=dt,
+                         h=model.get_spacing(),
+                         coordinates=data.receiver_coords,
+                         ndim=len(damp.shape),
+                         dtype=damp.dtype,
+                         nbpml=model.nbpml)
 
-class ForwardOperator(Operator):
-    def __init__(self, m, src, damp, rec, u, v, A, B, th, ph, time_order=2, spc_order=4, **kwargs):
         def Bhaskarasin(angle):
             if angle == 0:
                 return 0
@@ -135,23 +81,15 @@ def Bhaskaracos(angle):
         else:
             ang0 = Function('ang0')(x, y)
             ang1 = Function('ang1')(x, y)
-        assert(m.shape == damp.shape)
-        u.pad_time = False
-        v.pad_time = False
-
-        # Set time and space orders
-        u.time_order = time_order
-        u.space_order = spc_order
-        v.time_order = time_order
-        v.space_order = spc_order
+
         s, h = symbols('s h')
 
-        ang0 = Bhaskaracos(th)
-        ang1 = Bhaskarasin(th)
+        ang0 = Bhaskaracos(theta)
+        ang1 = Bhaskarasin(theta)
         # Derive stencil from symbolic equation
         if len(m.shape) == 3:
-            ang2 = Bhaskaracos(ph)
-            ang3 = Bhaskarasin(ph)
+            ang2 = Bhaskaracos(phi)
+            ang3 = Bhaskarasin(phi)
 
             Gy1p = (ang3 * u.dxl - ang2 * u.dyl)
             Gyy1 = (first_derivative(Gy1p, ang3, dim=x, side=1, order=spc_order/2) -
@@ -179,11 +117,11 @@ def Bhaskaracos(angle):
             Gzz2 = (first_derivative(Gz2r, ang1, ang2, dim=x, side=-1, order=spc_order/2) +
                     first_derivative(Gz2r, ang1, ang3, dim=y, side=-1, order=spc_order/2) +
                     first_derivative(Gz2r, ang0, dim=z, side=-1, order=spc_order/2))
-            parm = [m, damp, A, B, th, ph, u, v]
+            parm = [m, damp, epsilon, delta, theta, phi, u, v]
         else:
             Gyy2 = 0
             Gyy1 = 0
-            parm = [m, damp, A, B, th, u, v]
+            parm = [m, damp, epsilon, delta, theta, u, v]
             Gx1p = (ang0 * u.dxl - ang1 * u.dyl)
             Gz1r = (ang1 * v.dxl + ang0 * v.dyl)
             Gxx1 = (first_derivative(Gx1p * ang0, dim=x, side=1, order=spc_order/2) -
@@ -197,10 +135,12 @@ def Bhaskaracos(angle):
             Gzz2 = (first_derivative(Gz2r * ang1, dim=x, side=-1, order=spc_order/2) +
                     first_derivative(Gz2r * ang0, dim=y, side=-1, order=spc_order/2))
 
-        stencilp = 1.0 / (2.0 * m + s * damp) * (4.0 * m * u + (s * damp - 2.0 * m) * u.backward +
-                                                 2.0 * s**2 * (A * Hp + B * Hzr))
-        stencilr = 1.0 / (2.0 * m + s * damp) * (4.0 * m * v + (s * damp - 2.0 * m) * v.backward +
-                                                 2.0 * s**2 * (B * Hp + Hzr))
+        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))
+        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))
         Hp = -(.5 * Gxx1 + .5 * Gxx2 + .5 * Gyy1 + .5 * Gyy2)
         Hzr = -(.5 * Gzz1 + .5 * Gzz2)
         factorized = {"Hp": Hp, "Hzr": Hzr}
@@ -209,14 +149,21 @@ def Bhaskaracos(angle):
         first_stencil = Eq(u.forward, stencilp)
         second_stencil = Eq(v.forward, stencilr)
         stencils = [first_stencil, second_stencil]
-        super(ForwardOperator, self).__init__(src.nt, m.shape, stencils=stencils, substitutions=subs,
-                                              spc_border=spc_order/2, time_order=time_order, forward=True,
+        super(ForwardOperator, self).__init__(src.nt, m.shape,
+                                              stencils=stencils,
+                                              substitutions=subs,
+                                              spc_border=spc_order/2,
+                                              time_order=time_order,
+                                              forward=True,
                                               dtype=m.dtype,
-                                              input_params=parm, factorized=factorized, **kwargs)
+                                              input_params=parm,
+                                              factorized=factorized,
+                                              **kwargs)
 
         # Insert source and receiver terms post-hoc
         self.input_params += [src, src.coordinates, rec, rec.coordinates]
-        self.propagator.time_loop_stencils_a = src.add(m, u) + src.add(m, v) + rec.read(u, v)
+        self.output_params += [v, rec]
+        self.propagator.time_loop_stencils_a = src.add(m, u) + src.add(m, v) + rec.read2(u, v)
         self.propagator.add_devito_param(src)
         self.propagator.add_devito_param(src.coordinates)
         self.propagator.add_devito_param(rec)
@@ -228,23 +175,30 @@ def __init__(self, m, rec, damp, srca, time_order=4, spc_order=12):
         assert(m.shape == damp.shape)
 
         input_params = [m, rec, damp, srca]
-        v = TimeData("v", m.shape, rec.nt, time_order=time_order, save=True, dtype=m.dtype)
+        v = TimeData("v", m.shape, rec.nt, time_order=time_order,
+                     save=True, dtype=m.dtype)
         output_params = [v]
         dim = len(m.shape)
         total_dim = self.total_dim(dim)
         space_dim = self.space_dim(dim)
         lhs = v.indexed[total_dim]
         stencil, subs = self._init_taylor(dim, time_order, spc_order)[1]
-        stencil = self.smart_sympy_replace(dim, time_order, stencil, Function('p'), v, fw=False)
+        stencil = self.smart_sympy_replace(dim, time_order, stencil,
+                                           Function('p'), v, fw=False)
         main_stencil = Eq(lhs, stencil)
-        stencil_args = [m.indexed[space_dim], rec.dt, rec.h, damp.indexed[space_dim]]
+        stencil_args = [m.indexed[space_dim], rec.dt, rec.h,
+                        damp.indexed[space_dim]]
         stencils = [main_stencil]
         substitutions = [dict(zip(subs, stencil_args))]
 
-        super(AdjointOperator, self).__init__(rec.nt, m.shape, stencils=stencils,
-                                              substitutions=substitutions, spc_border=spc_order/2,
-                                              time_order=time_order, forward=False, dtype=m.dtype,
-                                              input_params=input_params, output_params=output_params)
+        super(AdjointOperator, self).__init__(rec.nt, m.shape,
+                                              stencils=stencils,
+                                              substitutions=substitutions,
+                                              spc_border=spc_order/2,
+                                              time_order=time_order,
+                                              forward=False, dtype=m.dtype,
+                                              input_params=input_params,
+                                              output_params=output_params)
 
         # Insert source and receiver terms post-hoc
         self.propagator.time_loop_stencils_a = rec.add(m, v) + srca.read(v)
@@ -255,7 +209,8 @@ def __init__(self, u, m, rec, damp, time_order=4, spc_order=12):
         assert(m.shape == damp.shape)
 
         input_params = [u, m, rec, damp]
-        v = TimeData("v", m.shape, rec.nt, time_order=time_order, save=False, dtype=m.dtype)
+        v = TimeData("v", m.shape, rec.nt, time_order=time_order,
+                     save=False, dtype=m.dtype)
         grad = DenseData("grad", m.shape, dtype=m.dtype)
         output_params = [grad, v]
         dim = len(m.shape)
@@ -263,8 +218,10 @@ def __init__(self, u, m, rec, damp, time_order=4, spc_order=12):
         space_dim = self.space_dim(dim)
         lhs = v.indexed[total_dim]
         stencil, subs = self._init_taylor(dim, time_order, spc_order)[1]
-        stencil = self.smart_sympy_replace(dim, time_order, stencil, Function('p'), v, fw=False)
-        stencil_args = [m.indexed[space_dim], rec.dt, rec.h, damp.indexed[space_dim]]
+        stencil = self.smart_sympy_replace(dim, time_order, stencil,
+                                           Function('p'), v, fw=False)
+        stencil_args = [m.indexed[space_dim], rec.dt, rec.h,
+                        damp.indexed[space_dim]]
         main_stencil = Eq(lhs, lhs + stencil)
         gradient_update = Eq(grad.indexed[space_dim],
                              grad.indexed[space_dim] -
@@ -274,10 +231,14 @@ def __init__(self, u, m, rec, damp, time_order=4, spc_order=12):
         stencils = [main_stencil, gradient_update, reset_v]
         substitutions = [dict(zip(subs, stencil_args)), {}, {}]
 
-        super(GradientOperator, self).__init__(rec.nt, m.shape, stencils=stencils,
-                                               substitutions=substitutions, spc_border=spc_order/2,
-                                               time_order=time_order, forward=False, dtype=m.dtype,
-                                               input_params=input_params, output_params=output_params)
+        super(GradientOperator, self).__init__(rec.nt, m.shape,
+                                               stencils=stencils,
+                                               substitutions=substitutions,
+                                               spc_border=spc_order/2,
+                                               time_order=time_order,
+                                               forward=False, dtype=m.dtype,
+                                               input_params=input_params,
+                                               output_params=output_params)
 
         # Insert source and receiver terms post-hoc
         self.propagator.time_loop_stencils_b = rec.add(m, v)
@@ -288,8 +249,10 @@ def __init__(self, dm, m, src, damp, rec, time_order=4, spc_order=12):
         assert(m.shape == damp.shape)
 
         input_params = [dm, m, src, damp, rec]
-        u = TimeData("u", m.shape, src.nt, time_order=time_order, save=False, dtype=m.dtype)
-        U = TimeData("U", m.shape, src.nt, time_order=time_order, save=False, dtype=m.dtype)
+        u = TimeData("u", m.shape, src.nt, time_order=time_order,
+                     save=False, dtype=m.dtype)
+        U = TimeData("U", m.shape, src.nt, time_order=time_order,
+                     save=False, dtype=m.dtype)
         output_params = [u, U]
         dim = len(m.shape)
         total_dim = self.total_dim(dim)
@@ -305,16 +268,22 @@ def __init__(self, dm, m, src, damp, rec, time_order=4, spc_order=12):
                 u.indexed[tuple((t - 2,) + space_dim)])*dm.indexed[space_dim])
         second_stencil_args = [m.indexed[space_dim], dt, h, damp.indexed[space_dim]]
         second_update = Eq(U.indexed[total_dim], second_stencil)
-        insert_second_source = Eq(U.indexed[total_dim], U.indexed[total_dim]+(dt*dt)/m.indexed[space_dim]*src2)
+        insert_second_source = Eq(U.indexed[total_dim], U.indexed[total_dim] +
+                                  (dt*dt)/m.indexed[space_dim]*src2)
         reset_u = Eq(u.indexed[tuple((t - 2,) + space_dim)], 0)
         stencils = [first_update, second_update, insert_second_source, reset_u]
         substitutions = [dict(zip(subs, first_stencil_args)),
                          dict(zip(subs, second_stencil_args)), {}, {}]
 
-        super(BornOperator, self).__init__(src.nt, m.shape, stencils=stencils,
-                                           substitutions=substitutions, spc_border=spc_order/2,
-                                           time_order=time_order, forward=True, dtype=m.dtype,
-                                           input_params=input_params, output_params=output_params)
+        super(BornOperator, self).__init__(src.nt, m.shape,
+                                           stencils=stencils,
+                                           substitutions=substitutions,
+                                           spc_border=spc_order/2,
+                                           time_order=time_order,
+                                           forward=True,
+                                           dtype=m.dtype,
+                                           input_params=input_params,
+                                           output_params=output_params)
 
         # Insert source and receiver terms post-hoc
         self.propagator.time_loop_stencils_b = src.add(m, u)
diff --git a/tests/test_adjointA.py b/tests/test_adjointA.py
old mode 100644
new mode 100755
index 4cc594b24c..450f1487ca
--- a/tests/test_adjointA.py
+++ b/tests/test_adjointA.py
@@ -6,21 +6,21 @@
 from examples.containers import IGrid, IShot
 
 
-class Test_AdjointA(object):
-    @pytest.fixture(params=[(60, 70), (50, 60, 70)])
-    def Acoustic(self, request, time_order, space_order):
+class TestAdjointA(object):
+    @pytest.fixture(params=[(60, 70), (60, 70, 80)])
+    def acoustic(self, request, time_order, space_order):
         model = IGrid()
         dimensions = request.param
         # dimensions are (x,z) and (x, y, z)
-        origin = tuple([0]*len(dimensions))
-        spacing = tuple([25]*len(dimensions))
+        origin = tuple([0.0]*len(dimensions))
+        spacing = tuple([15.0]*len(dimensions))
 
         # True velocity
-        true_vp = np.ones(dimensions) + 2.0
+        true_vp = np.ones(dimensions) + .5
         if len(dimensions) == 2:
-            true_vp[:, int(dimensions[0] / 2):dimensions[0]] = 4.5
+            true_vp[:, int(dimensions[0] / 2):dimensions[0]] = 2.5
         else:
-            true_vp[:, :, int(dimensions[0] / 2):dimensions[0]] = 4.5
+            true_vp[:, :, int(dimensions[0] / 2):dimensions[0]] = 2.5
         model.create_model(origin, spacing, true_vp)
         # Define seismic data.
         data = IShot()
@@ -28,7 +28,7 @@ def Acoustic(self, request, time_order, space_order):
         f0 = .010
         dt = model.get_critical_dt()
         t0 = 0.0
-        tn = 100.0
+        tn = 500.0
         nt = int(1+(tn-t0)/dt)
 
         # Set up the source as Ricker wavelet for f0
@@ -40,47 +40,48 @@ def source(t, f0):
         location = (origin[0] + dimensions[0] * spacing[0] * 0.5,
                     origin[-1] + 2 * spacing[-1])
         if len(dimensions) == 3:
-            location = (location[0], 0., location[1])
+            location = (location[0], origin[1] + dimensions[1] * spacing[1] * 0.5, location[1])
         data.set_source(time_series, dt, location)
-        receiver_coords = np.zeros((30, len(dimensions)))
-        receiver_coords[:, 0] = np.linspace(50, origin[0] + dimensions[0]*spacing[0] - 50, num=30)
+        receiver_coords = np.zeros((50, len(dimensions)))
+        receiver_coords[:, 0] = np.linspace(50, origin[0] + dimensions[0]*spacing[0] - 50, num=50)
         receiver_coords[:, -1] = location[-1]
-
+        if len(dimensions) == 3:
+            receiver_coords[:, -1] = location[1]
         data.set_receiver_pos(receiver_coords)
-        data.set_shape(nt, 30)
+        data.set_shape(nt, 50)
         # Adjoint test
-        wave_true = Acoustic_cg(model, data, None, None, t_order=time_order, s_order=space_order, nbpml=10)
+        wave_true = Acoustic_cg(model, data, t_order=time_order, s_order=space_order, nbpml=10)
         return wave_true
 
     @pytest.fixture(params=[2])
     def time_order(self, request):
         return request.param
 
-    @pytest.fixture(params=[4, 6, 8, 10])
+    @pytest.fixture(params=[2, 4, 6, 8, 10, 12])
     def space_order(self, request):
         return request.param
 
     @pytest.fixture
-    def forward(self, Acoustic):
-        (rec, u) = Acoustic.Forward()
-        return (rec, u)
+    def forward(self, acoustic):
+        rec, u = acoustic.Forward(save=False)
+        return rec
 
-    def test_adjoint(self, Acoustic, forward):
-        rec = forward[0]
-        srca, v = Acoustic.Adjoint(rec)
-        nt = Acoustic.nt
+    def test_adjoint(self, acoustic, forward):
+        rec = forward
+        srca = acoustic.Adjoint(rec)
+        nt = srca.shape[0]
         # Actual adjoint test
         term1 = 0
         for ti in range(0, nt):
-            term1 = term1 + srca[ti] * Acoustic.data.get_source(ti)
+            term1 = term1 + srca[ti] * acoustic.data.get_source(ti)
         term2 = linalg.norm(rec)**2
         print(term1, term2, term1 - term2, term1 / term2)
-        assert np.isclose(term1 / term2, 1.0)
+        assert np.isclose(term1 / term2, 1.0, atol=0.001)
 
 if __name__ == "__main__":
-    t = Test_AdjointA()
+    t = TestAdjointA()
     request = type('', (), {})()
     request.param = (60, 70, 80)
-    ac = t.Acoustic(request, 2, 12)
+    ac = t.acoustic(request, 2, 12)
     fw = t.forward(ac)
     t.test_adjoint(ac, fw)
diff --git a/tests/test_gradient.py b/tests/test_gradient.py
new file mode 100755
index 0000000000..8d4b0a9626
--- /dev/null
+++ b/tests/test_gradient.py
@@ -0,0 +1,117 @@
+import numpy as np
+import pytest
+from numpy import linalg
+
+from examples.Acoustic_codegen import Acoustic_cg
+from examples.containers import IGrid, IShot
+
+
+class TestGradient(object):
+    @pytest.fixture(params=[(70, 80), (60, 70, 80)])
+    def acoustic(self, request, time_order, space_order):
+        model = IGrid()
+        model0 = IGrid()
+        dimensions = request.param
+        # dimensions are (x,z) and (x, y, z)
+        origin = tuple([0]*len(dimensions))
+        spacing = tuple([15]*len(dimensions))
+
+        # velocity models
+        def smooth10(vel):
+            out = np.zeros(dimensions)
+            out[:] = vel[:]
+            for a in range(5, dimensions[-1]-6):
+                if len(dimensions) == 2:
+                    out[:, a] = np.sum(vel[:, a - 5:a + 5], axis=1) / 10
+                else:
+                    out[:, :, a] = np.sum(vel[:, :, a - 5:a + 5], axis=2) / 10
+            return out
+
+        # True velocity
+        true_vp = np.ones(dimensions) + .5
+        if len(dimensions) == 2:
+            true_vp[:, int(dimensions[1] / 3):] = 2.5
+        else:
+            true_vp[:, :, int(dimensions[2] / 3):] = 2.5
+        # Smooth velocity
+        initial_vp = smooth10(true_vp)
+        dm = true_vp**-2 - initial_vp**-2
+        model.create_model(origin, spacing, true_vp)
+        model0.create_model(origin, spacing, initial_vp)
+        # Define seismic data.
+        data = IShot()
+        f0 = .010
+        dt = model.get_critical_dt()
+        t0 = 0.0
+        tn = 500.0
+        nt = int(1+(tn-t0)/dt)
+        # Set up the source as Ricker wavelet for f0
+
+        def source(t, f0):
+            r = (np.pi * f0 * (t - 1./f0))
+            return (1-2.*r**2)*np.exp(-r**2)
+
+        time_series = source(np.linspace(t0, tn, nt), f0)
+        location = (origin[0] + dimensions[0] * spacing[0] * 0.5,
+                    origin[-1] + 2 * spacing[-1])
+        if len(dimensions) == 3:
+            location = (location[0], origin[1] + dimensions[1] * spacing[1] * 0.5, location[1])
+        data.set_source(time_series, dt, location)
+        receiver_coords = np.zeros((50, len(dimensions)))
+        receiver_coords[:, 0] = np.linspace(50, origin[0] + dimensions[0]*spacing[0] - 50, num=50)
+        receiver_coords[:, -1] = location[-1]
+        if len(dimensions) == 3:
+            receiver_coords[:, 1] = location[1]
+        data.set_receiver_pos(receiver_coords)
+        data.set_shape(nt, 50)
+        # Adjoint test
+        wave_true = Acoustic_cg(model, data, None, t_order=time_order, s_order=space_order, nbpml=10)
+        wave_0 = Acoustic_cg(model0, data, None, t_order=time_order, s_order=space_order, nbpml=10)
+        return wave_true, wave_0, dm, initial_vp
+
+    @pytest.fixture(params=[2])
+    def time_order(self, request):
+        return request.param
+
+    @pytest.fixture(params=[2])
+    def space_order(self, request):
+        return request.param
+
+    def test_grad(self, acoustic):
+        rec = acoustic[0].Forward()[0]
+        rec0, u0 = acoustic[1].Forward(save=True)
+        F0 = .5*linalg.norm(rec0 - rec)**2
+        gradient = acoustic[1].Gradient(rec0 - rec, u0)
+        # Actual Gradient test
+        G = np.dot(gradient.reshape(-1), acoustic[1].model.pad(acoustic[2]).reshape(-1))
+        # FWI Gradient test
+        H = [1, 0.1, 0.01, .001, 0.0001, 0.00001, 0.000001]
+        error1 = np.zeros(7)
+        error2 = np.zeros(7)
+        for i in range(0, 7):
+            acoustic[1].model.set_vp(np.sqrt((acoustic[3]**-2 + H[i] * acoustic[2])**(-1)))
+            d = acoustic[1].Forward()[0]
+            error1[i] = np.absolute(.5*linalg.norm(d - rec)**2 - F0)
+            error2[i] = np.absolute(.5*linalg.norm(d - rec)**2 - F0 - H[i] * G)
+            # print(F0, .5*linalg.norm(d - rec)**2, error1[i], H[i] *G, error2[i])
+            # print('For h = ', H[i], '\nFirst order errors is : ', error1[i],
+            #       '\nSecond order errors is ', error2[i])
+
+        hh = np.zeros(7)
+        for i in range(0, 7):
+            hh[i] = H[i] * H[i]
+
+        # Test slope of the  tests
+        p1 = np.polyfit(np.log10(H), np.log10(error1), 1)
+        p2 = np.polyfit(np.log10(H), np.log10(error2), 1)
+        print(p1)
+        print(p2)
+        assert np.isclose(p1[0], 1.0, rtol=0.05)
+        assert np.isclose(p2[0], 2.0, rtol=0.05)
+
+if __name__ == "__main__":
+    t = TestGradient()
+    request = type('', (), {})()
+    request.param = (60, 70, 80)
+    ac = t.acoustic(request, 2, 12)
+    t.test_grad(ac)