Function space SBP operators

.gitignore

@@ -9,3 +9,6 @@ docs/build
 docs/src/code_of_conduct.md
 docs/src/contributing.md
 docs/src/license.md
+
+run
+run/*

Project.toml

@@ -27,12 +27,14 @@ Unrolled = "9602ed7d-8fef-5bc8-8597-8f21381861e8"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 
 [extensions]
 SummationByPartsOperatorsBandedMatricesExt = "BandedMatrices"
 SummationByPartsOperatorsDiffEqCallbacksExt = "DiffEqCallbacks"
 SummationByPartsOperatorsForwardDiffExt = "ForwardDiff"
+SummationByPartsOperatorsOptimForwardDiffExt = ["Optim", "ForwardDiff"]
 SummationByPartsOperatorsStructArraysExt = "StructArrays"
 
 [compat]
@@ -46,6 +48,7 @@ InteractiveUtils = "1"
 LinearAlgebra = "1"
 LoopVectorization = "0.12.22"
 MuladdMacro = "0.2"
+Optim = "1"
 PolynomialBases = "0.4.15"
 PrecompileTools = "1.0.1"
 RecursiveArrayTools = "2.11, 3"
@@ -64,4 +67,5 @@ julia = "1.6"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"

ext/SummationByPartsOperatorsOptimForwardDiffExt.jl

@@ -0,0 +1,132 @@
+module SummationByPartsOperatorsOptimForwardDiffExt
+
+if isdefined(Base, :get_extension)
+    using Optim: Options, LBFGS, optimize, minimizer
+    using ForwardDiff: ForwardDiff
+else
+    using ..Optim: Options, LBFGS, optimize, minimizer
+    using ..ForwardDiff: ForwardDiff
+end
+
+using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023
+using LinearAlgebra: Diagonal, LowerTriangular, diag, norm, cond
+using SparseArrays: spzeros
+
+function inner_H1(f, g, f_derivative, g_derivative, nodes)
+    return sum(f.(nodes) .* g.(nodes) + f_derivative.(nodes) .* g_derivative.(nodes))
+end
+norm_H1(f, f_derivative, nodes) = sqrt(inner_H1(f, f, f_derivative, f_derivative, nodes))
+
+call_orthonormal_basis_function(A, basis_functions, k, x) = sum([basis_functions[i](x)*A[k, i] for i in 1:k])
+
+# This will orthonormalize the basis functions using the Gram-Schmidt process to reduce the condition
+# number of the Vandermonde matrix. The matrix A transfers the old basis functions to the new orthonormalized by
+# g(x) = A * f(x), where f(x) is the vector of old basis functions and g(x) is the vector of the new orthonormalized
+# basis functions. Analogously, we have g'(x) = A * f'(x).
+function orthonormalize_gram_schmidt(basis_functions, basis_functions_derivatives, nodes)
+    K = length(basis_functions)
+
+    A = LowerTriangular(zeros(eltype(nodes), K, K))
+
+    basis_functions_orthonormalized = Vector{Function}(undef, K)
+    basis_functions_orthonormalized_derivatives = Vector{Function}(undef, K)
+
+    for k = 1:K
+        A[k, k] = 1
+        for j = 1:k-1
+            g(x) = call_orthonormal_basis_function(A, basis_functions, j, x)
+            g_derivative(x) = call_orthonormal_basis_function(A, basis_functions_derivatives, j, x)
+            inner_product = inner_H1(basis_functions[k], g, basis_functions_derivatives[k], g_derivative, nodes)
+            norm_squared = inner_H1(g, g, g_derivative, g_derivative, nodes)
+            A[k, :] = A[k, :] - inner_product/norm_squared * A[j, :]
+        end
+
+        basis_functions_orthonormalized[k] = x -> call_orthonormal_basis_function(A, basis_functions, k, x)
+        basis_functions_orthonormalized_derivatives[k] = x -> call_orthonormal_basis_function(A, basis_functions_derivatives, k, x)
+        # Normalization
+        r = norm_H1(basis_functions_orthonormalized[k], basis_functions_orthonormalized_derivatives[k], nodes)
+        A[k, :] = A[k, :] / r
+    end
+    return basis_functions_orthonormalized, basis_functions_orthonormalized_derivatives
+end
+
+function vandermonde_matrix(functions, nodes)
+    N = length(nodes)
+    K = length(functions)
+    V = zeros(eltype(nodes), N, K)
+    for i in 1:N
+        for j in 1:K
+            V[i, j] = functions[j](nodes[i])
+        end
+    end
+    return V
+end
+
+function create_S(sigma, N)
+    S = zeros(eltype(sigma), N, N)
+    k = 1
+    for i in 1:N
+        for j in (i + 1):N
+            S[i, j] = sigma[k]
+            k += 1
+        end
+    end
+    return S - S'
+end
+
+sig(x) = 1 / (1 + exp(-x))
+
+function create_P(rho, x_length)
+    P = Diagonal(sig.(rho))
+    P *= x_length / sum(P)
+    return P
+end
+
+function SummationByPartsOperators.construct_function_space_operator(basis_functions, x_min, x_max, nodes,
+                                                                     ::GlaubitzNordströmÖffner2023)
+    K = length(basis_functions)
+    N = length(nodes)
+    basis_functions_derivatives = [x -> ForwardDiff.derivative(basis_functions[i], x) for i in 1:K]
+    basis_functions_orthonormalized, basis_functions_orthonormalized_derivatives = orthonormalize_gram_schmidt(basis_functions,
+                                                                                                               basis_functions_derivatives,
+                                                                                                               nodes)
+    V = vandermonde_matrix(basis_functions_orthonormalized, nodes)
+    V_x = vandermonde_matrix(basis_functions_orthonormalized_derivatives, nodes)
+    # Here, W satisfies W'*W = I
+    W = [V; -V_x]
+
+    B = spzeros(eltype(nodes), N, N)
+    B[1, 1] = -1
+    B[N, N] = 1
+
+    R = B * V / 2
+    x_length = x_max - x_min
+    p = (W, R, x_length)
+    @views function optimization_function(x, p)
+        W, R, x_length = p
+        (N, _) = size(R)
+        L = Integer(N*(N - 1)/2)
+        sigma = x[1:L]
+        rho = x[(L + 1):end]
+        S = create_S(sigma, N)
+        P = create_P(rho, x_length)
+        X = [S P]
+        # TODO: Which norm to use? XW + R is a NxK matrix. This uses the Frobenius norm.
+        return norm(X * W + R)^2
+    end
+    L = Integer(N*(N - 1)/2)
+    x0 = zeros(L + N)
+    result = optimize(x -> optimization_function(x, p), x0, LBFGS(), Options(g_tol = 1e-14, iterations = 10000),
+                      autodiff = :forward)
+    x = minimizer(result)
+    sigma = x[1:L]
+    rho = x[(L + 1):end]
+    S = create_S(sigma, N)
+    P = create_P(rho, x_length)
+    weights = diag(P)
+    Q = S + B/2
+    D = inv(P) * Q
+    return weights, D
+end
+
+end # module

src/SummationByPartsOperators.jl

@@ -100,6 +100,9 @@ function __init__()
     @require BandedMatrices="aae01518-5342-5314-be14-df237901396f" include("../ext/SummationByPartsOperatorsBandedMatricesExt.jl")
     @require DiffEqCallbacks="459566f4-90b8-5000-8ac3-15dfb0a30def" include("../ext/SummationByPartsOperatorsDiffEqCallbacksExt.jl")
     @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("../ext/SummationByPartsOperatorsForwardDiffExt.jl")
+    @require Optim="429524aa-4258-5aef-a3af-852621145aeb" begin
+      @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("../ext/SummationByPartsOperatorsOptimForwardDiffExt.jl")
+    end
     @require StructArrays="09ab397b-f2b6-538f-b94a-2f83cf4a842a" include("../ext/SummationByPartsOperatorsStructArraysExt.jl")
   end
 end
@@ -109,6 +112,7 @@ include("fourier_operators.jl")
 include("fourier_operators_2d.jl")
 include("legendre_operators.jl")
 include("matrix_operators.jl")
+include("function_space_operators.jl")
 include("upwind_operators.jl")
 include("SBP_coefficients/MattssonNordström2004.jl")
 include("SBP_coefficients/MattssonSvärdNordström2004.jl")
@@ -139,7 +143,7 @@ export PeriodicDerivativeOperator, PeriodicDissipationOperator,
        FourierPolynomialDerivativeOperator, FourierRationalDerivativeOperator,
        FourierDerivativeOperator2D,
        LegendreDerivativeOperator, LegendreSecondDerivativeOperator,
-       MatrixDerivativeOperator,
+       MatrixDerivativeOperator, FunctionSpaceOperator,
        UpwindOperators, PeriodicUpwindOperators
 export FilterCallback, ConstantFilter, ExponentialFilter
 export SafeMode, FastMode, ThreadedMode
@@ -169,6 +173,7 @@ export MattssonNordström2004, MattssonSvärdNordström2004, MattssonSvärdShoey
        SharanBradyLivescu2022
 export Tadmor1989, MadayTadmor1989, Tadmor1993,
        TadmorWaagan2012Standard, TadmorWaagan2012Convergent
+export GlaubitzNordströmÖffner2023
 
 export BurgersPeriodicSemidiscretization, BurgersNonperiodicSemidiscretization,
        CubicPeriodicSemidiscretization, CubicNonperiodicSemidiscretization,

src/function_space_operators.jl

@@ -0,0 +1,49 @@
+
+"""
+    GlaubitzNordströmÖffner2023()
+
+Function space SBP (FSBP) operators given in
+- Glaubitz, Nordström, Öffner (2023)
+  Summation-by-parts operators for general function spaces.
+  SIAM Journal on Numerical Analysis 61, 2, pp. 733-754.
+
+See also
+- Glaubitz, Nordström, Öffner (2024)
+  An optimization-based construction procedure for function space based
+  summation-by-parts operators on arbitrary grids.
+  arXiv, arXiv:2405.08770v1.
+"""
+struct GlaubitzNordströmÖffner2023 <: SourceOfCoefficients end
+
+function Base.show(io::IO, source::GlaubitzNordströmÖffner2023)
+  if get(io, :compact, false)
+    summary(io, source)
+  else
+      print(io,
+          "Glaubitz, Nordström, Öffner (2023) \n",
+          "  Summation-by-parts operators for general function spaces \n",
+          "  SIAM Journal on Numerical Analysis 61, 2, pp. 733-754. \n",
+          "See also \n",
+          "  Glaubitz, Nordström, Öffner (2024) \n",
+          "  An optimization-based construction procedure for function \n",
+          "    space based summation-by-parts operators on arbitrary grids \n",
+          "  arXiv, arXiv:2405.08770v1.")
+  end
+end
+
+"""
+    FunctionSpaceOperator(basis_functions, x_min, x_max, nodes, source; accuracy_order = 0)
+
+Construct a `MatrixDerivativeOperator` that represents a derivative operator in a function space.
+"""
+function FunctionSpaceOperator(basis_functions,
+                               x_min::T, x_max::T, nodes::Vector{T},
+                               source::SourceOfCoefficients;
+                               accuracy_order = 0) where {T, SourceOfCoefficients}
+
+  weights, D = construct_function_space_operator(basis_functions, x_min, x_max, nodes, source)
+  return MatrixDerivativeOperator(x_min, x_max, nodes, weights, D, accuracy_order, source)
+end
+
+# This function is extended in the package extension SummationByPartsOperatorsOptimExt
+function construct_function_space_operator end

test/Project.toml

@@ -4,6 +4,7 @@ BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -16,6 +17,7 @@ Aqua = "0.8"
 BandedMatrices = "0.15, 0.16, 0.17, 1"
 DiffEqCallbacks = "2, 3"
 ForwardDiff = "0.10"
+Optim = "1"
 OrdinaryDiffEq = "5, 6"
 SpecialFunctions = "1, 2"
 StaticArrays = "1.0"

test/matrix_operators_test.jl

@@ -1,6 +1,7 @@
 using Test
 using LinearAlgebra
 using SummationByPartsOperators
+using Optim
 
 # check construction of interior part of upwind operators
 @testset "Check against some upwind operators" begin
@@ -113,3 +114,38 @@ using SummationByPartsOperators
     @test SummationByPartsOperators.upper_bandwidth(Dm) == size(Dm, 1) - 1
   end
 end
+
+@testset "Function space operators" begin
+  N = 5
+  nodes = collect(range(-1, 1, length=N))
+  x_min = -1.0
+  x_max = 1.0
+  x = [x_min, -0.5, 0.0, 0.5, x_max]
+  source = GlaubitzNordströmÖffner2023()
+  for compact in (true, false)
+    show(IOContext(devnull, :compact=>compact), source)
+  end
+  B = zeros(N, N)
+  B[1, 1] = -1.0
+  B[N, N] = 1.0
+  let basis_functions = [x -> x^i for i in 0:3]
+    D = FunctionSpaceOperator(basis_functions, x_min, x_max, nodes, source)
+
+    @test ≈(D * ones(N), zeros(N); atol = 1e-13)
+    @test D * x ≈ ones(N)
+    @test D * (x .^ 2) ≈ 2 * x
+    @test D * (x .^ 3) ≈ 3 * (x .^ 2)
+    M = mass_matrix(D)
+    @test M * D.D + D.D' * M ≈ B
+  end
+
+  let basis_functions = [one, identity, exp]
+    D = FunctionSpaceOperator(basis_functions, x_min, x_max, nodes, source)
+
+    @test ≈(D * ones(N), zeros(N); atol = 1e-13)
+    @test D * x ≈ ones(N)
+    @test D * exp.(x) ≈ exp.(x)
+    M = mass_matrix(D)
+    @test M * D.D + D.D' * M ≈ B
+  end
+end