--- a/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl	Mon Feb 14 09:47:35 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl	Mon Feb 14 11:55:38 2022 +0100
@@ -42,13 +42,23 @@
 end
 
 function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, inner_stencil, closure_stencils, dir)
-    return SecondDerivativeVariable{dir, dimension(grid)}(inner_stencil, Tuple(closure_stencils), size(grid), coeff)
+    Δxᵢ = spacing(grid)[dir]
+    scaled_inner_stencil = scale(inner_stencil, 1/Δxᵢ^2)
+    scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2)
+    return SecondDerivativeVariable{dir, dimension(grid)}(scaled_inner_stencil, scaled_closure_stencils, size(grid), coeff)
 end
 
 function SecondDerivativeVariable(grid::EquidistantGrid{1}, coeff::AbstractVector, inner_stencil, closure_stencils)
     return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, 1)
 end
 
+function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, stencil_set, dir)
+    inner_stencil    = parse_nested_stencil(eltype(coeff), stencil_set["D2variable"]["inner_stencil"])
+    closure_stencils = parse_nested_stencil.(eltype(coeff), stencil_set["D2variable"]["closure_stencils"])
+
+    return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, dir)
+end
+
 derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir
 
 closure_size(op::SecondDerivativeVariable) = length(op.closure_stencils)

--- a/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl	Mon Feb 14 09:47:35 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl	Mon Feb 14 11:55:38 2022 +0100
@@ -5,7 +5,7 @@
 using Sbplib.SbpOperators
 using Sbplib.SbpOperators: NestedStencil, CenteredNestedStencil
 
-
+using LinearAlgebra
 
 @testset "SecondDerivativeVariable" begin
     interior_stencil = CenteredNestedStencil((1/2, 1/2, 0.),(-1/2, -1., -1/2),( 0., 1/2, 1/2))
@@ -107,6 +107,61 @@
             @test apply_to_functions(2,v=(x,y)->x,   c=(x,y)->  1.) == zeros(11,9)
             @test apply_to_functions(2,v=(x,y)->x,   c=(x,y)-> -x ) == zeros(11,9)
             @test apply_to_functions(2,v=(x,y)->x^2, c=(x,y)->  1.) == zeros(11,9)
+
+
+
+            # TBD: This should be moved somewhere else right?
+            @testset "real operators" begin
+                c(x,y) = exp(x) + exp(1.5(1-y))
+                v(x,y) = sin(x) + cos(1.5(1-y))
+
+                Dxv(x,y) = cos(x)*exp(x) - (exp(x) + exp(1.5 - 1.5y))*sin(x)
+                Dyv(x,y) = -1.5(1.5exp(x) + 1.5exp(1.5 - 1.5y))*cos(1.5 - 1.5y) - 2.25exp(1.5 - 1.5y)*sin(1.5 - 1.5y)
+
+                # g₁ = EquidistantGrid((30,37), (0.,0.), (1.,2.))
+                g₁ = EquidistantGrid((60,67), (0.,0.), (1.,2.))
+                g₂ = refine(g₁,2)
+
+                c̄₁ = evalOn(g₁, c)
+                c̄₂ = evalOn(g₂, c)
+
+                v̄₁ = evalOn(g₁, v)
+                v̄₂ = evalOn(g₂, v)
+
+
+                function convergence_rate_estimate(D₁, D₂, Dv_true)
+                    Dv̄₁ = D₁*v̄₁
+                    Dv̄₂ = D₂*v̄₂
+
+                    Dv₁ = evalOn(g₁,Dv_true)
+                    Dv₂ = evalOn(g₂,Dv_true)
+
+                    e₁ = norm(Dv̄₁ - Dv₁)/norm(Dv₁)
+                    e₂ = norm(Dv̄₂ - Dv₂)/norm(Dv₂)
+
+                    return log2(e₁/e₂)
+                end
+
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 2)
+                Dx₁ = SecondDerivativeVariable(g₁, c̄₁, stencil_set, 1)
+                Dx₂ = SecondDerivativeVariable(g₂, c̄₂, stencil_set, 1)
+                @test convergence_rate_estimate(Dx₁, Dx₂, Dxv) ≈ 1.5 rtol = 1e-1
+
+                Dy₁ = SecondDerivativeVariable(g₁, c̄₁, stencil_set, 2)
+                Dy₂ = SecondDerivativeVariable(g₂, c̄₂, stencil_set, 2)
+                @test convergence_rate_estimate(Dy₁, Dy₂, Dyv) ≈ 1.5 rtol = 1e-1
+
+
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
+                Dx₁ = SecondDerivativeVariable(g₁, c̄₁, stencil_set, 1)
+                Dx₂ = SecondDerivativeVariable(g₂, c̄₂, stencil_set, 1)
+                @test convergence_rate_estimate(Dx₁, Dx₂, Dxv) ≈ 2.5 rtol = 1e-1
+
+                Dy₁ = SecondDerivativeVariable(g₁, c̄₁, stencil_set, 2)
+                Dy₂ = SecondDerivativeVariable(g₂, c̄₂, stencil_set, 2)
+                @test convergence_rate_estimate(Dy₁, Dy₂, Dyv) ≈ 2.5 rtol = 2e-1
+
+            end
         end
     end
 end