--- a/src/SbpOperators/SbpOperators.jl	Tue Mar 15 21:40:31 2022 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Tue Mar 15 21:42:52 2022 +0100
@@ -7,6 +7,7 @@
 export Laplace
 export laplace
 export normal_derivative
+export first_derivative
 export second_derivative
 
 using Sbplib.RegionIndices
@@ -24,6 +25,7 @@
 include("readoperator.jl")
 include("volumeops/volume_operator.jl")
 include("volumeops/constant_interior_scaling_operator.jl")
+include("volumeops/derivatives/first_derivative.jl")
 include("volumeops/derivatives/second_derivative.jl")
 include("volumeops/derivatives/second_derivative_variable.jl")
 include("volumeops/laplace/laplace.jl")

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Tue Mar 15 21:42:52 2022 +0100
@@ -0,0 +1,31 @@
+"""
+    first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
+
+Creates the first-derivative operator `D1` as a `TensorMapping`
+
+`D1` approximates the first-derivative d/dξ on `grid` along the coordinate dimension specified by
+`direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
+for the points in the closure regions.
+
+On a one-dimensional `grid`, `D1` is a `VolumeOperator`. On a multi-dimensional `grid`, `D1` is the outer product of the
+one-dimensional operator with the `IdentityMapping`s in orthogonal coordinate dirrections.
+
+See also: [`volume_operator`](@ref).
+"""
+function first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
+    h_inv = inverse_spacing(grid)[direction]
+    return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd, direction)
+end
+first_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = first_derivative(grid,inner_stencil,closure_stencils,1)
+
+"""
+    first_derivative(grid, stencil_set, direction)
+
+Creates a `first_derivative` operator on `grid` along coordinate dimension `direction` given a parsed TOML
+`stencil_set`.
+"""
+function first_derivative(grid::EquidistantGrid, stencil_set, direction)
+    inner_stencil = parse_stencil(stencil_set["D1"]["inner_stencil"])
+    closure_stencils = parse_stencil.(stencil_set["D1"]["closure_stencils"])
+    first_derivative(grid,inner_stencil,closure_stencils,direction);
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Tue Mar 15 21:42:52 2022 +0100
@@ -0,0 +1,81 @@
+using Test
+
+
+using Sbplib.SbpOperators
+using Sbplib.Grids
+using Sbplib.LazyTensors
+
+using Sbplib.SbpOperators: closure_size, Stencil
+
+"""
+    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 "first_derivative" begin
+    @testset "Constructors" begin
+        stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+
+        g₁ = EquidistantGrid(11, 0., 1.)
+        g₂ = EquidistantGrid((11,14), (0.,1.), (1.,3.))
+
+        @test first_derivative(g₁, stencil_set, 1) isa TensorMapping{Float64,1,1}
+        @test first_derivative(g₂, stencil_set, 2) isa TensorMapping{Float64,2,2}
+
+        interior_stencil = CenteredStencil(-1,0,1)
+        closure_stencils = [Stencil(-1,1, center=1)]
+
+        @test first_derivative(g₁, interior_stencil, closure_stencils, 1) isa TensorMapping{Float64,1,1}
+        @test first_derivative(g₂, interior_stencil, closure_stencils, 2) isa TensorMapping{Float64,2,2}
+    end
+
+    @testset "Accuracy conditions" begin
+        N = 20
+        g = EquidistantGrid(N, 0//1,2//1)
+        @testset for order ∈ [2,4]
+            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
+            D₁ = first_derivative(g, stencil_set, 1)
+
+            @testset "boundary x^$k" for k ∈ 0:order÷2
+                v = evalOn(g, x->monomial(x,k))
+
+                @testset for i ∈ 1:closure_size(D₁)
+                    x, = points(g)[i]
+                    @test (D₁*v)[i] == monomial(x,k-1)
+                end
+
+                @testset for i ∈ (N-closure_size(D₁)+1):N
+                    x, = points(g)[i]
+                    @test (D₁*v)[i] == monomial(x,k-1)
+                end
+            end
+
+            @testset "interior x^$k" for k ∈ 0:order
+                v = evalOn(g, x->monomial(x,k))
+
+                x, = points(g)[10]
+                @test (D₁*v)[10] == monomial(x,k-1)
+            end
+        end
+    end
+
+    @testset "Accuracy on function" begin
+        g = EquidistantGrid(30, 0.,1.)
+        v = evalOn(g, x->exp(x))
+        @testset for (order, tol) ∈ [(2, 6e-3),(4, 2e-4)]
+            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
+            D₁ = first_derivative(g, stencil_set, 1)
+
+            @test D₁*v ≈ v rtol=tol
+        end
+    end
+end
+