From 0c50276e8d9f6f18d03ea01573adf29e4588ec2e Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Fri, 14 Jun 2024 16:12:44 +0100 Subject: [PATCH] Function space SBP operators (#268) * first draft for function space operators * some spaces * order alphabetically * don't force basis_function to be passed as vector * allow passing accuracy_order * first (basically) working version * orthonormalization via Gram-Schmidt * add some basic tests * some smaller optimizations * fix import of Options * another fix... * avoid unnecessary anonymous function * use version of spzeros that is compatible with older Julia versions * FunctionSpaceOperator -> function_space_operator * add proper docstring * extend docstring * put constructor function in extension * allow passing options * clarify that the operator is of derivative_order 1 * only Optim has to be loaded * only for Julia 1.9 * remove fallback imports with Requires * use PreallocationTools.jl to reduce allocations * add kwarg derivative-order * add analytical gradient (still much slower than ForwardDiff.jl) * fuse function and gradient evaluation * avoid some type conversions * use div instead of sum * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * change order of loops * Update ext/SummationByPartsOperatorsOptimForwardDiffExt.jl Co-authored-by: Hendrik Ranocha * import dot * move optimization_function_and_grad to fix type instability introduced by closure bug * move optimization_function_and_grad! outside (that was it!) * Update ext/SummationByPartsOperatorsOptimForwardDiffExt.jl Co-authored-by: Hendrik Ranocha * some cleanups * ignore stale dependency PreallocationTools in Aqua * add test with non-equidistant nodes * fix test * remove kwargs x_min and x_max * add experimental feature warning * Update test/matrix_operators_test.jl Co-authored-by: Hendrik Ranocha * remove l_hat2 since it this case is impossible * remove PreallocationTools again * throw ArgumentError instead of assertion for wrong derivative_order --------- Co-authored-by: Hendrik Ranocha Co-authored-by: Hendrik Ranocha --- .gitignore | 3 + Project.toml | 4 + ...tionByPartsOperatorsOptimForwardDiffExt.jl | 210 ++++++++++++++++++ src/SummationByPartsOperators.jl | 4 +- src/function_space_operators.jl | 63 ++++++ test/Project.toml | 2 + test/matrix_operators_test.jl | 60 +++++ 7 files changed, 345 insertions(+), 1 deletion(-) create mode 100644 ext/SummationByPartsOperatorsOptimForwardDiffExt.jl create mode 100644 src/function_space_operators.jl diff --git a/.gitignore b/.gitignore index e33a4292..75ac9fe8 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,6 @@ docs/build docs/src/code_of_conduct.md docs/src/contributing.md docs/src/license.md + +run +run/* diff --git a/Project.toml b/Project.toml index 287659d2..5299f448 100644 --- a/Project.toml +++ b/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" diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl new file mode 100644 index 00000000..b4801d41 --- /dev/null +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -0,0 +1,210 @@ +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 SparseArrays: spzeros + +function SummationByPartsOperators.function_space_operator(basis_functions, nodes::Vector{T}, + source::SourceOfCoefficients; + derivative_order = 1, accuracy_order = 0, + options = Options(g_tol = 1e-14, iterations = 10000)) where {T, SourceOfCoefficients} + + if derivative_order != 1 + throw(ArgumentError("Derivative order $derivative_order not implemented.")) + end + sort!(nodes) + weights, D = construct_function_space_operator(basis_functions, nodes, source; options = options) + return MatrixDerivativeOperator(first(nodes), last(nodes), 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, 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 = last(nodes) - first(nodes) + 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) + daij_dsigmak = zeros(eltype(nodes), N, K, L) + daij_drhok = zeros(eltype(nodes), N, K, N) + p = (V, V_x, R, x_length, S, A, SV, PV_x, 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, A, SV, PV_x, daij_dsigmak, daij_drhok = p + (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_hat = N - int_D_plus_one_half + if 1 <= l_hat <= i - 1 + daij_dsigmak[i, j, k] -= V[l_hat, 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 diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index 38955bee..52a6c2eb 100644 --- a/src/SummationByPartsOperators.jl +++ b/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, diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl new file mode 100644 index 00000000..081cd167 --- /dev/null +++ b/src/function_space_operators.jl @@ -0,0 +1,63 @@ + +""" + 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, 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`, where `x_min` is taken as the minimal value in +`nodes` and `x_max` the maximal value. Note that the `nodes` will be sorted internally. 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. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. +""" +function function_space_operator end diff --git a/test/Project.toml b/test/Project.toml index 4c20ba58..29c0f20e 100644 --- a/test/Project.toml +++ b/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" diff --git a/test/matrix_operators_test.jl b/test/matrix_operators_test.jl index 3a517a40..f7a31092 100644 --- a/test/matrix_operators_test.jl +++ b/test/matrix_operators_test.jl @@ -1,6 +1,7 @@ using Test using LinearAlgebra using SummationByPartsOperators +using Optim: Optim # to enable loading the function space operator optimization code # check construction of interior part of upwind operators @testset "Check against some upwind operators" begin @@ -113,3 +114,62 @@ using SummationByPartsOperators @test SummationByPartsOperators.upper_bandwidth(Dm) == size(Dm, 1) - 1 end end + +if VERSION >= v"1.9" + @testset "Function space operators" begin + N = 5 + x_min = -1.0 + x_max = 1.0 + nodes = collect(range(x_min, x_max, length=N)) + 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 = function_space_operator(basis_functions, nodes, source) + # Only first-derivative operators are implemented yet + @test_throws ArgumentError function_space_operator(basis_functions, nodes, source; derivative_order = 2) + + @test grid(D) ≈ nodes + @test all(isapprox.(D * ones(N), zeros(N); atol = 1e-13)) + @test D * nodes ≈ ones(N) + @test D * (nodes .^ 2) ≈ 2 * nodes + @test D * (nodes .^ 3) ≈ 3 * (nodes .^ 2) + M = mass_matrix(D) + @test M * D.D + D.D' * M ≈ B + end + + let basis_functions = [one, identity, exp] + D = function_space_operator(basis_functions, nodes, source) + + @test grid(D) ≈ nodes + @test all(isapprox.(D * ones(N), zeros(N); atol = 1e-13)) + @test D * nodes ≈ ones(N) + @test D * exp.(nodes) ≈ exp.(nodes) + M = mass_matrix(D) + @test M * D.D + D.D' * M ≈ B + end + + # test non-equidistant nodes generated by `nodes = [0.0, rand(8)..., 1.0]` + nodes = [0.0, 0.01585580467018155, 0.18010381213204507, 0.270467434432868, + 0.37699483985320303, 0.5600831197666554, 0.5698824835924449, 0.623949064816263, + 0.8574665549914025, 1.0] + N = length(nodes) + B = zeros(N, N) + B[1, 1] = -1.0 + B[N, N] = 1.0 + let basis_functions = [one, identity, exp] + D = function_space_operator(basis_functions, nodes, source) + + @test grid(D) ≈ nodes + @test all(isapprox.(D * ones(N), zeros(N); atol = 1e-11)) + @test D * nodes ≈ ones(N) + @test D * exp.(nodes) ≈ exp.(nodes) + M = mass_matrix(D) + @test M * D.D + D.D' * M ≈ B + end + end +end