--- a/EquidistantGrid.jl	Thu Feb 14 12:46:58 2019 +0100
+++ b/EquidistantGrid.jl	Thu Feb 14 16:25:22 2019 +0100
@@ -8,17 +8,21 @@
     size::NTuple{Dim, Int} # First coordinate direction stored first
     limit_lower::NTuple{Dim, T}
     limit_upper::NTuple{Dim, T}
-    spacing::NTuple{Dim, T}
+    inverse_spacing::NTuple{Dim, T} # The reciprocal of the grid spacing
 
     # General constructor
     function EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T
         @assert all(size.>0)
         @assert all(limit_upper.-limit_lower .!= 0)
-        spacing = abs.(limit_upper.-limit_lower)./(size.-1)
-        return new{Dim,T}(size, limit_lower, limit_upper, spacing)
+        inverse_spacing = (size.-1)./abs.(limit_upper.-limit_lower)
+        return new{Dim,T}(size, limit_lower, limit_upper, inverse_spacing)
     end
 end
 
+function Base.eachindex(grid::EquidistantGrid)
+    CartesianIndices(grid.size)
+end
+
 # Returns the number of dimensions of an EquidistantGrid.
 #
 # @Input: grid - an EquidistantGrid
@@ -27,8 +31,10 @@
     return length(grid.size)
 end
 
-function Base.eachindex(grid::EquidistantGrid)
-    CartesianIndices(grid.size)
+# Returns the spacing of the grid
+#
+function spacing(grid::EquidistantGrid)
+    return 1.0./grid.inverse_spacing
 end
 
 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered

--- a/diffOp.jl	Thu Feb 14 12:46:58 2019 +0100
+++ b/diffOp.jl	Thu Feb 14 16:25:22 2019 +0100
@@ -113,13 +113,13 @@
     return uᵢ
 end
 
-function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
+@inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
     # 2nd x-derivative
     @inbounds vx = view(v, :, Int(I[2]))
-    @inbounds uᵢ = L.a*apply(L.op, L.grid.spacing[1], vx , I[1])
+    @inbounds uᵢ = L.a*apply(L.op, L.grid.inverse_spacing[1], vx , I[1])
     # 2nd y-derivative
     @inbounds vy = view(v, Int(I[1]), :)
-    @inbounds uᵢ += L.a*apply(L.op, L.grid.spacing[2], vy, I[2])
+    @inbounds uᵢ += L.a*apply(L.op, L.grid.inverse_spacing[2], vy, I[2])
     return uᵢ
 end
 

--- a/sbpD2.jl	Thu Feb 14 12:46:58 2019 +0100
+++ b/sbpD2.jl	Thu Feb 14 16:25:22 2019 +0100
@@ -2,16 +2,16 @@
 
 # Apply for different regions Lower/Interior/Upper or Unknown region
 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Lower})
-    return @inbounds apply(op.closureStencils[Int(i)], v, Int(i))/h^2
+    return @inbounds h*h*apply(op.closureStencils[Int(i)], v, Int(i))
 end
 
 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Interior})
-    return @inbounds apply(op.innerStencil, v, Int(i))/h^2
+    return @inbounds h*h*apply(op.innerStencil, v, Int(i))
 end
 
 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Upper})
     N = length(v)
-    return @inbounds Int(op.parity)*apply_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))/h^2
+    return @inbounds h*h*Int(op.parity)*apply_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
 end
 
 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, index::Index{Unknown})

--- a/stencil.jl	Thu Feb 14 12:46:58 2019 +0100
+++ b/stencil.jl	Thu Feb 14 16:25:22 2019 +0100
@@ -17,13 +17,14 @@
 end
 
 Base.@propagate_inbounds @inline function apply(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
-    w = zero(eltype(v))
-    @simd for k ∈ 1:N
+    w = s.weights[1]*v[i+ s.range[1]]
+    @simd for k ∈ 2:N
         w += s.weights[k]*v[i+ s.range[1] + k-1]
     end
     return w
 end
 
+# TODO: Fix loop unrolling here as well. Then we can also remove Base.getindex(::Stencil)
 Base.@propagate_inbounds @inline function apply_backwards(s::Stencil, v::AbstractVector, i::Int)
     w = zero(eltype(v))
     for j ∈ s.range[2]:-1:s.range[1]