--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Sun Jun 09 17:04:01 2024 -0700
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Sun Jun 09 17:51:57 2024 -0700
@@ -53,14 +53,13 @@
 end
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
 
-
 """
-sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; tuning)
+    sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
 
 The operators required to construct the SAT for imposing a Dirichlet condition.
-`tuning` specifies the strength of the penalty. See
+`H_tuning` and `R_tuning` are used to specify the strength of the penalty.
 
-See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref).
+See also: [`sat`](@ref),[`DirichletCondition`](@ref),[`positivity_decomposition`](@ref).
 """
 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
     id = boundary(bc)
@@ -75,12 +74,11 @@
 end
 
 """
-sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
+    sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
 
 The operators required to construct the SAT for imposing a Neumann condition
 
-
-See also: [`sat`,`NeumannCondition`](@ref).
+See also: [`sat`](@ref),[`NeumannCondition`](@ref).
 """
 function sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
     id = boundary(bc)
@@ -94,13 +92,22 @@
     return penalty_tensor, d
 end
 
+"""
+    positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
 
+Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive definite with respect to
+the boundary quadrature. Here `d` is the normal derivative and `e` is the boundary restriction operator.
+`B` can then be used to form a symmetric and energy stable penalty for a Dirichlet condition.
+The parameters `H_tuning` and `R_tuning` are used to specify the strength of the penalty and
+must be greater than 1. For details we refer to https://doi.org/10.1016/j.jcp.2020.109294
+"""
 function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
+    @assert(H_tuning ≥ 1.)
+    @assert(R_tuning ≥ 1.)
     Nτ_H, τ_R = positivity_limits(Δ,g,bc)
     return H_tuning*Nτ_H + R_tuning*τ_R
 end
 
-
 # TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
 # change bc::BoundaryCondition to id::BoundaryIdentifier
 function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition)