diff -r 4c17a7d6ae5e -r ba023fc09961 src/SbpOperators/SbpOperators.jl
--- a/src/SbpOperators/SbpOperators.jl	Tue Mar 15 06:53:06 2022 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Tue Mar 15 07:29:41 2022 +0100
@@ -7,6 +7,7 @@
 export Laplace
 export laplace
 export normal_derivative
+export first_derivative
 export second_derivative
 
 using Sbplib.RegionIndices
@@ -22,6 +23,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/laplace/laplace.jl")
 include("volumeops/inner_products/inner_product.jl")

diff -r 4c17a7d6ae5e -r ba023fc09961 src/SbpOperators/volumeops/derivatives/first_derivative.jl
--- a/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Tue Mar 15 06:53:06 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Tue Mar 15 07:29:41 2022 +0100
@@ -1,5 +1,3 @@
-export first_derivative
-
 """
     first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
 
@@ -12,11 +10,22 @@
 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: [`SbpOperators.volume_operator`](@ref).
+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, closure_stencils) = first_derivative(grid,inner_stencil,closure_stencils,1)
+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

diff -r 4c17a7d6ae5e -r ba023fc09961 test/SbpOperators/volumeops/derivatives/first_derivative_test.jl
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Tue Mar 15 06:53:06 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Tue Mar 15 07:29:41 2022 +0100
@@ -3,3 +3,50 @@
 using Sbplib.SbpOperators
 using Sbplib.Grids
 
+using Sbplib.SbpOperators: closure_size
+
+"""
+    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 "accuracy" 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 accuracy $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 accuracy $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
+end
+