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

@@ -12,6 +12,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 PolynomialBases = "c74db56a-226d-5e98-8bb0-a6049094aeea"
+PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -27,12 +28,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,7 +49,9 @@ InteractiveUtils = "1"
 LinearAlgebra = "1"
 LoopVectorization = "0.12.22"
 MuladdMacro = "0.2"
+Optim = "1"
 PolynomialBases = "0.4.15"
+PreallocationTools = "0.4"
 PrecompileTools = "1.0.1"
 RecursiveArrayTools = "2.11, 3"
 Reexport = "0.2, 1.0"
@@ -64,4 +69,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,222 @@
+module SummationByPartsOperatorsOptimForwardDiffExt
+
+using Optim: Optim, Options, LBFGS, optimize, minimizer
+using ForwardDiff: ForwardDiff
+
+using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023, MatrixDerivativeOperator
+using LinearAlgebra: Diagonal, LowerTriangular, dot, diag, norm, mul!
+using PreallocationTools: DiffCache, get_tmp
+using SparseArrays: spzeros
+
+function SummationByPartsOperators.function_space_operator(basis_functions,
+                                                           x_min::T, x_max::T, nodes::Vector{T},
+                                                           source::SourceOfCoefficients;
+                                                           derivative_order = 1, accuracy_order = 0,
+                                                           options = Options(g_tol = 1e-14, iterations = 10000)) where {T, SourceOfCoefficients}
+
+    @assert derivative_order == 1 "Only first derivative operators are supported"
+    sort!(nodes)
+    weights, D = construct_function_space_operator(basis_functions, x_min, x_max, nodes, source; options = options)
+    return MatrixDerivativeOperator(x_min, x_max, nodes, weights, D, accuracy_order, source)
+end
+
+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)
+    set_S!(S, sigma, N)
+    return S
+end
+
+function set_S!(S, sigma, N)
+    k = 1
+    for i in 1:N
+        for j in (i + 1):N
+            S[i, j] = sigma[k]
+            S[j, i] = -sigma[k]
+            k += 1
+        end
+    end
+end
+
+sig(x) = 1 / (1 + exp(-x))
+sig_deriv(x) = sig(x) * (1 - sig(x))
+
+function create_P(rho, x_length)
+    P = Diagonal(sig.(rho))
+    P *= x_length / sum(P)
+    return P
+end
+
+function construct_function_space_operator(basis_functions, x_min, x_max, nodes,
+                                           ::GlaubitzNordströmÖffner2023;
+                                           options = Options(g_tol = 1e-14, iterations = 10000))
+    K = length(basis_functions)
+    N = length(nodes)
+    L = div(N * (N - 1), 2)
+    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
+    S = zeros(eltype(nodes), N, N)
+    A = zeros(eltype(nodes), N, K)
+    SV = zeros(eltype(nodes), N, K)
+    PV_x = zeros(eltype(nodes), N, K)
+    S_cache = DiffCache(S)
+    A_cache = DiffCache(A)
+    SV_cache = DiffCache(SV)
+    PV_x_cache = DiffCache(PV_x)
+    daij_dsigmak = zeros(eltype(nodes), N, K, L)
+    daij_drhok = zeros(eltype(nodes), N, K, N)
+    p = (V, V_x, R, x_length, S_cache, A_cache, SV_cache, PV_x_cache, daij_dsigmak, daij_drhok)
+
+    x0 = zeros(L + N)
+    fg!(F, G, x) = optimization_function_and_grad!(F, G, x, p)
+    result = optimize(Optim.only_fg!(fg!), x0, LBFGS(), options)
+
+    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
+
+@views function optimization_function_and_grad!(F, G, x, p)
+    V, V_x, R, x_length, S_cache, A_cache, SV_cache, PV_x_cache, daij_dsigmak, daij_drhok = p
+    S = get_tmp(S_cache, x)
+    A = get_tmp(A_cache, x)
+    SV = get_tmp(SV_cache, x)
+    PV_x = get_tmp(PV_x_cache, x)
+    (N, _) = size(R)
+    L = div(N * (N - 1), 2)
+    sigma = x[1:L]
+    rho = x[(L + 1):end]
+    set_S!(S, sigma, N)
+    P = create_P(rho, x_length)
+    mul!(SV, S, V)
+    mul!(PV_x, P, V_x)
+    @. A = SV - PV_x + R
+    if !isnothing(G)
+        fill!(daij_dsigmak, zero(eltype(daij_dsigmak)))
+        for k in axes(daij_dsigmak, 3)
+            for j in axes(daij_dsigmak, 2)
+                for i in axes(daij_dsigmak, 1)
+                    l_tilde = k + i - N * (i - 1) + div(i * (i - 1), 2)
+                    # same as above, but needs more type conversions
+                    # l_tilde = Int(k + i - (i - 1) * (N - i/2))
+                    if i + 1 <= l_tilde <= N
+                        daij_dsigmak[i, j, k] += V[l_tilde, j]
+                    else
+                        C = N^2 - 3*N + 2*i - 2*k + 1/4
+                        if C >= 0
+                            D = sqrt(C)
+                            D_plus_one_half = D + 0.5
+                            D_plus_one_half_trunc = trunc(D_plus_one_half)
+                            if D_plus_one_half == D_plus_one_half_trunc
+                                int_D_plus_one_half = trunc(Int, D_plus_one_half_trunc)
+                                l_hat1 = N + int_D_plus_one_half - 1
+                                l_hat2 = N - int_D_plus_one_half
+                                if 1 <= l_hat1 <= i - 1
+                                    daij_dsigmak[i, j, k] -= V[l_hat1, j]
+                                end
+                                if 1 <= l_hat2 <= i - 1
+                                    daij_dsigmak[i, j, k] -= V[l_hat2, j]
+                                end
+                            end
+                        end
+                    end
+                end
+            end
+        end
+        sig_rho = sig.(rho)
+        sig_deriv_rho = sig_deriv.(rho)
+        sum_sig_rho = sum(sig_rho)
+        for k in axes(daij_drhok, 3)
+            for j in axes(daij_drhok, 2)
+                for i in axes(daij_drhok, 1)
+                    factor1 = x_length * V_x[i, j] / sum_sig_rho^2
+                    factor =  factor1 * sig_deriv_rho[k]
+                    if k == i
+                        daij_drhok[i, j, k] = -factor * (sum_sig_rho - sig_rho[k])
+                    else
+                        daij_drhok[i, j, k] = factor * sig_rho[i]
+                    end
+                end
+            end
+        end
+        for k in axes(daij_dsigmak, 3)
+            G[k] = 2 * dot(daij_dsigmak[:, :, k], A)
+        end
+        for k in axes(daij_drhok, 3)
+            G[L + k] = 2 * dot(daij_drhok[:, :, k], A)
+        end
+    end
+    if !isnothing(F)
+        return norm(A)^2
+    end
+end
+
+end # module

src/SummationByPartsOperators.jl

@@ -109,6 +109,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")
@@ -154,7 +155,7 @@ export periodic_central_derivative_operator, periodic_derivative_operator, deriv
        dissipation_operator, var_coef_derivative_operator,
        fourier_derivative_operator,
        legendre_derivative_operator, legendre_second_derivative_operator,
-       upwind_operators
+       upwind_operators, function_space_operator
 export UniformMesh1D, UniformPeriodicMesh1D
 export couple_continuously, couple_discontinuously
 export mul!
@@ -169,6 +170,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,60 @@
+
+"""
+    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.
+
+See [`function_space_operator`](@ref).
+"""
+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
+
+# This function is extended in the package extension SummationByPartsOperatorsOptimExt
+"""
+    function_space_operator(basis_functions, x_min, x_max, nodes, source;
+                            derivative_order = 1, accuracy_order = 0,
+                            options = Optim.Options(g_tol = 1e-14, iterations = 10000))
+
+Construct an operator that represents a first-derivative operator in a function space spanned by
+the `basis_functions`, which is an iterable of functions. The operator is constructed on the
+interval `[x_min, x_max]` with the nodes `nodes`. The `accuracy_order` is the order of the
+accuracy of the operator, which can optionally be passed, but does not have any effect on the
+operator. The operator is constructed solving an optimization problem with Optim.jl. You
+can specify the options for the optimization problem with the `options` argument, see also
+the [documentation of Optim.jl](https://julianlsolvers.github.io/Optim.jl/stable/user/config/).
+
+The operator that is returned follows the general interface. Currently, it is wrapped in a
+[`MatrixDerivativeOperator`](@ref), but this might change in the future.
+In order to use this function, the package `Optim` must be loaded.
+
+See also [`GlaubitzNordströmÖffner2023`](@ref).
+
+!!! compat "Julia 1.9"
+    This function requires at least Julia 1.9.
+"""
+function 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/aqua.jl

@@ -4,5 +4,5 @@ using SummationByPartsOperators
 Aqua.test_all(SummationByPartsOperators;
   ambiguities = false, # a lot of false positives from dependencies
   unbound_args = false, # TODO: a strange problem I do not understand right now
-  stale_deps = (; ignore = [:PrecompileTools]),
+  stale_deps = (; ignore = [:PrecompileTools, :PreallocationTools]),
 )