--- a/src/SbpOperators/SbpOperators.jl	Mon Mar 21 10:04:15 2022 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Mon Mar 21 15:12:59 2022 +0100
@@ -9,6 +9,7 @@
 export normal_derivative
 export first_derivative
 export second_derivative
+export dissipation
 
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors
@@ -25,6 +26,7 @@
 include("volumeops/constant_interior_scaling_operator.jl")
 include("volumeops/derivatives/first_derivative.jl")
 include("volumeops/derivatives/second_derivative.jl")
+include("volumeops/derivatives/dissipation.jl")
 include("volumeops/laplace/laplace.jl")
 include("volumeops/inner_products/inner_product.jl")
 include("volumeops/inner_products/inverse_inner_product.jl")

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl	Mon Mar 21 15:12:59 2022 +0100
@@ -0,0 +1,61 @@
+function dissipation(g::EquidistantGrid, p, direction)
+    h_inv = inverse_spacing(g)[direction]
+
+    # D = volume_operator(g,CenteredStencil(1),(CenteredStencil(1)), )
+    return nothing, nothing
+end
+
+dissipation(g::EquidistantGrid{1}, p) = dissipation(g, p, 1)
+
+function dissipation_interior_weights(p)
+   if p == 0
+       return (1,)
+   end
+
+   return (0, dissipation_interior_weights(p-1)...) .- (dissipation_interior_weights(p-1)..., 0)
+end
+
+function dissipation_interior_stencil(p)
+    w = dissipation_interior_weights(p)
+    Stencil(w..., center=midpoint(w))
+end
+
+function dissipation_transpose_interior_stencil(p)
+    w = dissipation_interior_weights(p)
+    Stencil(w..., center=midpoint_transpose(w))
+end
+
+
+midpoint(weights) = length(weights)÷2 + 1
+midpoint_transpose(weights) = length(weights)+1 - midpoint(weights)
+
+dissipation_lower_closure_size(weights) = midpoint(weights) - 1
+dissipation_upper_closure_size(weights) = length(weights) - midpoint(weights)
+
+dissipation_lower_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=i                       ), dissipation_lower_closure_size(interior_weights))
+dissipation_upper_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=length(interior_weights)-dissipation_upper_closure_size(interior_weights)+i), dissipation_upper_closure_size(interior_weights))
+
+dissipation_transpose_lower_closure_stencils(interior_weights) =         ntuple(i->dissipation_transpose_lower_closure_stencil(interior_weights, i), length(interior_weights))
+dissipation_transpose_upper_closure_stencils(interior_weights) = reverse(ntuple(i->dissipation_transpose_upper_closure_stencil(interior_weights, i), length(interior_weights)))
+
+
+function dissipation_transpose_lower_closure_stencil(interior_weights, i)
+    w = ntuple(k->interior_weights[i], dissipation_lower_closure_size(interior_weights))
+
+    for k ∈ i:-1:1
+        w = (w..., interior_weights[k])
+    end
+
+    return Stencil(w..., center = i)
+end
+
+function dissipation_transpose_upper_closure_stencil(interior_weights, i)
+    j = length(interior_weights)+1-i
+    w = ntuple(k->interior_weights[j], dissipation_upper_closure_size(interior_weights))
+
+    for k ∈ j:1:length(interior_weights)
+        w = (interior_weights[k], w...)
+    end
+
+    return Stencil(w..., center = length(interior_weights)-midpoint(interior_weights)+1)
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Mon Mar 21 15:12:59 2022 +0100
@@ -0,0 +1,189 @@
+using Test
+
+using Sbplib.SbpOperators
+using Sbplib.Grids
+using Sbplib.LazyTensors
+
+using Sbplib.SbpOperators: Stencil
+
+using Sbplib.SbpOperators: dissipation_interior_weights
+using Sbplib.SbpOperators: dissipation_interior_stencil, dissipation_transpose_interior_stencil
+using Sbplib.SbpOperators: midpoint, midpoint_transpose
+using Sbplib.SbpOperators: dissipation_lower_closure_size, dissipation_upper_closure_size
+using Sbplib.SbpOperators: dissipation_lower_closure_stencils,dissipation_upper_closure_stencils
+using Sbplib.SbpOperators: dissipation_transpose_lower_closure_stencils, dissipation_transpose_upper_closure_stencils
+
+"""
+    monomial(x,k)
+
+Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``.
+Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)``
+"""
+function monomial(x,k)
+    if k < 0
+        return zero(x)
+    end
+    x^k/factorial(k)
+end
+
+@testset "dissipation" begin
+    g = EquidistantGrid(20, 0., 11.)
+    D,Dᵀ = dissipation(g, 1)
+
+    @test D isa LazyTensor{Float64,1,1} where T
+    @test Dᵀ isa LazyTensor{Float64,1,1} where T
+
+     @testset "Accuracy conditions" begin
+        N = 20
+        g = EquidistantGrid(N, 0//1,2//1)
+        @testset "D_$p" for p ∈ [1,2,3,4]
+            D,Dᵀ = dissipation(g, p)
+
+            @testset "x^$k" for k ∈ 0:1
+                v = evalOn(g, x->monomial(x,k))
+
+                x, = points(g)[10]
+                @test (D*v)[10] == monomial(x,k-1)
+            end
+
+            # Test Dᵀ works backwards and interior forwards
+        end
+    end
+end
+
+@testset "dissipation_interior_weights" begin
+    @test dissipation_interior_weights(1) == (-1, 1)
+    @test dissipation_interior_weights(2) == (1,-2, 1)
+    @test dissipation_interior_weights(3) == (-1, 3,-3, 1)
+    @test dissipation_interior_weights(4) == (1, -4, 6, -4, 1)
+end
+
+@testset "dissipation_interior_stencil" begin
+    @test dissipation_interior_stencil(1) == Stencil(-1,1, center=2)
+    @test dissipation_interior_stencil(2) == Stencil(1,-2,1, center=2)
+    @test dissipation_interior_stencil(3) == Stencil(-1,3,-3,1, center=3)
+    @test dissipation_interior_stencil(4) == Stencil(1, -4, 6, -4, 1, center=3)
+end
+
+@testset "dissipation_transpose_interior_stencil" begin
+    @test dissipation_transpose_interior_stencil(1) == Stencil(-1,1, center=1)
+    @test dissipation_transpose_interior_stencil(2) == Stencil(1,-2,1, center=2)
+    @test dissipation_transpose_interior_stencil(3) == Stencil(-1,3,-3,1, center=2)
+    @test dissipation_transpose_interior_stencil(4) == Stencil(1, -4, 6, -4, 1, center=3)
+end
+
+@testset "midpoint" begin
+    @test midpoint((1,1)) == 2
+    @test midpoint((1,1,1)) == 2
+    @test midpoint((1,1,1,1)) == 3
+    @test midpoint((1,1,1,1,1)) == 3
+end
+
+@testset "midpoint_transpose" begin
+    @test midpoint_transpose((1,1)) == 1
+    @test midpoint_transpose((1,1,1)) == 2
+    @test midpoint_transpose((1,1,1,1)) == 2
+    @test midpoint_transpose((1,1,1,1,1)) == 3
+end
+
+@testset "dissipation_lower_closure_size" begin
+    @test dissipation_lower_closure_size((1,1)) == 1
+    @test dissipation_lower_closure_size((1,1,1)) == 1
+    @test dissipation_lower_closure_size((1,1,1,1)) == 2
+    @test dissipation_lower_closure_size((1,1,1,1,1)) == 2
+end
+
+@testset "dissipation_upper_closure_size" begin
+    @test dissipation_upper_closure_size((1,1)) == 0
+    @test dissipation_upper_closure_size((1,1,1)) == 1
+    @test dissipation_upper_closure_size((1,1,1,1)) == 1
+    @test dissipation_upper_closure_size((1,1,1,1,1)) == 2
+end
+
+@testset "dissipation_lower_closure_stencils" begin
+    cases = (
+        (-1,1) => (
+            Stencil(-1, 1, center=1),
+        ),
+        (1,-2,1) => (
+            Stencil( 1,-2, 1, center=1),
+        ),
+        (-1,3,-3,1) => (
+            Stencil(-1,3,-3,1, center=1),
+            Stencil(-1,3,-3,1, center=2),
+        ),
+        (1, -4, 6, -4, 1) => (
+            Stencil(1, -4, 6, -4, 1, center=1),
+            Stencil(1, -4, 6, -4, 1, center=2),
+        )
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_lower_closure_stencils(w) == closure_stencils
+    end
+end
+
+@testset "dissipation_upper_closure_stencils" begin
+    cases = (
+        (-1,1) => (),
+        (1,-2,1) => (
+            Stencil( 1,-2, 1, center=3),
+        ),
+        (-1,3,-3,1) => (
+            Stencil(-1,3,-3,1, center=4),
+        ),
+        (1, -4, 6, -4, 1) => (
+            Stencil(1, -4, 6, -4, 1, center=4),
+            Stencil(1, -4, 6, -4, 1, center=5),
+        )
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_upper_closure_stencils(w) == closure_stencils
+    end
+end
+
+
+@testset "dissipation_transpose_lower_closure_stencils" begin
+    cases = (
+        (-1,1) => (
+            Stencil(-1,-1,    center=1),
+            Stencil( 1, 1,-1, center=2),
+        ),
+        (1,-2,1) => (
+            Stencil( 1, 1,    center=1),
+            Stencil(-2,-2, 1,  center=2),
+            Stencil( 1, 1,-2, 1, center=3),
+        ),
+        (-1,3,-3,1) => (
+            Stencil(-1,-1,-1,          center=1),
+            Stencil( 3, 3, 3,-1,       center=2),
+            Stencil(-3,-3,-3, 3,-1,    center=3),
+            Stencil( 1, 1, 1,-3, 3,-1, center=4),
+        ),
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_transpose_lower_closure_stencils(w) == closure_stencils
+    end
+end
+
+@testset "dissipation_transpose_upper_closure_stencils" begin
+    cases = (
+        (-1,1) => (
+            Stencil( 1,-1, center = 1),
+            Stencil(    1, center = 1),
+        ),
+        (1,-2,1) => (
+            Stencil( 1, -2, 1, 1, center=2),
+            Stencil(     1,-2,-2, center=2),
+            Stencil(        1, 1, center=2),
+        ),
+        (-1,3,-3,1) => (
+            Stencil( 1,-3, 3,-1,-1, center=2),
+            Stencil(    1,-3, 3, 3, center=2),
+            Stencil(       1,-3,-3, center=2),
+            Stencil(          1, 1, center=2),
+        ),
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_transpose_upper_closure_stencils(w) == closure_stencils
+    end
+end