diff --git a/src/TimeModeling/Modeling/misfit_fg.jl b/src/TimeModeling/Modeling/misfit_fg.jl
index 4ad179cb7..585e21b01 100644
--- a/src/TimeModeling/Modeling/misfit_fg.jl
+++ b/src/TimeModeling/Modeling/misfit_fg.jl
@@ -48,7 +48,7 @@ function multi_src_fg(model_full::AbstractModel, source::judiVector, dObs::judiV
     @juditime "Python call to J_adjoint" begin
         argout = rlock_pycall(ac."J_adjoint", Tuple{Float32, PyArray, IT...}, modelPy,
                 src_coords, qIn, rec_coords, dObserved, t_sub=options.subsampling_factor,
-                space_order=options.space_order, checkpointing=options.optimal_checkpointing,
+                checkpointing=options.optimal_checkpointing,
                 freq_list=freqs, ic=options.IC, is_residual=false, born_fwd=lin, nlind=nlind,
                 dft_sub=options.dft_subsampling_factor[1], f0=options.f0, return_obj=true,
                 misfit=mfunc, illum=illum)
diff --git a/src/TimeModeling/Modeling/python_interface.jl b/src/TimeModeling/Modeling/python_interface.jl
index 50c42d636..7e49d9d3a 100644
--- a/src/TimeModeling/Modeling/python_interface.jl
+++ b/src/TimeModeling/Modeling/python_interface.jl
@@ -62,7 +62,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Geometry, srcData::Arr
     rec_coords = setup_grid(recGeometry, modelPy.shape)
 
     # Devito call
-    return wrapcall_data(ac."forward_rec", modelPy, src_coords, qIn, rec_coords, fw=fw, space_order=options.space_order, f0=options.f0, illum=illum)
+    return wrapcall_data(ac."forward_rec", modelPy, src_coords, qIn, rec_coords, fw=fw, f0=options.f0, illum=illum)
 end
 
 # u = F*Ps'*q
@@ -76,7 +76,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Geometry, srcData::Arr
     src_coords = setup_grid(srcGeometry, modelPy.shape)
 
     # Devito call
-    return wrapcall_wf(ac."forward_no_rec", modelPy, src_coords, qIn, fw=fw, space_order=options.space_order, illum=illum)
+    return wrapcall_wf(ac."forward_no_rec", modelPy, src_coords, qIn, fw=fw, illum=illum)
 end
 
 # d_obs = Pr*F*u
@@ -89,7 +89,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Nothing, srcData::Arra
     rec_coords = setup_grid(recGeometry, modelPy.shape)
 
     # Devito call
-    return wrapcall_data(ac."forward_wf_src", modelPy, srcData, rec_coords, fw=fw, space_order=options.space_order, f0=options.f0, illum=illum)
+    return wrapcall_data(ac."forward_wf_src", modelPy, srcData, rec_coords, fw=fw, f0=options.f0, illum=illum)
 end
 
 # u_out = F*u_in
@@ -99,7 +99,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Nothing, srcData::Arra
     dtComp = convert(Float32, modelPy."critical_dt")
 
     # Devito call
-    return wrapcall_wf(ac."forward_wf_src_norec", modelPy, srcData, fw=fw, space_order=options.space_order, illum=illum)
+    return wrapcall_wf(ac."forward_wf_src_norec", modelPy, srcData, fw=fw, illum=illum)
 end
 
 # d_lin = J*dm
@@ -117,7 +117,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Geometry, srcData::Arr
 
     # Devito call
     return wrapcall_data(ac."born_rec", modelPy, src_coords, qIn, rec_coords, fw=fw,
-                  space_order=options.space_order, ic=options.IC, f0=options.f0, illum=illum)
+                  ic=options.IC, f0=options.f0, illum=illum)
 end
 
 # dm = J'*d_lin
@@ -136,7 +136,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Geometry, srcData::Arr
     length(options.frequencies) == 0 ? freqs = nothing : freqs = options.frequencies
     return wrapcall_grad(ac."J_adjoint", modelPy,
                   src_coords, qIn, rec_coords, dIn, fw=fw, t_sub=options.subsampling_factor,
-                  space_order=options.space_order, checkpointing=options.optimal_checkpointing,
+                  checkpointing=options.optimal_checkpointing,
                   freq_list=freqs, ic=options.IC, is_residual=true,
                   dft_sub=options.dft_subsampling_factor[1], f0=options.f0, illum=illum)
 end
@@ -157,7 +157,7 @@ function devito_interface(modelPy::PyObject, weights::Array, srcData::Array, rec
 
     # Devito call
     return wrapcall_data(ac."forward_rec_w", modelPy, weights, qIn, rec_coords,
-                         fw=fw, space_order=options.space_order, f0=options.f0, illum=illum)
+                         fw=fw, f0=options.f0, illum=illum)
 end
 
 # dw = Pw*F'*Pr'*d_obs - adjoint modeling w/ extended source
@@ -173,7 +173,7 @@ function devito_interface(modelPy::PyObject, recGeometry::Geometry, recData::Arr
 
     # Devito call
     return wrapcall_weights(ac."adjoint_w", modelPy, rec_coords, dIn, qIn,
-                            fw=fw, space_order=options.space_order, f0=options.f0, illum=illum)
+                            fw=fw, f0=options.f0, illum=illum)
 end
 
 # Jacobian of extended source modeling: d_lin = J*dm
@@ -190,7 +190,7 @@ function devito_interface(modelPy::PyObject, weights::Array, srcData::Array, rec
 
     # Devito call
     return wrapcall_data(ac."born_rec_w", modelPy, weights, qIn, rec_coords,
-                         fw=fw, space_order=options.space_order, ic=options.IC, f0=options.f0, illum=illum)
+                         fw=fw, ic=options.IC, f0=options.f0, illum=illum)
 end
 
 # Adjoint Jacobian of extended source modeling: dm = J'*d_lin
@@ -208,7 +208,7 @@ function devito_interface(modelPy::PyObject, weights::Array, srcData::Array, rec
     length(options.frequencies) == 0 ? freqs = nothing : freqs = options.frequencies
     return wrapcall_grad(ac."J_adjoint", modelPy,
                   nothing, qIn, rec_coords, dIn, fw=fw, t_sub=options.subsampling_factor,
-                  space_order=options.space_order, checkpointing=options.optimal_checkpointing,
+                  checkpointing=options.optimal_checkpointing,
                   freq_list=freqs, ic=options.IC, ws=weights, is_residual=true,
                   dft_sub=options.dft_subsampling_factor[1], f0=options.f0, illum=illum)
 end
diff --git a/src/pysource/FD_utils.py b/src/pysource/FD_utils.py
index 81738952b..2dd7977cf 100644
--- a/src/pysource/FD_utils.py
+++ b/src/pysource/FD_utils.py
@@ -9,8 +9,7 @@ def laplacian(v, irho):
     """
     irho = irho or 1
     if isinstance(irho, Differentiable):
-        so = irho.space_order // 2
-        Lap = div(irho * grad(v, shift=.5, order=so), shift=-.5, order=so)
+        Lap = div(irho * grad(v, shift=.5), shift=-.5)
     else:
         Lap = irho * v.laplace
 
@@ -93,16 +92,12 @@ def sa_tti(u, v, model):
     model: Model
         Model structure
     """
-    # Space order
-    so = u.space_order // 2
     # Matrix of Thomsen params
     A, B, C = thomsen_mat(model)
     # Rotation Matrix
     R = R_mat(model)
 
-    PI = R.T * (A * R * grad(u, shift=.5, order=so) +
-                B * R * grad(v, shift=.5, order=so))
-    MI = R.T * (B * R * grad(u, shift=.5, order=so) +
-                C * R * grad(v, shift=.5, order=so))
+    PI = R.T * (A * R * grad(u, shift=.5) + B * R * grad(v, shift=.5))
+    MI = R.T * (B * R * grad(u, shift=.5) + C * R * grad(v, shift=.5))
 
-    return div(PI, shift=-.5, order=so), div(MI, shift=-.5, order=so)
+    return div(PI, shift=-.5), div(MI, shift=-.5)
diff --git a/src/pysource/interface.py b/src/pysource/interface.py
index 319a8b0eb..0b00945f4 100644
--- a/src/pysource/interface.py
+++ b/src/pysource/interface.py
@@ -14,7 +14,7 @@
 
 
 # Forward wrappers Pr*F*Ps'*q
-def forward_rec(model, src_coords, wavelet, rec_coords, space_order=8, f0=0.015,
+def forward_rec(model, src_coords, wavelet, rec_coords, f0=0.015,
                 illum=False, fw=True):
     """
     Modeling of a point source with receivers Pr*F*Ps^T*q.
@@ -29,8 +29,6 @@ def forward_rec(model, src_coords, wavelet, rec_coords, space_order=8, f0=0.015,
         Source signature
     rec_coords: Array
         Coordiantes of the receiver(s)
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     f0: float
         peak frequency
     illum: bool
@@ -44,12 +42,12 @@ def forward_rec(model, src_coords, wavelet, rec_coords, space_order=8, f0=0.015,
         Shot record
     """
     rec, _, I, _ = forward(model, src_coords, rec_coords, wavelet, save=False,
-                           space_order=space_order, f0=f0, illum=illum, fw=fw)
+                           f0=f0, illum=illum, fw=fw)
     return rec.data, getattr(I, "data", None)
 
 
 #  Pr*F*Pw'*w
-def forward_rec_w(model, weight, wavelet, rec_coords, space_order=8, f0=0.015,
+def forward_rec_w(model, weight, wavelet, rec_coords, f0=0.015,
                   illum=False, fw=True):
     """
     Forward modeling of an extended source with receivers  Pr*F*Pw^T*w
@@ -64,8 +62,6 @@ def forward_rec_w(model, weight, wavelet, rec_coords, space_order=8, f0=0.015,
         Source signature
     rec_coords: Array
         Coordiantes of the receiver(s)
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     f0: float
         peak frequency
     illum: bool
@@ -79,12 +75,12 @@ def forward_rec_w(model, weight, wavelet, rec_coords, space_order=8, f0=0.015,
         Shot record
     """
     rec, _, I, _ = forward(model, None, rec_coords, wavelet, save=False, ws=weight,
-                           space_order=space_order, f0=f0, illum=illum, fw=fw)
+                           f0=f0, illum=illum, fw=fw)
     return rec.data, getattr(I, "data", None)
 
 
 # F*Ps'*q
-def forward_no_rec(model, src_coords, wavelet, space_order=8, f0=0.015, illum=False,
+def forward_no_rec(model, src_coords, wavelet, f0=0.015, illum=False,
                    fw=True):
     """
     Forward modeling of a point source without receiver.
@@ -97,8 +93,6 @@ def forward_no_rec(model, src_coords, wavelet, space_order=8, f0=0.015, illum=Fa
         Coordiantes of the source(s)
     wavelet: Array
         Source signature
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     f0: float
         peak frequency
     illum: bool
@@ -111,13 +105,13 @@ def forward_no_rec(model, src_coords, wavelet, space_order=8, f0=0.015, illum=Fa
     Array
         Wavefield
     """
-    _, u, I, _ = forward(model, src_coords, None, wavelet, space_order=space_order,
+    _, u, I, _ = forward(model, src_coords, None, wavelet,
                          save=True, f0=f0, illum=illum, fw=fw)
     return u.data, getattr(I, "data", None)
 
 
 # Pr*F*u
-def forward_wf_src(model, u, rec_coords, space_order=8, f0=0.015, illum=False, fw=True):
+def forward_wf_src(model, u, rec_coords, f0=0.015, illum=False, fw=True):
     """
     Forward modeling of a full wavefield source Pr*F*u.
 
@@ -129,8 +123,6 @@ def forward_wf_src(model, u, rec_coords, space_order=8, f0=0.015, illum=False, f
         Time-space dependent wavefield
     rec_coords: Array
         Coordiantes of the receiver(s)
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     f0: float
         peak frequency
     illum: bool
@@ -144,13 +136,13 @@ def forward_wf_src(model, u, rec_coords, space_order=8, f0=0.015, illum=False, f
         Shot record
     """
     wsrc = src_wavefield(model, u, fw=True)
-    rec, _, I, _ = forward(model, None, rec_coords, None, space_order=space_order,
+    rec, _, I, _ = forward(model, None, rec_coords, None,
                            qwf=wsrc, illum=illum, f0=f0, fw=fw)
     return rec.data, getattr(I, "data", None)
 
 
 # F*u
-def forward_wf_src_norec(model, u, space_order=8, f0=0.015, illum=False, fw=True):
+def forward_wf_src_norec(model, u, f0=0.015, illum=False, fw=True):
     """
     Forward modeling of a full wavefield source without receiver F*u.
 
@@ -160,8 +152,6 @@ def forward_wf_src_norec(model, u, space_order=8, f0=0.015, illum=False, fw=True
         Physical model
     u: TimeFunction or Array
         Time-space dependent wavefield
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     f0: float
         peak frequency
     illum: bool
@@ -175,13 +165,13 @@ def forward_wf_src_norec(model, u, space_order=8, f0=0.015, illum=False, fw=True
         Wavefield
     """
     wf_src = src_wavefield(model, u, fw=True)
-    _, u, I, _ = forward(model, None, None, None, space_order=space_order, save=True,
+    _, u, I, _ = forward(model, None, None, None, save=True,
                          qwf=wf_src, f0=f0, illum=illum, fw=fw)
     return u.data, getattr(I, "data", None)
 
 
 # Pw*F'*Pr'*d_obs
-def adjoint_w(model, rec_coords, data, wavelet, space_order=8, f0=0.015, illum=False,
+def adjoint_w(model, rec_coords, data, wavelet, f0=0.015, illum=False,
               fw=True):
     """
     Adjoint/backward modeling of a shot record (receivers as source) for an
@@ -197,8 +187,6 @@ def adjoint_w(model, rec_coords, data, wavelet, space_order=8, f0=0.015, illum=F
         Shot gather
     wavelet: Array
         Time signature of the forward source for stacking along time
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     f0: float
         peak frequency
     illum: bool
@@ -212,13 +200,13 @@ def adjoint_w(model, rec_coords, data, wavelet, space_order=8, f0=0.015, illum=F
         spatial distribution
     """
     w, _, I, _ = forward(model, rec_coords, None, data, wr=wavelet,
-                         space_order=space_order, f0=f0, illum=illum, fw=fw)
+                         f0=f0, illum=illum, fw=fw)
     return w.data, getattr(I, "data", None)
 
 
 # Linearized modeling ∂/∂m (Pr*F*Ps'*q)
 def born_rec(model, src_coords, wavelet, rec_coords,
-             space_order=8, ic="as", f0=0.015, illum=False, fw=True):
+             ic="as", f0=0.015, illum=False, fw=True):
     """
     Linearized (Born) modeling of a point source for a model perturbation
     (square slowness) dm.
@@ -233,8 +221,6 @@ def born_rec(model, src_coords, wavelet, rec_coords,
         Source signature
     rec_coords: Array
         Coordiantes of the receiver(s)
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     ic: String
         Imaging conditions ("as", "isic" or "fwi"), defaults to "as"
     f0: float
@@ -250,13 +236,13 @@ def born_rec(model, src_coords, wavelet, rec_coords,
         Shot record
     """
     rec, _, I, _ = born(model, src_coords, rec_coords, wavelet, save=False,
-                        space_order=space_order, ic=ic, f0=f0, illum=illum, fw=fw)
+                        ic=ic, f0=f0, illum=illum, fw=fw)
     return rec.data, getattr(I, "data", None)
 
 
 # ∂/∂m (Pr*F*Pw'*w)
 def born_rec_w(model, weight, wavelet, rec_coords,
-               space_order=8, ic="as", f0=0.015, illum=False, fw=True):
+               ic="as", f0=0.015, illum=False, fw=True):
     """
     Linearized (Born) modeling of an extended source for a model
     perturbation (square slowness) dm with an extended source
@@ -271,8 +257,6 @@ def born_rec_w(model, weight, wavelet, rec_coords,
         Source signature
     rec_coords: Array
         Coordiantes of the receiver(s)
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     ic: String
         Imaging conditions ("as", "isic" or "fwi"), defaults to "as"
     f0: float
@@ -288,11 +272,11 @@ def born_rec_w(model, weight, wavelet, rec_coords,
         Shot record
     """
     rec, _, I, _ = born(model, None, rec_coords, wavelet, save=False, ws=weight,
-                        space_order=space_order, ic=ic, f0=f0, illum=illum, fw=fw)
+                        ic=ic, f0=f0, illum=illum, fw=fw)
     return rec.data, getattr(I, "data", None)
 
 
-def J_adjoint(model, src_coords, wavelet, rec_coords, recin, space_order=8,
+def J_adjoint(model, src_coords, wavelet, rec_coords, recin,
               is_residual=False, checkpointing=False, n_checkpoints=None, t_sub=1,
               return_obj=False, freq_list=[], dft_sub=None, ic="as", illum=False,
               ws=None, f0=0.015, born_fwd=False, nlind=False, misfit=None, fw=True):
@@ -315,8 +299,6 @@ def J_adjoint(model, src_coords, wavelet, rec_coords, recin, space_order=8,
         Coordiantes of the receiver(s)
     recin: Array
         Receiver data
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     checkpointing: Bool
         Whether or not to use checkpointing
     n_checkpoints: Int
@@ -345,25 +327,25 @@ def J_adjoint(model, src_coords, wavelet, rec_coords, recin, space_order=8,
     """
     if checkpointing:
         return J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin,
-                                       space_order=8, is_residual=is_residual, ws=ws,
+                                       is_residual=is_residual, ws=ws,
                                        n_checkpoints=n_checkpoints, ic=ic, f0=f0,
                                        nlind=nlind, return_obj=return_obj, illum=illum,
                                        born_fwd=born_fwd, misfit=misfit, fw=fw)
     elif freq_list is not None:
         return J_adjoint_freq(model, src_coords, wavelet, rec_coords, recin, ws=ws,
-                              space_order=space_order, dft_sub=dft_sub, f0=f0, ic=ic,
+                              dft_sub=dft_sub, f0=f0, ic=ic,
                               freq_list=freq_list, is_residual=is_residual, nlind=nlind,
                               return_obj=return_obj, misfit=misfit, born_fwd=born_fwd,
                               illum=illum, fw=fw)
     else:
         return J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin,
                                   is_residual=is_residual, ic=ic, ws=ws, t_sub=t_sub,
-                                  return_obj=return_obj, space_order=space_order,
+                                  return_obj=return_obj,
                                   born_fwd=born_fwd, f0=f0, nlind=nlind,
                                   illum=illum, misfit=misfit, fw=fw)
 
 
-def J_adjoint_freq(model, src_coords, wavelet, rec_coords, recin, space_order=8,
+def J_adjoint_freq(model, src_coords, wavelet, rec_coords, recin,
                    freq_list=[], is_residual=False, return_obj=False, nlind=False,
                    dft_sub=None, ic="as", ws=None, born_fwd=False, f0=0.015,
                    misfit=None, illum=False, fw=True):
@@ -384,8 +366,6 @@ def J_adjoint_freq(model, src_coords, wavelet, rec_coords, recin, space_order=8,
         Coordiantes of the receiver(s)
     recin: Array
         Receiver data
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     freq_list: List
         List of frequencies for on-the-fly DFT
     dft_sub: Int
@@ -415,20 +395,20 @@ def J_adjoint_freq(model, src_coords, wavelet, rec_coords, recin, space_order=8,
     """
     ffunc = op_fwd_J[born_fwd]
     rec, u, Iu, _ = ffunc(model, src_coords, rec_coords, wavelet, save=False,
-                          space_order=space_order, freq_list=freq_list, ic=ic, ws=ws,
+                          freq_list=freq_list, ic=ic, ws=ws,
                           dft_sub=dft_sub, nlind=nlind, illum=illum, f0=f0, fw=fw)
     # Residual and gradient
     f, residual = Loss(rec, recin, model.critical_dt,
                        is_residual=is_residual, misfit=misfit)
 
-    g, Iv, _ = gradient(model, residual, rec_coords, u, space_order=space_order, ic=ic,
+    g, Iv, _ = gradient(model, residual, rec_coords, u, ic=ic,
                         freq=freq_list, dft_sub=dft_sub, f0=f0, illum=illum, fw=fw)
     if return_obj:
         return f, g.data, getattr(Iu, "data", None), getattr(Iv, "data", None)
     return g.data, getattr(Iu, "data", None), getattr(Iv, "data", None)
 
 
-def J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin, space_order=8,
+def J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin,
                        is_residual=False, return_obj=False, born_fwd=False, illum=False,
                        ic="as", ws=None, t_sub=1, nlind=False, f0=0.015, misfit=None,
                        fw=True):
@@ -449,8 +429,6 @@ def J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin, space_orde
         Coordiantes of the receiver(s)
     recin: Array
         Receiver data
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     ic: String
         Imaging conditions ("as", "isic" or "fwi"), defaults to "as"
     ws : Array
@@ -476,14 +454,14 @@ def J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin, space_orde
     """
     ffunc = op_fwd_J[born_fwd]
     rec, u, Iu, _ = ffunc(model, src_coords, rec_coords, wavelet, save=True, nlind=nlind,
-                          f0=f0, ws=ws, space_order=space_order, illum=illum, ic=ic,
+                          f0=f0, ws=ws, illum=illum, ic=ic,
                           t_sub=t_sub, fw=fw)
 
     # Residual and gradient
     f, residual = Loss(rec, recin, model.critical_dt,
                        is_residual=is_residual, misfit=misfit)
 
-    g, Iv, _ = gradient(model, residual, rec_coords, u, space_order=space_order, ic=ic,
+    g, Iv, _ = gradient(model, residual, rec_coords, u, ic=ic,
                         f0=f0, illum=illum, fw=fw)
 
     if return_obj:
@@ -491,7 +469,7 @@ def J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin, space_orde
     return g.data, getattr(Iu, "data", None), getattr(Iv, "data", None)
 
 
-def J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin, space_order=8,
+def J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin,
                             is_residual=False, n_checkpoints=None, born_fwd=False,
                             return_obj=False, ic="as", ws=None, nlind=False, f0=0.015,
                             misfit=None, illum=False, fw=True):
@@ -511,8 +489,6 @@ def J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin, space
         Coordiantes of the receiver(s)
     recin: Array
         Receiver data
-    space_order: Int (optional)
-        Spatial discretization order, defaults to 8
     checkpointing: Bool
         Whether or not to use checkpointing
     n_checkpoints: Int
@@ -545,9 +521,9 @@ def J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin, space
     ffunc = op_fwd_J[born_fwd]
     # Optimal checkpointing
     op_f, u, rec_g, kwu = ffunc(model, src_coords, rec_coords, wavelet, fw=fw,
-                                save=False, space_order=space_order, return_op=True,
+                                save=False, return_op=True,
                                 ic=ic, nlind=nlind, ws=ws, f0=f0, illum=illum)
-    op, g, kwg = gradient(model, recin, rec_coords, u, space_order=space_order,
+    op, g, kwg = gradient(model, recin, rec_coords, u,
                           return_op=True, ic=ic, f0=f0, save=False, illum=illum,
                           fw=fw)
 
@@ -588,7 +564,7 @@ def J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin, space
 op_fwd_J = {False: forward, True: born}
 
 
-def wri_func(model, src_coords, wavelet, rec_coords, recin, yin, space_order=8,
+def wri_func(model, src_coords, wavelet, rec_coords, recin, yin,
              ic="as", ws=None, t_sub=1, grad="m", grad_corr=False,
              alpha_op=False, w_fun=None, eps=0, freq_list=[], wfilt=None, f0=0.015):
     """
@@ -606,7 +582,7 @@ def wri_func(model, src_coords, wavelet, rec_coords, recin, yin, space_order=8,
     # F(m0) * q if y is not an input and compute y = r(m0)
     if yin is None or grad_corr:
         y, u0, _, _ = forward(model, src_coords, rec_coords, wavelet, save=grad_corr,
-                              space_order=space_order, ws=ws, f0=f0)
+                              ws=ws, f0=f0)
         ydat = recin[:] - y.data[:]
     else:
         ydat = yin
diff --git a/src/pysource/propagators.py b/src/pysource/propagators.py
index 4481df0bd..dbaf59681 100644
--- a/src/pysource/propagators.py
+++ b/src/pysource/propagators.py
@@ -12,7 +12,7 @@
 
 
 # Forward propagation
-def forward(model, src_coords, rcv_coords, wavelet, space_order=8, save=False,
+def forward(model, src_coords, rcv_coords, wavelet, save=False,
             qwf=None, return_op=False, freq_list=None, dft_sub=None,
             norm_wf=False, w_fun=None, ws=None, wr=None, t_sub=1, f0=0.015,
             illum=False, fw=True, **kwargs):
@@ -20,6 +20,8 @@ def forward(model, src_coords, rcv_coords, wavelet, space_order=8, save=False,
     Low level propagator, to be used through `interface.py`
     Compute forward wavefield u = A(m)^{-1}*f and related quantities (u(xrcv))
     """
+    # Space order
+    space_order = model.space_order
     # Number of time steps
     nt = as_tuple(qwf)[0].shape[0] if wavelet is None else wavelet.shape[0]
 
@@ -92,12 +94,14 @@ def adjoint(*args, **kwargs):
     return forward(*args, fw=fw, **kwargs)
 
 
-def gradient(model, residual, rcv_coords, u, return_op=False, space_order=8, fw=True,
+def gradient(model, residual, rcv_coords, u, return_op=False, fw=True,
              w=None, freq=None, dft_sub=None, ic="as", f0=0.015, save=True, illum=False):
     """
     Low level propagator, to be used through `interface.py`
     Compute the action of the adjoint Jacobian onto a residual J'* δ d.
     """
+    # Space order
+    space_order = model.space_order
     # Setting adjoint wavefieldgradient
     v = wavefield(model, space_order, fw=not fw)
     try:
@@ -144,7 +148,7 @@ def gradient(model, residual, rcv_coords, u, return_op=False, space_order=8, fw=
     return gradm, I, summary
 
 
-def born(model, src_coords, rcv_coords, wavelet, space_order=8, save=False,
+def born(model, src_coords, rcv_coords, wavelet, save=False,
          qwf=None, return_op=False, ic="as", freq_list=None, dft_sub=None,
          ws=None, t_sub=1, nlind=False, f0=0.015, illum=False, fw=True):
     """
@@ -152,6 +156,8 @@ def born(model, src_coords, rcv_coords, wavelet, space_order=8, save=False,
     Compute linearized wavefield U = J(m)* δ m
     and related quantities.
     """
+    # Space order
+    space_order = model.space_order
     nt = wavelet.shape[0]
 
     # Wavefields
@@ -219,12 +225,14 @@ def born(model, src_coords, rcv_coords, wavelet, space_order=8, save=False,
 
 
 # Forward propagation
-def forward_grad(model, src_coords, rcv_coords, wavelet, v, space_order=8,
+def forward_grad(model, src_coords, rcv_coords, wavelet, v,
                  q=None, ws=None, ic="as", w=None, freq=None, f0=0.015, **kwargs):
     """
     Low level propagator, to be used through `interface.py`
     Compute forward wavefield u = A(m)^{-1}*f and related quantities (u(xrcv))
     """
+    # Space order
+    space_order = model.space_order
     # Number of time steps
     nt = as_tuple(q)[0].shape[0] if wavelet is None else wavelet.shape[0]
 
diff --git a/src/pysource/sensitivity.py b/src/pysource/sensitivity.py
index fbeb3f903..e06a0fbd8 100644
--- a/src/pysource/sensitivity.py
+++ b/src/pysource/sensitivity.py
@@ -205,8 +205,7 @@ def isic_src(model, u, **kwargs):
     ics = kwargs.get('icsign', 1)
     dus = []
     for uu in u:
-        so = uu.space_order // 2
-        du_aux = div(dm * irho * grad(uu, shift=.5, order=so), shift=-.5, order=so)
+        du_aux = div(dm * irho * grad(uu, shift=.5), shift=-.5)
         dus.append(dm * irho * uu.dt2 * m - ics * du_aux)
     if model.is_tti:
         return (-dus[0], -dus[1])
@@ -224,9 +223,7 @@ def inner_grad(u, v):
     v: TimeFunction or Function
         Second wavefield
     """
-    so = u.space_order // 2
-    return sum([a*b for a, b in zip(grad(u, order=so),
-                                    grad(v, order=so))])
+    return sum([a*b for a, b in zip(grad(u), grad(v))])
 
 
 fwi_src = lambda *ar, **kw: isic_src(*ar, icsign=-1, **kw)