From d888354722700a4269eb9cb10c1361cb04eaaa80 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 24 May 2024 19:01:16 +0200 Subject: [PATCH 01/45] first draft for function space operators --- .gitignore | 3 ++ Project.toml | 4 ++ ext/SummationByPartsOperatorsOptimExt.jl | 18 +++++++++ src/SummationByPartsOperators.jl | 5 ++- src/function_space_operators.jl | 49 ++++++++++++++++++++++++ 5 files changed, 78 insertions(+), 1 deletion(-) create mode 100644 ext/SummationByPartsOperatorsOptimExt.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..10db5ed7 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" +SummationByPartsOperatorsOptimExt = "Optim" 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" @@ -65,3 +68,4 @@ BandedMatrices = "aae01518-5342-5314-be14-df237901396f" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" diff --git a/ext/SummationByPartsOperatorsOptimExt.jl b/ext/SummationByPartsOperatorsOptimExt.jl new file mode 100644 index 00000000..fe050f74 --- /dev/null +++ b/ext/SummationByPartsOperatorsOptimExt.jl @@ -0,0 +1,18 @@ +module SummationByPartsOperatorsOptimExt + +if isdefined(Base, :get_extension) + using Optim: LBFGS, optimize +else + using ..Optim: LBFGS, optimize +end + +using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023 + +function SummationByPartsOperators.construct_function_space_operator(basis, x_min, x_max, nodes, ::GlaubitzNordströmÖffner2023) + N = length(nodes) + weights = ones(N) + D = zeros(N, N) + return weights, D +end + +end # module diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index 38955bee..023a8d35 100644 --- a/src/SummationByPartsOperators.jl +++ b/src/SummationByPartsOperators.jl @@ -101,6 +101,7 @@ function __init__() @require DiffEqCallbacks="459566f4-90b8-5000-8ac3-15dfb0a30def" include("../ext/SummationByPartsOperatorsDiffEqCallbacksExt.jl") @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("../ext/SummationByPartsOperatorsForwardDiffExt.jl") @require StructArrays="09ab397b-f2b6-538f-b94a-2f83cf4a842a" include("../ext/SummationByPartsOperatorsStructArraysExt.jl") + @require Optim="429524aa-4258-5aef-a3af-852621145aeb" include("../ext/SummationByPartsOperatorsOptimExt.jl") end end @@ -109,6 +110,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 +141,7 @@ export PeriodicDerivativeOperator, PeriodicDissipationOperator, FourierPolynomialDerivativeOperator, FourierRationalDerivativeOperator, FourierDerivativeOperator2D, LegendreDerivativeOperator, LegendreSecondDerivativeOperator, - MatrixDerivativeOperator, + MatrixDerivativeOperator, FunctionSpaceOperator, UpwindOperators, PeriodicUpwindOperators export FilterCallback, ConstantFilter, ExponentialFilter export SafeMode, FastMode, ThreadedMode @@ -169,6 +171,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..22fbdb84 --- /dev/null +++ b/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", + " 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", + " space based summation-by-parts operators on arbitrary grids \n", + " arXiv, arXiv:2405.08770v1.") + end +end + +""" + FunctionSpaceOperator(basis, x_min, x_max, nodes, source) + +Construct a `MatrixDerivativeOperator` that represents a derivative operator in a function space. +""" +function FunctionSpaceOperator(basis::Vector{BasisFunctions}, + x_min::T, x_max::T, nodes::Vector{T}, + source::SourceOfCoefficients) where {BasisFunctions, T, SourceOfCoefficients} + + weights, D = construct_function_space_operator(basis, x_min, x_max, nodes, source) + accuracy_order = length(basis) # TODO: this might needs to be adjusted + 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 From 7fe321759802bf60b0a840fa3abedeae13f7a5ab Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 24 May 2024 19:07:27 +0200 Subject: [PATCH 02/45] some spaces --- ext/SummationByPartsOperatorsOptimExt.jl | 1 + src/function_space_operators.jl | 6 +++--- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimExt.jl b/ext/SummationByPartsOperatorsOptimExt.jl index fe050f74..064aa6a9 100644 --- a/ext/SummationByPartsOperatorsOptimExt.jl +++ b/ext/SummationByPartsOperatorsOptimExt.jl @@ -9,6 +9,7 @@ end using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023 function SummationByPartsOperators.construct_function_space_operator(basis, x_min, x_max, nodes, ::GlaubitzNordströmÖffner2023) + # TODO: implement the construction of the function space operator by solving an optimization problem N = length(nodes) weights = ones(N) D = zeros(N, N) diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index 22fbdb84..a084656e 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -21,12 +21,12 @@ function Base.show(io::IO, source::GlaubitzNordströmÖffner2023) else print(io, "Glaubitz, Nordström, Öffner (2023) \n", - " Summation-by-parts operators for general function spaces", + " 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", - " space based summation-by-parts operators on arbitrary grids \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 From 72e1522cb90df77c95fb9074835ff0ce05e61308 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 24 May 2024 19:24:05 +0200 Subject: [PATCH 03/45] order alphabetically --- Project.toml | 2 +- src/SummationByPartsOperators.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 10db5ed7..3b3b7535 100644 --- a/Project.toml +++ b/Project.toml @@ -67,5 +67,5 @@ julia = "1.6" BandedMatrices = "aae01518-5342-5314-be14-df237901396f" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" Optim = "429524aa-4258-5aef-a3af-852621145aeb" +StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index 023a8d35..4bb2535e 100644 --- a/src/SummationByPartsOperators.jl +++ b/src/SummationByPartsOperators.jl @@ -100,8 +100,8 @@ 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 StructArrays="09ab397b-f2b6-538f-b94a-2f83cf4a842a" include("../ext/SummationByPartsOperatorsStructArraysExt.jl") @require Optim="429524aa-4258-5aef-a3af-852621145aeb" include("../ext/SummationByPartsOperatorsOptimExt.jl") + @require StructArrays="09ab397b-f2b6-538f-b94a-2f83cf4a842a" include("../ext/SummationByPartsOperatorsStructArraysExt.jl") end end From 5f85d67239eba859010d270eab29886c8062a380 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Sat, 25 May 2024 14:11:43 +0200 Subject: [PATCH 04/45] don't force basis_function to be passed as vector --- ext/SummationByPartsOperatorsOptimExt.jl | 2 +- src/function_space_operators.jl | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimExt.jl b/ext/SummationByPartsOperatorsOptimExt.jl index 064aa6a9..aa952aee 100644 --- a/ext/SummationByPartsOperatorsOptimExt.jl +++ b/ext/SummationByPartsOperatorsOptimExt.jl @@ -8,7 +8,7 @@ end using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023 -function SummationByPartsOperators.construct_function_space_operator(basis, x_min, x_max, nodes, ::GlaubitzNordströmÖffner2023) +function SummationByPartsOperators.construct_function_space_operator(basis_functions, x_min, x_max, nodes, ::GlaubitzNordströmÖffner2023) # TODO: implement the construction of the function space operator by solving an optimization problem N = length(nodes) weights = ones(N) diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index a084656e..64c2bc1f 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -32,16 +32,16 @@ function Base.show(io::IO, source::GlaubitzNordströmÖffner2023) end """ - FunctionSpaceOperator(basis, x_min, x_max, nodes, source) + FunctionSpaceOperator(basis_functions, x_min, x_max, nodes, source) Construct a `MatrixDerivativeOperator` that represents a derivative operator in a function space. """ -function FunctionSpaceOperator(basis::Vector{BasisFunctions}, +function FunctionSpaceOperator(basis_functions, x_min::T, x_max::T, nodes::Vector{T}, - source::SourceOfCoefficients) where {BasisFunctions, T, SourceOfCoefficients} + source::SourceOfCoefficients) where {T, SourceOfCoefficients} - weights, D = construct_function_space_operator(basis, x_min, x_max, nodes, source) - accuracy_order = length(basis) # TODO: this might needs to be adjusted + weights, D = construct_function_space_operator(basis_functions, x_min, x_max, nodes, source) + accuracy_order = length(basis_functions) # TODO: this might needs to be adjusted return MatrixDerivativeOperator(x_min, x_max, nodes, weights, D, accuracy_order, source) end From f7a00768a4977d9ce76dbf83202396b8f923b8fa Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Sat, 25 May 2024 14:14:39 +0200 Subject: [PATCH 05/45] allow passing accuracy_order --- src/function_space_operators.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index 64c2bc1f..6c1357ee 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -32,16 +32,16 @@ function Base.show(io::IO, source::GlaubitzNordströmÖffner2023) end """ - FunctionSpaceOperator(basis_functions, x_min, x_max, nodes, source) + 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) where {T, SourceOfCoefficients} + source::SourceOfCoefficients; + accuracy_order = 0) where {T, SourceOfCoefficients} weights, D = construct_function_space_operator(basis_functions, x_min, x_max, nodes, source) - accuracy_order = length(basis_functions) # TODO: this might needs to be adjusted return MatrixDerivativeOperator(x_min, x_max, nodes, weights, D, accuracy_order, source) end From 944046bb1c614b16733dce70b439c682533a10d8 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 10:59:53 +0200 Subject: [PATCH 06/45] first (basically) working version --- Project.toml | 2 +- ext/SummationByPartsOperatorsOptimExt.jl | 19 ---- ...tionByPartsOperatorsOptimForwardDiffExt.jl | 89 +++++++++++++++++++ src/SummationByPartsOperators.jl | 4 +- 4 files changed, 93 insertions(+), 21 deletions(-) delete mode 100644 ext/SummationByPartsOperatorsOptimExt.jl create mode 100644 ext/SummationByPartsOperatorsOptimForwardDiffExt.jl diff --git a/Project.toml b/Project.toml index 3b3b7535..5299f448 100644 --- a/Project.toml +++ b/Project.toml @@ -34,7 +34,7 @@ StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" SummationByPartsOperatorsBandedMatricesExt = "BandedMatrices" SummationByPartsOperatorsDiffEqCallbacksExt = "DiffEqCallbacks" SummationByPartsOperatorsForwardDiffExt = "ForwardDiff" -SummationByPartsOperatorsOptimExt = "Optim" +SummationByPartsOperatorsOptimForwardDiffExt = ["Optim", "ForwardDiff"] SummationByPartsOperatorsStructArraysExt = "StructArrays" [compat] diff --git a/ext/SummationByPartsOperatorsOptimExt.jl b/ext/SummationByPartsOperatorsOptimExt.jl deleted file mode 100644 index aa952aee..00000000 --- a/ext/SummationByPartsOperatorsOptimExt.jl +++ /dev/null @@ -1,19 +0,0 @@ -module SummationByPartsOperatorsOptimExt - -if isdefined(Base, :get_extension) - using Optim: LBFGS, optimize -else - using ..Optim: LBFGS, optimize -end - -using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023 - -function SummationByPartsOperators.construct_function_space_operator(basis_functions, x_min, x_max, nodes, ::GlaubitzNordströmÖffner2023) - # TODO: implement the construction of the function space operator by solving an optimization problem - N = length(nodes) - weights = ones(N) - D = zeros(N, N) - return weights, D -end - -end # module diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl new file mode 100644 index 00000000..afd2907b --- /dev/null +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -0,0 +1,89 @@ +module SummationByPartsOperatorsOptimForwardDiffExt + +if isdefined(Base, :get_extension) + using Optim: Options, LBFGS, optimize, minimizer + using ForwardDiff: ForwardDiff +else + using ..Optim: LBFGS, optimize, minimizer + using ..ForwardDiff: ForwardDiff +end + +using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023 +using LinearAlgebra: Diagonal, diag, norm, cond + +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 + +function create_P(rho, x_length) + sig(x) = 1 / (1 + exp(-x)) + 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] + # TODO: Orthogonalize basis w.r.t. H1 inner product via Gram-Schmidt + V = vandermonde_matrix(basis_functions, nodes) + V_x = vandermonde_matrix(basis_functions_derivatives, nodes) + W = [V; -V_x] + # display(cond(W)) + + B = zeros((N, N)) + B[1, 1] = -1.0 + B[N, N] = 1.0 + + R = B * V / 2 + x_length = x_max - x_min + p = (W, R, x_length) + 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) + # display(result) + 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 diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index 4bb2535e..994d74c1 100644 --- a/src/SummationByPartsOperators.jl +++ b/src/SummationByPartsOperators.jl @@ -100,7 +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" include("../ext/SummationByPartsOperatorsOptimExt.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 From 18ce61261df661ae8aa7cd2cb6c7a95d5d230f68 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 16:06:37 +0200 Subject: [PATCH 07/45] orthonormalization via Gram-Schmidt --- ...tionByPartsOperatorsOptimForwardDiffExt.jl | 62 ++++++++++++++++--- 1 file changed, 52 insertions(+), 10 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index afd2907b..c4b4567e 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -5,11 +5,49 @@ if isdefined(Base, :get_extension) using ForwardDiff: ForwardDiff else using ..Optim: LBFGS, optimize, minimizer - using ..ForwardDiff: ForwardDiff + using ..ForwardDiff: Options, ForwardDiff end using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023 -using LinearAlgebra: Diagonal, diag, norm, cond +using LinearAlgebra: Diagonal, LowerTriangular, diag, norm, cond + +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{eltype(basis_functions)}(undef, K) + basis_functions_orthonormalized_derivatives = Vector{eltype(basis_functions)}(undef, K) + + for k = 1:K + A[k, k] = 1.0 + 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) @@ -35,22 +73,26 @@ function create_S(sigma, N) return S - S' end +sig(x) = 1 / (1 + exp(-x)) + function create_P(rho, x_length) - sig(x) = 1 / (1 + exp(-x)) 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) +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] - # TODO: Orthogonalize basis w.r.t. H1 inner product via Gram-Schmidt - V = vandermonde_matrix(basis_functions, nodes) - V_x = vandermonde_matrix(basis_functions_derivatives, nodes) + 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] - # display(cond(W)) B = zeros((N, N)) B[1, 1] = -1.0 @@ -73,8 +115,8 @@ function SummationByPartsOperators.construct_function_space_operator(basis_funct 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) - # display(result) + 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] From b576ece0cc319abf28e245016b491146c1979e11 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 17:06:25 +0200 Subject: [PATCH 08/45] add some basic tests --- ...tionByPartsOperatorsOptimForwardDiffExt.jl | 4 +-- test/Project.toml | 2 ++ test/matrix_operators_test.jl | 36 +++++++++++++++++++ 3 files changed, 40 insertions(+), 2 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index c4b4567e..e0618be5 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -27,8 +27,8 @@ function orthonormalize_gram_schmidt(basis_functions, basis_functions_derivative A = LowerTriangular(zeros(eltype(nodes), K, K)) - basis_functions_orthonormalized = Vector{eltype(basis_functions)}(undef, K) - basis_functions_orthonormalized_derivatives = Vector{eltype(basis_functions)}(undef, K) + basis_functions_orthonormalized = Vector{Function}(undef, K) + basis_functions_orthonormalized_derivatives = Vector{Function}(undef, K) for k = 1:K A[k, k] = 1.0 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..46134a9a 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 # 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 From f3f720853a6aaacf497fb13deb2db412dcd0b40d Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 17:19:05 +0200 Subject: [PATCH 09/45] some smaller optimizations --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index e0618be5..dc6a906a 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -10,6 +10,7 @@ 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)) @@ -31,7 +32,7 @@ function orthonormalize_gram_schmidt(basis_functions, basis_functions_derivative basis_functions_orthonormalized_derivatives = Vector{Function}(undef, K) for k = 1:K - A[k, k] = 1.0 + 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) @@ -94,14 +95,14 @@ function SummationByPartsOperators.construct_function_space_operator(basis_funct # Here, W satisfies W'*W = I W = [V; -V_x] - B = zeros((N, N)) - B[1, 1] = -1.0 - B[N, N] = 1.0 + 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) - function optimization_function(x, p) + @views function optimization_function(x, p) W, R, x_length = p (N, _) = size(R) L = Integer(N*(N - 1)/2) From cd92e5394502acf3074e66d4407a682afe22d5d3 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 17:46:32 +0200 Subject: [PATCH 10/45] fix import of Options --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index dc6a906a..5fb78cb0 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -4,8 +4,8 @@ if isdefined(Base, :get_extension) using Optim: Options, LBFGS, optimize, minimizer using ForwardDiff: ForwardDiff else - using ..Optim: LBFGS, optimize, minimizer - using ..ForwardDiff: Options, ForwardDiff + using ..Optim: LOptions, BFGS, optimize, minimizer + using ..ForwardDiff: ForwardDiff end using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023 From 8c77cd5571667f8f296f171c1e00ab3c97eab951 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 17:47:56 +0200 Subject: [PATCH 11/45] another fix... --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 5fb78cb0..beff963f 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -4,7 +4,7 @@ if isdefined(Base, :get_extension) using Optim: Options, LBFGS, optimize, minimizer using ForwardDiff: ForwardDiff else - using ..Optim: LOptions, BFGS, optimize, minimizer + using ..Optim: Options, LBFGS, optimize, minimizer using ..ForwardDiff: ForwardDiff end From d62d73b94772bb39c4a10a67867a36d81a1d057c Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 17:58:50 +0200 Subject: [PATCH 12/45] avoid unnecessary anonymous function --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index beff963f..b8ac29a4 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -34,8 +34,8 @@ function orthonormalize_gram_schmidt(basis_functions, basis_functions_derivative 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) + 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, :] From 9d7cc4a9d903746f23f7fefcc3ce41f6bd2c541b Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 18:18:24 +0200 Subject: [PATCH 13/45] use version of spzeros that is compatible with older Julia versions --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index b8ac29a4..abf0bc4d 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -35,7 +35,7 @@ function orthonormalize_gram_schmidt(basis_functions, basis_functions_derivative 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) + 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, :] @@ -95,7 +95,7 @@ function SummationByPartsOperators.construct_function_space_operator(basis_funct # Here, W satisfies W'*W = I W = [V; -V_x] - B = spzeros(eltype(nodes), (N, N)) + B = spzeros(eltype(nodes), N, N) B[1, 1] = -1 B[N, N] = 1 From a2ef507933aa4964d4264d1a86427f0242d519b9 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 20:45:49 +0200 Subject: [PATCH 14/45] FunctionSpaceOperator -> function_space_operator --- src/SummationByPartsOperators.jl | 4 ++-- src/function_space_operators.jl | 10 +++++----- test/matrix_operators_test.jl | 4 ++-- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index 994d74c1..78e6ddaf 100644 --- a/src/SummationByPartsOperators.jl +++ b/src/SummationByPartsOperators.jl @@ -143,7 +143,7 @@ export PeriodicDerivativeOperator, PeriodicDissipationOperator, FourierPolynomialDerivativeOperator, FourierRationalDerivativeOperator, FourierDerivativeOperator2D, LegendreDerivativeOperator, LegendreSecondDerivativeOperator, - MatrixDerivativeOperator, FunctionSpaceOperator, + MatrixDerivativeOperator, UpwindOperators, PeriodicUpwindOperators export FilterCallback, ConstantFilter, ExponentialFilter export SafeMode, FastMode, ThreadedMode @@ -158,7 +158,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! diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index 6c1357ee..fe42228e 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -32,14 +32,14 @@ function Base.show(io::IO, source::GlaubitzNordströmÖffner2023) end """ - FunctionSpaceOperator(basis_functions, x_min, x_max, nodes, source; accuracy_order = 0) + function_space_operator(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} +function function_space_operator(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) diff --git a/test/matrix_operators_test.jl b/test/matrix_operators_test.jl index 46134a9a..803a50d0 100644 --- a/test/matrix_operators_test.jl +++ b/test/matrix_operators_test.jl @@ -129,7 +129,7 @@ end 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) + D = function_space_operator(basis_functions, x_min, x_max, nodes, source) @test ≈(D * ones(N), zeros(N); atol = 1e-13) @test D * x ≈ ones(N) @@ -140,7 +140,7 @@ end end let basis_functions = [one, identity, exp] - D = FunctionSpaceOperator(basis_functions, x_min, x_max, nodes, source) + D = function_space_operator(basis_functions, x_min, x_max, nodes, source) @test ≈(D * ones(N), zeros(N); atol = 1e-13) @test D * x ≈ ones(N) From b4b4c7d24f16c11627d15ddf2329f1100e9f6603 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 20:52:57 +0200 Subject: [PATCH 15/45] add proper docstring --- src/function_space_operators.jl | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index fe42228e..b3250712 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -12,6 +12,8 @@ See also 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 @@ -34,7 +36,15 @@ end """ function_space_operator(basis_functions, x_min, x_max, nodes, source; accuracy_order = 0) -Construct a `MatrixDerivativeOperator` that represents a derivative operator in a function space. +Construct an operator that represents a 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 that is returned follows the general interface. Currently, it is wrapped in a +[`MatrixDerivativeOperator`](@ref), but this might change in the future. + +See also [`GlaubitzNordströmÖffner2023`](@ref). """ function function_space_operator(basis_functions, x_min::T, x_max::T, nodes::Vector{T}, From 1b67472f663e8382bbf48a5d15f37dc3706c9c58 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 20:54:17 +0200 Subject: [PATCH 16/45] extend docstring --- src/function_space_operators.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index b3250712..5ab1bd25 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -43,6 +43,7 @@ accuracy of the operator, which can optionally be passed, but does not have any operator. 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 packages `Optim` and `ForwardDiff` must be loaded. See also [`GlaubitzNordströmÖffner2023`](@ref). """ From d4266fe62f3306c0912913be1259b2e1512f9a52 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 20:59:00 +0200 Subject: [PATCH 17/45] put constructor function in extension --- ...ummationByPartsOperatorsOptimForwardDiffExt.jl | 15 ++++++++++++--- src/function_space_operators.jl | 12 ++---------- 2 files changed, 14 insertions(+), 13 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index abf0bc4d..cdd4f7da 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -8,10 +8,19 @@ else using ..ForwardDiff: ForwardDiff end -using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023 +using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023, MatrixDerivativeOperator using LinearAlgebra: Diagonal, LowerTriangular, diag, norm, cond using SparseArrays: spzeros +function SummationByPartsOperators.function_space_operator(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 + function inner_H1(f, g, f_derivative, g_derivative, nodes) return sum(f.(nodes) .* g.(nodes) + f_derivative.(nodes) .* g_derivative.(nodes)) end @@ -82,8 +91,8 @@ function create_P(rho, x_length) return P end -function SummationByPartsOperators.construct_function_space_operator(basis_functions, x_min, x_max, nodes, - ::GlaubitzNordströmÖffner2023) +function 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] diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index 5ab1bd25..50c3b3da 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -33,6 +33,7 @@ function Base.show(io::IO, source::GlaubitzNordströmÖffner2023) end end +# This function is extended in the package extension SummationByPartsOperatorsOptimExt """ function_space_operator(basis_functions, x_min, x_max, nodes, source; accuracy_order = 0) @@ -47,14 +48,5 @@ In order to use this function, the packages `Optim` and `ForwardDiff` must be lo See also [`GlaubitzNordströmÖffner2023`](@ref). """ -function function_space_operator(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 +function function_space_operator end -# This function is extended in the package extension SummationByPartsOperatorsOptimExt -function construct_function_space_operator end From 1b8d85e101dbdc6fce402a2947c7496b280b0c92 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 21:09:07 +0200 Subject: [PATCH 18/45] allow passing options --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 10 ++++++---- src/function_space_operators.jl | 8 ++++++-- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index cdd4f7da..ee5ddf26 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -15,9 +15,10 @@ using SparseArrays: spzeros function SummationByPartsOperators.function_space_operator(basis_functions, x_min::T, x_max::T, nodes::Vector{T}, source::SourceOfCoefficients; - accuracy_order = 0) where {T, SourceOfCoefficients} + accuracy_order = 0, + options = Options(g_tol = 1e-14, iterations = 10000)) where {T, SourceOfCoefficients} - weights, D = construct_function_space_operator(basis_functions, x_min, x_max, nodes, source) + 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 @@ -92,7 +93,8 @@ function create_P(rho, x_length) end function construct_function_space_operator(basis_functions, x_min, x_max, nodes, - ::GlaubitzNordströmÖffner2023) + ::GlaubitzNordströmÖffner2023; + options = Options(g_tol = 1e-14, iterations = 10000)) K = length(basis_functions) N = length(nodes) basis_functions_derivatives = [x -> ForwardDiff.derivative(basis_functions[i], x) for i in 1:K] @@ -125,7 +127,7 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, 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), + result = optimize(x -> optimization_function(x, p), x0, LBFGS(), options, autodiff = :forward) x = minimizer(result) sigma = x[1:L] diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index 50c3b3da..e128dd40 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -35,13 +35,17 @@ end # This function is extended in the package extension SummationByPartsOperatorsOptimExt """ - function_space_operator(basis_functions, x_min, x_max, nodes, source; accuracy_order = 0) + function_space_operator(basis_functions, x_min, x_max, nodes, source; + accuracy_order = 0, options = Optim.Options(g_tol = 1e-14, iterations = 10000)) Construct an operator that represents a 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. +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 packages `Optim` and `ForwardDiff` must be loaded. From 76f7f3f0c6eb5915c96fc07267190099a0c5d08b Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 May 2024 21:20:54 +0200 Subject: [PATCH 19/45] clarify that the operator is of derivative_order 1 --- src/function_space_operators.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index e128dd40..2d230759 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -38,7 +38,7 @@ end function_space_operator(basis_functions, x_min, x_max, nodes, source; accuracy_order = 0, options = Optim.Options(g_tol = 1e-14, iterations = 10000)) -Construct an operator that represents a derivative operator in a function space spanned by +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 From d52d55cd42302a5e7900dd520257774e4a2e587c Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 29 May 2024 13:57:43 +0200 Subject: [PATCH 20/45] only Optim has to be loaded --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 3 ++- src/function_space_operators.jl | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index ee5ddf26..e0a91bd7 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -122,13 +122,14 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, 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. + # Use the Frobenius norm since it is strictly convex and cheap to compute 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, autodiff = :forward) + display(result) x = minimizer(result) sigma = x[1:L] rho = x[(L + 1):end] diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index 2d230759..6f79e651 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -48,7 +48,7 @@ the [documentation of Optim.jl](https://julianlsolvers.github.io/Optim.jl/stable 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 packages `Optim` and `ForwardDiff` must be loaded. +In order to use this function, the package `Optim` must be loaded. See also [`GlaubitzNordströmÖffner2023`](@ref). """ From aa868bd5f65432d478b3f8925552f3a92a35adf7 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 29 May 2024 14:10:00 +0200 Subject: [PATCH 21/45] only for Julia 1.9 --- src/SummationByPartsOperators.jl | 3 -- src/function_space_operators.jl | 3 ++ test/matrix_operators_test.jl | 62 ++++++++++++++++---------------- 3 files changed, 35 insertions(+), 33 deletions(-) diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index 78e6ddaf..52a6c2eb 100644 --- a/src/SummationByPartsOperators.jl +++ b/src/SummationByPartsOperators.jl @@ -100,9 +100,6 @@ 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 diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index 6f79e651..c0a0f27c 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -51,6 +51,9 @@ The operator that is returned follows the general interface. Currently, it is wr 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 diff --git a/test/matrix_operators_test.jl b/test/matrix_operators_test.jl index 803a50d0..12321583 100644 --- a/test/matrix_operators_test.jl +++ b/test/matrix_operators_test.jl @@ -115,37 +115,39 @@ using Optim 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 = function_space_operator(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 +if VERSION >= v"1.9" + @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 = function_space_operator(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 = function_space_operator(basis_functions, x_min, x_max, nodes, source) + let basis_functions = [one, identity, exp] + D = function_space_operator(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 + @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 end From ad62dcc808e119558e54924a44f548c645639f28 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 29 May 2024 14:34:10 +0200 Subject: [PATCH 22/45] remove fallback imports with Requires --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index e0a91bd7..0bfc345e 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -1,12 +1,7 @@ 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 Optim: Options, LBFGS, optimize, minimizer +using ForwardDiff: ForwardDiff using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023, MatrixDerivativeOperator using LinearAlgebra: Diagonal, LowerTriangular, diag, norm, cond @@ -129,7 +124,7 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, x0 = zeros(L + N) result = optimize(x -> optimization_function(x, p), x0, LBFGS(), options, autodiff = :forward) - display(result) + x = minimizer(result) sigma = x[1:L] rho = x[(L + 1):end] From 2144cd15775dd940b927cf7ca7ee1f926198bd46 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Thu, 30 May 2024 11:24:54 +0200 Subject: [PATCH 23/45] use PreallocationTools.jl to reduce allocations --- Project.toml | 1 + ...mmationByPartsOperatorsOptimForwardDiffExt.jl | 16 ++++++++++++---- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 5299f448..634b0371 100644 --- a/Project.toml +++ b/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" diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 0bfc345e..8d103ae3 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -5,6 +5,7 @@ using ForwardDiff: ForwardDiff using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023, MatrixDerivativeOperator using LinearAlgebra: Diagonal, LowerTriangular, diag, norm, cond +using PreallocationTools: DiffCache, get_tmp using SparseArrays: spzeros function SummationByPartsOperators.function_space_operator(basis_functions, @@ -69,6 +70,10 @@ end function create_S(sigma, N) S = zeros(eltype(sigma), N, N) + set_S!(S, sigma, N) +end + +function set_S!(S, sigma, N) k = 1 for i in 1:N for j in (i + 1):N @@ -107,14 +112,17 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, R = B * V / 2 x_length = x_max - x_min - p = (W, R, x_length) + S = zeros(eltype(nodes), N, N) + S_cache = DiffCache(S) + p = (W, R, x_length, S_cache) @views function optimization_function(x, p) - W, R, x_length = p + W, R, x_length, S_cache = p + S = get_tmp(S_cache, x) (N, _) = size(R) L = Integer(N*(N - 1)/2) sigma = x[1:L] rho = x[(L + 1):end] - S = create_S(sigma, N) + S = set_S!(S, sigma, N) P = create_P(rho, x_length) X = [S P] # Use the Frobenius norm since it is strictly convex and cheap to compute @@ -124,7 +132,7 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, x0 = zeros(L + N) result = optimize(x -> optimization_function(x, p), x0, LBFGS(), options, autodiff = :forward) - + display(result) x = minimizer(result) sigma = x[1:L] rho = x[(L + 1):end] From 29f6830a64d4cc0ce48e1f85d35aec6c98e1a2e9 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Thu, 30 May 2024 11:29:43 +0200 Subject: [PATCH 24/45] add kwarg derivative-order --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 3 ++- src/function_space_operators.jl | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 8d103ae3..23d7bd4f 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -11,9 +11,10 @@ using SparseArrays: spzeros function SummationByPartsOperators.function_space_operator(basis_functions, x_min::T, x_max::T, nodes::Vector{T}, source::SourceOfCoefficients; - accuracy_order = 0, + 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" 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 diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index c0a0f27c..40b47f04 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -36,7 +36,8 @@ end # This function is extended in the package extension SummationByPartsOperatorsOptimExt """ function_space_operator(basis_functions, x_min, x_max, nodes, source; - accuracy_order = 0, options = Optim.Options(g_tol = 1e-14, iterations = 10000)) + 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 From d4ec79efbcff890303508b0b517e6a9765590da9 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 31 May 2024 14:00:19 +0200 Subject: [PATCH 25/45] add analytical gradient (still much slower than ForwardDiff.jl) --- ...tionByPartsOperatorsOptimForwardDiffExt.jl | 101 ++++++++++++++++-- 1 file changed, 90 insertions(+), 11 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 23d7bd4f..0c4c9ca1 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -4,7 +4,7 @@ using Optim: Options, LBFGS, optimize, minimizer using ForwardDiff: ForwardDiff using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023, MatrixDerivativeOperator -using LinearAlgebra: Diagonal, LowerTriangular, diag, norm, cond +using LinearAlgebra: Diagonal, LowerTriangular, diag, norm, mul! using PreallocationTools: DiffCache, get_tmp using SparseArrays: spzeros @@ -72,6 +72,7 @@ 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) @@ -79,13 +80,14 @@ function set_S!(S, sigma, N) 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 - return S - S' 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)) @@ -98,6 +100,7 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, options = Options(g_tol = 1e-14, iterations = 10000)) K = length(basis_functions) N = length(nodes) + L = Integer(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, @@ -105,7 +108,7 @@ function construct_function_space_operator(basis_functions, x_min, x_max, 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] + # W = [V; -V_x] B = spzeros(eltype(nodes), N, N) B[1, 1] = -1 @@ -114,25 +117,101 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, 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) - p = (W, R, x_length, S_cache) + 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) @views function optimization_function(x, p) - W, R, x_length, S_cache = p + V, V_x, R, x_length, S_cache, A_cache, SV_cache, PV_x_cache = 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 = Integer(N*(N - 1)/2) sigma = x[1:L] rho = x[(L + 1):end] - S = set_S!(S, sigma, N) + set_S!(S, sigma, N) P = create_P(rho, x_length) - X = [S P] + mul!(SV, S, V) + mul!(PV_x, P, V_x) + A .= SV .- PV_x .+ R # Use the Frobenius norm since it is strictly convex and cheap to compute - return norm(X * W + R)^2 + return norm(A)^2 end - L = Integer(N*(N - 1)/2) + @views function optimization_function_grad!(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) + 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 + fill!(daij_dsigmak, zero(eltype(daij_dsigmak))) + for i in 1:N + for j in 1:K + for k in 1:L + l_tilde = Integer(k + i - (i - 1) * (N - i/2)) + C = N^2 - 3*N + 2*i - 2*k + 1/4 + if C >= 0 + D = sqrt(C) + if isinteger(1/2 + D) + l_hat1 = Integer(N - 1/2 + D) + l_hat2 = Integer(N - 1/2 - D) + 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 + if i + 1 <= l_tilde <= N + daij_dsigmak[i, j, k] += V[l_tilde, j] + end + end + end + end + sig_rho = sig.(rho) + sig_deriv_rho = sig_deriv.(rho) + sum_sig_rho = sum(sig_rho) + for i in 1:N + for j in 1:K + factor1 = x_length * V_x[i, j] / sum_sig_rho^2 + for k in 1:N + 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 1:L + G[k] = 2 * sum(daij_dsigmak[:, :, k] .* A) + end + for k in 1:N + G[L + k] = 2 * sum(daij_drhok[:, :, k] .* A) + end + end + x0 = zeros(L + N) - result = optimize(x -> optimization_function(x, p), x0, LBFGS(), options, - autodiff = :forward) + opti(x) = optimization_function(x, p) + opti_grad!(G, x) = optimization_function_grad!(G, x, p) + result = optimize(opti, opti_grad!, x0, LBFGS(), options; inplace = true) + # result = optimize(opti, x0, LBFGS(), options; autodiff = :forward) display(result) x = minimizer(result) sigma = x[1:L] From 578a7d21e36b18c06064df19804653df04f3242a Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 31 May 2024 14:55:33 +0200 Subject: [PATCH 26/45] fuse function and gradient evaluation --- ...tionByPartsOperatorsOptimForwardDiffExt.jl | 128 ++++++++++-------- 1 file changed, 69 insertions(+), 59 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 0c4c9ca1..55929ce9 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -1,6 +1,6 @@ module SummationByPartsOperatorsOptimForwardDiffExt -using Optim: Options, LBFGS, optimize, minimizer +using Optim: Optim, Options, LBFGS, optimize, minimizer using ForwardDiff: ForwardDiff using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023, MatrixDerivativeOperator @@ -127,8 +127,26 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, 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) - @views function optimization_function(x, p) - V, V_x, R, x_length, S_cache, A_cache, SV_cache, PV_x_cache = p + # @views function optimization_function(x, p) + # V, V_x, R, x_length, S_cache, A_cache, SV_cache, PV_x_cache = 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 = Integer(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 + # # Use the Frobenius norm since it is strictly convex and cheap to compute + # return norm(A)^2 + # 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) @@ -141,76 +159,68 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, P = create_P(rho, x_length) mul!(SV, S, V) mul!(PV_x, P, V_x) - A .= SV .- PV_x .+ R - # Use the Frobenius norm since it is strictly convex and cheap to compute - return norm(A)^2 - end - @views function optimization_function_grad!(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) - 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 - fill!(daij_dsigmak, zero(eltype(daij_dsigmak))) - for i in 1:N - for j in 1:K - for k in 1:L - l_tilde = Integer(k + i - (i - 1) * (N - i/2)) - C = N^2 - 3*N + 2*i - 2*k + 1/4 - if C >= 0 - D = sqrt(C) - if isinteger(1/2 + D) - l_hat1 = Integer(N - 1/2 + D) - l_hat2 = Integer(N - 1/2 - D) - 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] + if !isnothing(G) + fill!(daij_dsigmak, zero(eltype(daij_dsigmak))) + for i in 1:N + for j in 1:K + for k in 1:L + l_tilde = Integer(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) + if isinteger(1/2 + D) + l_hat1 = Integer(N - 1/2 + D) + l_hat2 = Integer(N - 1/2 - D) + 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 - if i + 1 <= l_tilde <= N - daij_dsigmak[i, j, k] += V[l_tilde, j] - end end end - end - sig_rho = sig.(rho) - sig_deriv_rho = sig_deriv.(rho) - sum_sig_rho = sum(sig_rho) - for i in 1:N - for j in 1:K - factor1 = x_length * V_x[i, j] / sum_sig_rho^2 - for k in 1:N - 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] + sig_rho = sig.(rho) + sig_deriv_rho = sig_deriv.(rho) + sum_sig_rho = sum(sig_rho) + for i in 1:N + for j in 1:K + factor1 = x_length * V_x[i, j] / sum_sig_rho^2 + for k in 1:N + 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 1:L + G[k] = 2 * sum(daij_dsigmak[:, :, k] .* A) + end + for k in 1:N + G[L + k] = 2 * sum(daij_drhok[:, :, k] .* A) + end end - for k in 1:L - G[k] = 2 * sum(daij_dsigmak[:, :, k] .* A) - end - for k in 1:N - G[L + k] = 2 * sum(daij_drhok[:, :, k] .* A) + if !isnothing(F) + return norm(A)^2 end end x0 = zeros(L + N) - opti(x) = optimization_function(x, p) - opti_grad!(G, x) = optimization_function_grad!(G, x, p) - result = optimize(opti, opti_grad!, x0, LBFGS(), options; inplace = true) + # opti(x) = optimization_function(x, p) + # opti_grad!(G, x) = optimization_function_grad!(G, x, p) + fg!(F, G, x) = optimization_function_and_grad!(F, G, x, p) + result = optimize(Optim.only_fg!(fg!), x0, LBFGS(), options) + # result = optimize(opti, opti_grad!, x0, LBFGS(), options) # result = optimize(opti, x0, LBFGS(), options; autodiff = :forward) display(result) x = minimizer(result) From fff300b1f00ef4a71f56fa1d5b9d584730f866ca Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 31 May 2024 19:33:53 +0200 Subject: [PATCH 27/45] avoid some type conversions --- Project.toml | 1 + ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 8 +++++--- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 634b0371..34914199 100644 --- a/Project.toml +++ b/Project.toml @@ -51,6 +51,7 @@ 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" diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 55929ce9..6055a525 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -165,7 +165,9 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, for i in 1:N for j in 1:K for k in 1:L - l_tilde = Integer(k + i - (i - 1) * (N - i/2)) + ll = (i - 1)*(N - i) + sum(1:(i - 1)) # =(i - 1) * (N - i/2) + l_tilde = k + i - ll + # 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 @@ -173,8 +175,8 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, if C >= 0 D = sqrt(C) if isinteger(1/2 + D) - l_hat1 = Integer(N - 1/2 + D) - l_hat2 = Integer(N - 1/2 - D) + l_hat1 = Int(N - 1/2 + D) + l_hat2 = Int(N - 1/2 - D) if 1 <= l_hat1 <= i - 1 daij_dsigmak[i, j, k] -= V[l_hat1, j] end From 9af4c6e7bdf3ba7b00cd82a9cd95d8809e256846 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 31 May 2024 19:59:08 +0200 Subject: [PATCH 28/45] use div instead of sum --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 6055a525..9512bc1c 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -165,8 +165,8 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, for i in 1:N for j in 1:K for k in 1:L - ll = (i - 1)*(N - i) + sum(1:(i - 1)) # =(i - 1) * (N - i/2) - l_tilde = k + i - ll + 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] From 540db12df4d9b0f04ec3acd624b1d99f0d6f703b Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Fri, 31 May 2024 20:04:27 +0200 Subject: [PATCH 29/45] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 9512bc1c..363b34e3 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -100,7 +100,7 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, options = Options(g_tol = 1e-14, iterations = 10000)) K = length(basis_functions) N = length(nodes) - L = Integer(N*(N - 1)/2) + 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, @@ -152,7 +152,7 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, SV = get_tmp(SV_cache, x) PV_x = get_tmp(PV_x_cache, x) (N, _) = size(R) - L = Integer(N*(N - 1)/2) + L = div(N * (N - 1), 2) sigma = x[1:L] rho = x[(L + 1):end] set_S!(S, sigma, N) From c133295f4e1cab0c79efe2ba92c34cb232d6e7dd Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 31 May 2024 20:07:24 +0200 Subject: [PATCH 30/45] change order of loops --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 363b34e3..d66dcff3 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -162,9 +162,9 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, @. A = SV - PV_x + R if !isnothing(G) fill!(daij_dsigmak, zero(eltype(daij_dsigmak))) - for i in 1:N + for k in 1:L for j in 1:K - for k in 1:L + for i in 1:N 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)) @@ -192,10 +192,10 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, sig_rho = sig.(rho) sig_deriv_rho = sig_deriv.(rho) sum_sig_rho = sum(sig_rho) - for i in 1:N + for k in 1:N for j in 1:K - factor1 = x_length * V_x[i, j] / sum_sig_rho^2 - for k in 1:N + for i in 1:N + 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]) From f77df41aeae75f8c06bcabe17d552ae33b6648bf Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Fri, 31 May 2024 20:09:11 +0200 Subject: [PATCH 31/45] Update ext/SummationByPartsOperatorsOptimForwardDiffExt.jl Co-authored-by: Hendrik Ranocha --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index d66dcff3..2b9dd712 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -206,10 +206,10 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, end end for k in 1:L - G[k] = 2 * sum(daij_dsigmak[:, :, k] .* A) + G[k] = 2 * dot(daij_dsigmak[:, :, k], A) end for k in 1:N - G[L + k] = 2 * sum(daij_drhok[:, :, k] .* A) + G[L + k] = 2 * dot(daij_drhok[:, :, k], A) end end if !isnothing(F) From 519edcd9fc70f6c16801de7cfbd491437ca4eb8b Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 31 May 2024 20:10:15 +0200 Subject: [PATCH 32/45] import dot --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 2b9dd712..0929d7ca 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -4,7 +4,7 @@ using Optim: Optim, Options, LBFGS, optimize, minimizer using ForwardDiff: ForwardDiff using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023, MatrixDerivativeOperator -using LinearAlgebra: Diagonal, LowerTriangular, diag, norm, mul! +using LinearAlgebra: Diagonal, LowerTriangular, dot, diag, norm, mul! using PreallocationTools: DiffCache, get_tmp using SparseArrays: spzeros From feacf4e5d5c451800f4dd61d00375a9c6413ece3 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 31 May 2024 20:52:29 +0200 Subject: [PATCH 33/45] move optimization_function_and_grad to fix type instability introduced by closure bug --- ...tionByPartsOperatorsOptimForwardDiffExt.jl | 144 +++++++++--------- 1 file changed, 72 insertions(+), 72 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 0929d7ca..6aadfdb1 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -145,78 +145,6 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, # # Use the Frobenius norm since it is strictly convex and cheap to compute # return norm(A)^2 # 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 1:L - for j in 1:K - for i in 1:N - 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) - if isinteger(1/2 + D) - l_hat1 = Int(N - 1/2 + D) - l_hat2 = Int(N - 1/2 - D) - 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 1:N - for j in 1:K - for i in 1:N - 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 1:L - G[k] = 2 * dot(daij_dsigmak[:, :, k], A) - end - for k in 1:N - G[L + k] = 2 * dot(daij_drhok[:, :, k], A) - end - end - if !isnothing(F) - return norm(A)^2 - end - end - x0 = zeros(L + N) # opti(x) = optimization_function(x, p) # opti_grad!(G, x) = optimization_function_grad!(G, x, p) @@ -236,4 +164,76 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, 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) + if isinteger(1/2 + D) + l_hat1 = Int(N - 1/2 + D) + l_hat2 = Int(N - 1/2 - D) + 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 From 8ebe174e45deb48726cb7e029e83fd0784becba1 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 31 May 2024 20:53:18 +0200 Subject: [PATCH 34/45] move optimization_function_and_grad! outside (that was it!) --- ...tionByPartsOperatorsOptimForwardDiffExt.jl | 180 +++++++++--------- 1 file changed, 91 insertions(+), 89 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 0929d7ca..60c9bc64 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -95,6 +95,97 @@ function create_P(rho, x_length) return P end +# @views function optimization_function(x, p) +# V, V_x, R, x_length, S_cache, A_cache, SV_cache, PV_x_cache = 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 = Integer(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 +# # Use the Frobenius norm since it is strictly convex and cheap to compute +# return norm(A)^2 +# 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, K) = 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 1:L + for j in 1:K + for i in 1:N + 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) + if isinteger(1/2 + D) + l_hat1 = Int(N - 1/2 + D) + l_hat2 = Int(N - 1/2 - D) + 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 1:N + for j in 1:K + for i in 1:N + 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 1:L + G[k] = 2 * dot(daij_dsigmak[:, :, k], A) + end + for k in 1:N + G[L + k] = 2 * dot(daij_drhok[:, :, k], A) + end + end + if !isnothing(F) + return norm(A)^2 + end +end + function construct_function_space_operator(basis_functions, x_min, x_max, nodes, ::GlaubitzNordströmÖffner2023; options = Options(g_tol = 1e-14, iterations = 10000)) @@ -127,95 +218,6 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, 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) - # @views function optimization_function(x, p) - # V, V_x, R, x_length, S_cache, A_cache, SV_cache, PV_x_cache = 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 = Integer(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 - # # Use the Frobenius norm since it is strictly convex and cheap to compute - # return norm(A)^2 - # 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 1:L - for j in 1:K - for i in 1:N - 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) - if isinteger(1/2 + D) - l_hat1 = Int(N - 1/2 + D) - l_hat2 = Int(N - 1/2 - D) - 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 1:N - for j in 1:K - for i in 1:N - 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 1:L - G[k] = 2 * dot(daij_dsigmak[:, :, k], A) - end - for k in 1:N - G[L + k] = 2 * dot(daij_drhok[:, :, k], A) - end - end - if !isnothing(F) - return norm(A)^2 - end - end x0 = zeros(L + N) # opti(x) = optimization_function(x, p) From 047051b0c3bd1c575a3a1bb3bdb0cef343973dee Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Fri, 31 May 2024 21:09:01 +0200 Subject: [PATCH 35/45] Update ext/SummationByPartsOperatorsOptimForwardDiffExt.jl Co-authored-by: Hendrik Ranocha --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index c98ff975..0c6ea5a2 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -195,9 +195,12 @@ end C = N^2 - 3*N + 2*i - 2*k + 1/4 if C >= 0 D = sqrt(C) - if isinteger(1/2 + D) - l_hat1 = Int(N - 1/2 + D) - l_hat2 = Int(N - 1/2 - D) + 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 From 0bd5619069de330f32976c5015a38374cb62e789 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 31 May 2024 21:12:36 +0200 Subject: [PATCH 36/45] some cleanups --- ...tionByPartsOperatorsOptimForwardDiffExt.jl | 23 ------------------- 1 file changed, 23 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 0c6ea5a2..d22befde 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -129,12 +129,8 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, p = (V, V_x, R, x_length, S_cache, A_cache, SV_cache, PV_x_cache, daij_dsigmak, daij_drhok) x0 = zeros(L + N) - # opti(x) = optimization_function(x, p) - # opti_grad!(G, x) = optimization_function_grad!(G, x, p) fg!(F, G, x) = optimization_function_and_grad!(F, G, x, p) result = optimize(Optim.only_fg!(fg!), x0, LBFGS(), options) - # result = optimize(opti, opti_grad!, x0, LBFGS(), options) - # result = optimize(opti, x0, LBFGS(), options; autodiff = :forward) display(result) x = minimizer(result) sigma = x[1:L] @@ -147,25 +143,6 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, return weights, D end -# @views function optimization_function(x, p) -# V, V_x, R, x_length, S_cache, A_cache, SV_cache, PV_x_cache = 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 = Integer(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 -# # Use the Frobenius norm since it is strictly convex and cheap to compute -# return norm(A)^2 -# 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) From 914cab21119aba27754e15279058bb50d4b705c2 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Thu, 13 Jun 2024 16:32:53 +0200 Subject: [PATCH 37/45] ignore stale dependency PreallocationTools in Aqua --- test/aqua.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/aqua.jl b/test/aqua.jl index ec31e437..e4487da8 100644 --- a/test/aqua.jl +++ b/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]), ) From 1a2a37605f61f92090a44ac6ba90147840f4d368 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 14 Jun 2024 12:45:05 +0200 Subject: [PATCH 38/45] add test with non-equidistant nodes --- ...tionByPartsOperatorsOptimForwardDiffExt.jl | 3 +- test/matrix_operators_test.jl | 34 +++++++++++++++---- 2 files changed, 29 insertions(+), 8 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index d22befde..dc5cfb8b 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -15,6 +15,7 @@ function SummationByPartsOperators.function_space_operator(basis_functions, 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 @@ -131,7 +132,7 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, 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) - display(result) + x = minimizer(result) sigma = x[1:L] rho = x[(L + 1):end] diff --git a/test/matrix_operators_test.jl b/test/matrix_operators_test.jl index 12321583..c738e080 100644 --- a/test/matrix_operators_test.jl +++ b/test/matrix_operators_test.jl @@ -118,10 +118,9 @@ end if VERSION >= v"1.9" @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] + nodes = collect(range(x_min, x_max, length=N)) source = GlaubitzNordströmÖffner2023() for compact in (true, false) show(IOContext(devnull, :compact=>compact), source) @@ -132,10 +131,11 @@ if VERSION >= v"1.9" let basis_functions = [x -> x^i for i in 0:3] D = function_space_operator(basis_functions, x_min, x_max, nodes, source) + @test grid(D) ≈ nodes @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) + @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 @@ -143,9 +143,29 @@ if VERSION >= v"1.9" let basis_functions = [one, identity, exp] D = function_space_operator(basis_functions, x_min, x_max, nodes, source) + @test grid(D) ≈ nodes @test ≈(D * ones(N), zeros(N); atol = 1e-13) - @test D * x ≈ ones(N) - @test D * exp.(x) ≈ exp.(x) + @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, first(nodes), last(nodes), nodes, source) + + @test grid(D) ≈ nodes + @test ≈(D * ones(N), zeros(N); atol = 5e-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 From 03107f9407c00343c8b6d21e9af53069a0310f24 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 14 Jun 2024 14:56:58 +0200 Subject: [PATCH 39/45] fix test --- test/matrix_operators_test.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/matrix_operators_test.jl b/test/matrix_operators_test.jl index c738e080..212dce3b 100644 --- a/test/matrix_operators_test.jl +++ b/test/matrix_operators_test.jl @@ -132,7 +132,7 @@ if VERSION >= v"1.9" D = function_space_operator(basis_functions, x_min, x_max, nodes, source) @test grid(D) ≈ nodes - @test ≈(D * ones(N), zeros(N); atol = 1e-13) + @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) @@ -144,7 +144,7 @@ if VERSION >= v"1.9" D = function_space_operator(basis_functions, x_min, x_max, nodes, source) @test grid(D) ≈ nodes - @test ≈(D * ones(N), zeros(N); atol = 1e-13) + @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) @@ -163,7 +163,7 @@ if VERSION >= v"1.9" D = function_space_operator(basis_functions, first(nodes), last(nodes), nodes, source) @test grid(D) ≈ nodes - @test ≈(D * ones(N), zeros(N); atol = 5e-13) + @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) From 3e03c7370c898da0f4bc20e4ab87dc4d77555b7a Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 14 Jun 2024 15:38:48 +0200 Subject: [PATCH 40/45] remove kwargs x_min and x_max --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 11 +++++------ src/function_space_operators.jl | 13 +++++++------ test/matrix_operators_test.jl | 6 +++--- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index dc5cfb8b..c20b7ddb 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -8,16 +8,15 @@ 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}, +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} @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) + 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) @@ -96,7 +95,7 @@ function create_P(rho, x_length) return P end -function construct_function_space_operator(basis_functions, x_min, x_max, nodes, +function construct_function_space_operator(basis_functions, nodes, ::GlaubitzNordströmÖffner2023; options = Options(g_tol = 1e-14, iterations = 10000)) K = length(basis_functions) @@ -116,7 +115,7 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes, B[N, N] = 1 R = B * V / 2 - x_length = x_max - x_min + x_length = last(nodes) - first(nodes) S = zeros(eltype(nodes), N, N) A = zeros(eltype(nodes), N, K) SV = zeros(eltype(nodes), N, K) diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index 40b47f04..b255ee66 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -35,17 +35,18 @@ end # This function is extended in the package extension SummationByPartsOperatorsOptimExt """ - function_space_operator(basis_functions, x_min, x_max, nodes, source; + 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`. 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/). +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. diff --git a/test/matrix_operators_test.jl b/test/matrix_operators_test.jl index 212dce3b..15761959 100644 --- a/test/matrix_operators_test.jl +++ b/test/matrix_operators_test.jl @@ -129,7 +129,7 @@ if VERSION >= v"1.9" 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, x_min, x_max, nodes, source) + D = function_space_operator(basis_functions, nodes, source) @test grid(D) ≈ nodes @test all(isapprox.(D * ones(N), zeros(N); atol = 1e-13)) @@ -141,7 +141,7 @@ if VERSION >= v"1.9" end let basis_functions = [one, identity, exp] - D = function_space_operator(basis_functions, x_min, x_max, nodes, source) + D = function_space_operator(basis_functions, nodes, source) @test grid(D) ≈ nodes @test all(isapprox.(D * ones(N), zeros(N); atol = 1e-13)) @@ -160,7 +160,7 @@ if VERSION >= v"1.9" B[1, 1] = -1.0 B[N, N] = 1.0 let basis_functions = [one, identity, exp] - D = function_space_operator(basis_functions, first(nodes), last(nodes), nodes, source) + D = function_space_operator(basis_functions, nodes, source) @test grid(D) ≈ nodes @test all(isapprox.(D * ones(N), zeros(N); atol = 1e-11)) From 38081b1822de3efbadef4a23afe5e6dcacd87155 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 14 Jun 2024 15:41:24 +0200 Subject: [PATCH 41/45] add experimental feature warning --- src/function_space_operators.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/function_space_operators.jl b/src/function_space_operators.jl index b255ee66..081cd167 100644 --- a/src/function_space_operators.jl +++ b/src/function_space_operators.jl @@ -56,6 +56,8 @@ 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 - From 4e9512c099c0f8b0484034e293137aae1d1453ef Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Fri, 14 Jun 2024 15:42:42 +0200 Subject: [PATCH 42/45] Update test/matrix_operators_test.jl Co-authored-by: Hendrik Ranocha --- test/matrix_operators_test.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/matrix_operators_test.jl b/test/matrix_operators_test.jl index 15761959..1279de58 100644 --- a/test/matrix_operators_test.jl +++ b/test/matrix_operators_test.jl @@ -1,7 +1,7 @@ using Test using LinearAlgebra using SummationByPartsOperators -using Optim +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 From fd470ee623f3a9585a5be3c6558fd2e99cf9af98 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 14 Jun 2024 16:00:08 +0200 Subject: [PATCH 43/45] remove l_hat2 since it this case is impossible --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index c20b7ddb..920ad83d 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -176,13 +176,9 @@ end 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] + 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 From 7f54a303cfcb0c0f5be49700223c546de6ce83e3 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 14 Jun 2024 16:12:50 +0200 Subject: [PATCH 44/45] remove PreallocationTools again --- Project.toml | 2 -- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 13 ++----------- test/aqua.jl | 2 +- 3 files changed, 3 insertions(+), 14 deletions(-) diff --git a/Project.toml b/Project.toml index 34914199..5299f448 100644 --- a/Project.toml +++ b/Project.toml @@ -12,7 +12,6 @@ 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" @@ -51,7 +50,6 @@ 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" diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 920ad83d..8aa4be77 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -5,7 +5,6 @@ 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, nodes::Vector{T}, @@ -120,13 +119,9 @@ function construct_function_space_operator(basis_functions, nodes, 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) + 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) @@ -144,11 +139,7 @@ function construct_function_space_operator(basis_functions, nodes, 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) + 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] diff --git a/test/aqua.jl b/test/aqua.jl index e4487da8..ec31e437 100644 --- a/test/aqua.jl +++ b/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, :PreallocationTools]), + stale_deps = (; ignore = [:PrecompileTools]), ) From b62f40bfc103fda067b0618218acbad210bdd11a Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 14 Jun 2024 16:48:36 +0200 Subject: [PATCH 45/45] throw ArgumentError instead of assertion for wrong derivative_order --- ext/SummationByPartsOperatorsOptimForwardDiffExt.jl | 4 +++- test/matrix_operators_test.jl | 2 ++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl index 8aa4be77..b4801d41 100644 --- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl +++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl @@ -12,7 +12,9 @@ function SummationByPartsOperators.function_space_operator(basis_functions, node 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" + 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) diff --git a/test/matrix_operators_test.jl b/test/matrix_operators_test.jl index 1279de58..f7a31092 100644 --- a/test/matrix_operators_test.jl +++ b/test/matrix_operators_test.jl @@ -130,6 +130,8 @@ if VERSION >= v"1.9" 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))