--- a/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl	Fri Feb 18 07:19:08 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl	Fri Feb 18 07:39:07 2022 +0100
@@ -31,6 +31,26 @@
     return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, 1)
 end
 
+@doc raw"""
+    SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, stencil_set, dir)
+
+Create a `TensorMapping` for the second derivative with a variable coefficient
+`coeff` on `grid` from the stencils in `stencil_set`. The direction is
+determined by `dir`.
+
+`coeff` is a grid function on `grid`.
+
+# Example
+With
+```
+D = SecondDerivativeVariable(g, c, stencil_set, 2)
+```
+then `D*u` approximates
+```math
+\frac{\partial}{\partial y} c(x,y) \frac{\partial u}{\partial y},
+```
+on ``(0,1)⨯(0,1)`` represented by `g`.
+"""
 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"])

--- a/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl	Fri Feb 18 07:19:08 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl	Fri Feb 18 07:39:07 2022 +0100
@@ -123,9 +123,6 @@
             @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?
-            # REVIEW: Better name for testset
             @testset "standard diagonal operators" begin
                 c(x,y) = exp(x) + exp(1.5(1-y))
                 v(x,y) = sin(x) + cos(1.5(1-y))
@@ -133,7 +130,6 @@
                 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)