make scale_by_mass_matrix! and scale_by_inverse_mass_matrix! publ…

…ic API (#278)

* export scale_by_mass_matrix! and scale_by_inverse_mass_matrix!

* implement more methods for mass matrix scaling

* also for upwind operators

* mention in docs

* fix

Project.toml

@@ -1,7 +1,7 @@
 name = "SummationByPartsOperators"
 uuid = "9f78cca6-572e-554e-b819-917d2f1cf240"
 author = ["Hendrik Ranocha"]
-version = "0.5.62"
+version = "0.5.63"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"

docs/src/introduction.md

@@ -332,7 +332,9 @@ are necessary.
   are often needed. These can be obtained without forming `M = mass_matrix(D)`
   explicitly via [`left_boundary_weight`](@ref) and [`right_boundary_weight`](@ref).
 - Instead of forming a mass matrix explicitly, discrete integrals can be evaluated
-  efficiently using [`integrate`](@ref).
+  efficiently using [`integrate`](@ref) and vectors can be scaled by the mass matrix
+  or its inverse using [`scale_by_mass_matrix!`](@ref) and
+  [`scale_by_inverse_mass_matrix!`](@ref), respectively.
 - Dissipation operators based on the same discrete inner product as SBP derivative
   operators can be obtained via [`dissipation_operator`](@ref).
 

src/SBP_operators.jl

@@ -655,6 +655,48 @@ function integrate(func, u::AbstractVector, D::DerivativeOperator)
     Δx * res
 end
 
+function scale_by_mass_matrix!(u::AbstractVector, D::DerivativeOperator)
+    Base.require_one_based_indexing(u)
+    @boundscheck begin
+        length(u) == length(grid(D))
+    end
+    @unpack Δx = D
+    @unpack left_weights, right_weights = D.coefficients
+
+    @inbounds for i in (1 + length(left_weights)):(length(u) - length(right_weights))
+        u[i] *= Δx
+    end
+    @inbounds for i in Base.OneTo(length(left_weights))
+        u[i] *= left_weights[i] * Δx
+    end
+    @inbounds for i in Base.OneTo(length(right_weights))
+        u[end-i+1] *= right_weights[i] * Δx
+    end
+
+    return u
+end
+
+function scale_by_inverse_mass_matrix!(u::AbstractVector, D::DerivativeOperator)
+    Base.require_one_based_indexing(u)
+    @boundscheck begin
+        length(u) == length(grid(D))
+    end
+    @unpack Δx = D
+    @unpack left_weights, right_weights = D.coefficients
+
+    @inbounds for i in (1 + length(left_weights)):(length(u) - length(right_weights))
+        u[i] /= Δx
+    end
+    @inbounds for i in Base.OneTo(length(left_weights))
+        u[i] /= left_weights[i] * Δx
+    end
+    @inbounds for i in Base.OneTo(length(right_weights))
+        u[end-i+1] /= right_weights[i] * Δx
+    end
+
+    return u
+end
+
 
 """
     derivative_operator(source_of_coefficients,

src/SummationByPartsOperators.jl

@@ -147,6 +147,7 @@ export SafeMode, FastMode, ThreadedMode
 export derivative_order, accuracy_order, source_of_coefficients, grid, semidiscretize
 export mass_matrix
 export integrate, left_boundary_weight, right_boundary_weight,
+       scale_by_mass_matrix!, scale_by_inverse_mass_matrix!,
        derivative_left, derivative_right,
        mul_transpose_derivative_left!, mul_transpose_derivative_right!,
        evaluate_coefficients, evaluate_coefficients!,

src/fourier_operators.jl

@@ -149,6 +149,18 @@ function mass_matrix(D::FourierDerivativeOperator)
     Δx * I
 end
 
+function scale_by_mass_matrix!(u::AbstractVector, D::FourierDerivativeOperator)
+    @unpack Δx = D
+
+    u .*= Δx
+end
+
+function scale_by_inverse_mass_matrix!(u::AbstractVector, D::FourierDerivativeOperator)
+    @unpack Δx = D
+
+    u ./= Δx
+end
+
 
 
 """

src/general_operators.jl

@@ -73,6 +73,22 @@ function Base.summary(io::IO, D::AbstractDerivativeOperator)
               ", accuracy:", accuracy_order(D), ")")
 end
 
+"""
+    scale_by_mass_matrix!(u, D)
+
+Scale the vector `u` by the mass matrix associated to the
+derivative operator `D`.
+"""
+function scale_by_mass_matrix! end
+
+"""
+    scale_by_inverse_mass_matrix!(u, D)
+
+Scale the vector `u` by the inverse of the mass matrix associated to the
+derivative operator `D`.
+"""
+function scale_by_inverse_mass_matrix! end
+
 
 
 abstract type AbstractExecutionMode end

src/periodic_operators.jl

@@ -624,6 +624,18 @@ function mass_matrix(D::PeriodicDerivativeOperator)
     Δx * I
 end
 
+function scale_by_mass_matrix!(u::AbstractVector, D::PeriodicDerivativeOperator)
+    @unpack Δx = D
+
+    u .*= Δx
+end
+
+function scale_by_inverse_mass_matrix!(u::AbstractVector, D::PeriodicDerivativeOperator)
+    @unpack Δx = D
+
+    u ./= Δx
+end
+
 
 """
     periodic_central_derivative_operator(derivative_order, accuracy_order,

src/upwind_operators.jl

@@ -133,6 +133,13 @@ function integrate(func, u::AbstractVector, D::Union{UpwindOperators,PeriodicUpw
   integrate(func, u, D.minus)
 end
 
+function scale_by_mass_matrix!(u::AbstractVector, D::Union{UpwindOperators,PeriodicUpwindOperators})
+  scale_by_mass_matrix!(u, D.minus)
+end
+
+function scale_by_inverse_mass_matrix!(u::AbstractVector, D::Union{UpwindOperators,PeriodicUpwindOperators})
+  scale_by_inverse_mass_matrix!(u, D.minus)
+end
 
 """
     upwind_operators(source_type, args...; derivative_order = 1, kwargs...)

test/SBP_operators_test.jl

@@ -65,6 +65,14 @@ for source in D_test_list, T in (Float32,Float64)
         # boundary derivative
         @test derivative_left( D, x1, Val{0}()) ≈ x1[begin]
         @test derivative_right(D, x1, Val{0}()) ≈ x1[end]
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     acc_order = 4
@@ -122,6 +130,14 @@ for source in D_test_list, T in (Float32,Float64)
         # boundary derivative
         @test derivative_left( D, x1, Val{0}()) ≈ x1[begin]
         @test derivative_right(D, x1, Val{0}()) ≈ x1[end]
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     acc_order = 6
@@ -185,6 +201,14 @@ for source in D_test_list, T in (Float32,Float64)
         # boundary derivative
         @test derivative_left( D, x1, Val{0}()) ≈ x1[begin]
         @test derivative_right(D, x1, Val{0}()) ≈ x1[end]
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     acc_order = 8
@@ -246,6 +270,14 @@ for source in D_test_list, T in (Float32,Float64)
         # boundary derivative
         @test derivative_left( D, x1, Val{0}()) ≈ x1[begin]
         @test derivative_right(D, x1, Val{0}()) ≈ x1[end]
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     @test_throws Union{MethodError,ArgumentError} derivative_operator(source, der_order, 16, xmin, xmax, N)
@@ -308,6 +340,14 @@ for source in D_test_list, T in (Float32,Float64)
         @test derivative_right(D, x1, Val{1}()) ≈ one(T)
         @test derivative_left( D, x2, Val{1}()) ≈ 2xmin
         @test derivative_right(D, x2, Val{1}()) ≈ 2xmax
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     acc_order = 4
@@ -367,6 +407,14 @@ for source in D_test_list, T in (Float32,Float64)
         @test derivative_right(D, x2, Val{1}()) ≈ 2xmax
         @test derivative_left( D, x3, Val{1}()) ≈ 3xmin^2
         @test derivative_right(D, x3, Val{1}()) ≈ 3xmax^2
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     acc_order = 6
@@ -430,6 +478,14 @@ for source in D_test_list, T in (Float32,Float64)
         @test derivative_right(D, x3, Val{1}()) ≈ 3xmax^2
         @test derivative_left( D, x4, Val{1}()) ≈ 4xmin^3
         @test derivative_right(D, x4, Val{1}()) ≈ 4xmax^3
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     acc_order = 8
@@ -498,6 +554,14 @@ for source in D_test_list, T in (Float32,Float64)
         @test derivative_right(D, x4, Val{1}()) ≈ 4xmax^3
         @test derivative_left( D, x5, Val{1}()) ≈ 5xmin^4
         @test derivative_right(D, x5, Val{1}()) ≈ 5xmax^4
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     @test_throws Union{MethodError,ArgumentError} derivative_operator(source, der_order, 16, xmin, xmax, N)
@@ -556,6 +620,14 @@ end
         @test abs(derivative_right(D, x1, Val{2}())) < 600 * eps(T)
         @test derivative_left( D, x2, Val{2}()) ≈ 2
         @test derivative_right(D, x2, Val{2}()) ≈ 2
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     acc_order = 4
@@ -613,6 +685,14 @@ end
         @test derivative_right(D, x2, Val{2}()) ≈ 2
         @test derivative_left( D, x3, Val{2}()) ≈ 6xmin
         @test derivative_right(D, x3, Val{2}()) ≈ 6xmax
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     acc_order = 6
@@ -675,6 +755,14 @@ end
         @test isapprox(derivative_right(D, x2, Val{2}()), 2, atol=20_000*eps(T))
         @test isapprox(derivative_left( D, x3, Val{2}()), 6xmin, atol=50_000*eps(T))
         @test isapprox(derivative_right(D, x3, Val{2}()), 6xmax, atol=80_000*eps(T))
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     @test_throws Union{MethodError,ArgumentError} derivative_operator(source, der_order, 16, xmin, xmax, N)
@@ -742,6 +830,14 @@ end
         @test abs(derivative_right(D, x1, Val{3}())) < 10 * eps(T) / D.Δx^3
         @test abs(derivative_left( D, x2, Val{3}())) < 10 * eps(T) / D.Δx^3
         @test abs(derivative_right(D, x2, Val{3}())) < 10 * eps(T) / D.Δx^3
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     acc_order = 4
@@ -809,6 +905,14 @@ end
         @test abs(derivative_right(D, x2, Val{3}())) < 100_000*eps(T)
         @test isapprox(derivative_left( D, x3, Val{3}()), 6, atol=150_000*eps(T))
         @test isapprox(derivative_right(D, x3, Val{3}()), 6, atol=1_000_000*eps(T))
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     acc_order = 6
@@ -876,6 +980,14 @@ end
         @test abs(derivative_right(D, x2, Val{3}())) < 1_000_000*eps(T)
         @test isapprox(derivative_left( D, x3, Val{3}()), 6, atol=900_000*eps(T))
         @test isapprox(derivative_right(D, x3, Val{3}()), 6, atol=5_000_000*eps(T))
+        # mass matrix scaling
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 
     @test_throws Union{MethodError,ArgumentError} derivative_operator(source, der_order, 16, xmin, xmax, N)

test/fourier_operators_test.jl

@@ -40,6 +40,16 @@ for T in (Float32, Float64)
             @test maximum(abs, duplot - dufunc.(xplot)) < 5N*eps(T)
             @test abs(integrate(u, D)) < N*eps(T)
         end
+
+        # mass matrix scaling
+        x1 = grid(D)
+        M = @inferred mass_matrix(D)
+        u = sinpi.(x1)
+        v = copy(u)
+        scale_by_mass_matrix!(v, D)
+        @test v ≈ M * u
+        scale_by_inverse_mass_matrix!(v, D)
+        @test v ≈ u
     end
 end