-
Notifications
You must be signed in to change notification settings - Fork 12
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
* first draft for function space operators * some spaces * order alphabetically * don't force basis_function to be passed as vector * allow passing accuracy_order * first (basically) working version * orthonormalization via Gram-Schmidt * add some basic tests * some smaller optimizations * fix import of Options * another fix... * avoid unnecessary anonymous function * use version of spzeros that is compatible with older Julia versions * FunctionSpaceOperator -> function_space_operator * add proper docstring * extend docstring * put constructor function in extension * allow passing options * clarify that the operator is of derivative_order 1 * only Optim has to be loaded * only for Julia 1.9 * remove fallback imports with Requires * use PreallocationTools.jl to reduce allocations * add kwarg derivative-order * add analytical gradient (still much slower than ForwardDiff.jl) * fuse function and gradient evaluation * avoid some type conversions * use div instead of sum * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * change order of loops * Update ext/SummationByPartsOperatorsOptimForwardDiffExt.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * import dot * move optimization_function_and_grad to fix type instability introduced by closure bug * move optimization_function_and_grad! outside (that was it!) * Update ext/SummationByPartsOperatorsOptimForwardDiffExt.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * some cleanups * ignore stale dependency PreallocationTools in Aqua * add test with non-equidistant nodes * fix test * remove kwargs x_min and x_max * add experimental feature warning * Update test/matrix_operators_test.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * remove l_hat2 since it this case is impossible * remove PreallocationTools again * throw ArgumentError instead of assertion for wrong derivative_order --------- Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> Co-authored-by: Hendrik Ranocha <mail@ranocha.de>
- Loading branch information
3 people
committed
Jun 14, 2024
1 parent
96e695d
commit 0c50276
Showing
7 changed files
with
345 additions
and
1 deletion.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -9,3 +9,6 @@ docs/build | |
| docs/src/code_of_conduct.md | ||
| docs/src/contributing.md | ||
| docs/src/license.md | ||
|
|
||
| run | ||
| run/* | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,210 @@ | ||
| module SummationByPartsOperatorsOptimForwardDiffExt | ||
|
|
||
| using Optim: Optim, Options, LBFGS, optimize, minimizer | ||
| using ForwardDiff: ForwardDiff | ||
|
|
||
| using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023, MatrixDerivativeOperator | ||
| using LinearAlgebra: Diagonal, LowerTriangular, dot, diag, norm, mul! | ||
| using SparseArrays: spzeros | ||
|
|
||
| function SummationByPartsOperators.function_space_operator(basis_functions, nodes::Vector{T}, | ||
| source::SourceOfCoefficients; | ||
| derivative_order = 1, accuracy_order = 0, | ||
| options = Options(g_tol = 1e-14, iterations = 10000)) where {T, SourceOfCoefficients} | ||
|
|
||
| if derivative_order != 1 | ||
| throw(ArgumentError("Derivative order $derivative_order not implemented.")) | ||
| end | ||
| sort!(nodes) | ||
| weights, D = construct_function_space_operator(basis_functions, nodes, source; options = options) | ||
| return MatrixDerivativeOperator(first(nodes), last(nodes), nodes, weights, D, accuracy_order, source) | ||
| end | ||
|
|
||
| function inner_H1(f, g, f_derivative, g_derivative, nodes) | ||
| return sum(f.(nodes) .* g.(nodes) + f_derivative.(nodes) .* g_derivative.(nodes)) | ||
| end | ||
| norm_H1(f, f_derivative, nodes) = sqrt(inner_H1(f, f, f_derivative, f_derivative, nodes)) | ||
|
|
||
| call_orthonormal_basis_function(A, basis_functions, k, x) = sum([basis_functions[i](x)*A[k, i] for i in 1:k]) | ||
|
|
||
| # This will orthonormalize the basis functions using the Gram-Schmidt process to reduce the condition | ||
| # number of the Vandermonde matrix. The matrix A transfers the old basis functions to the new orthonormalized by | ||
| # g(x) = A * f(x), where f(x) is the vector of old basis functions and g(x) is the vector of the new orthonormalized | ||
| # basis functions. Analogously, we have g'(x) = A * f'(x). | ||
| function orthonormalize_gram_schmidt(basis_functions, basis_functions_derivatives, nodes) | ||
| K = length(basis_functions) | ||
|
|
||
| A = LowerTriangular(zeros(eltype(nodes), K, K)) | ||
|
|
||
| basis_functions_orthonormalized = Vector{Function}(undef, K) | ||
| basis_functions_orthonormalized_derivatives = Vector{Function}(undef, K) | ||
|
|
||
| for k = 1:K | ||
| A[k, k] = 1 | ||
| for j = 1:k-1 | ||
| g(x) = call_orthonormal_basis_function(A, basis_functions, j, x) | ||
| g_derivative(x) = call_orthonormal_basis_function(A, basis_functions_derivatives, j, x) | ||
| inner_product = inner_H1(basis_functions[k], g, basis_functions_derivatives[k], g_derivative, nodes) | ||
| norm_squared = inner_H1(g, g, g_derivative, g_derivative, nodes) | ||
| A[k, :] = A[k, :] - inner_product/norm_squared * A[j, :] | ||
| end | ||
|
|
||
| basis_functions_orthonormalized[k] = x -> call_orthonormal_basis_function(A, basis_functions, k, x) | ||
| basis_functions_orthonormalized_derivatives[k] = x -> call_orthonormal_basis_function(A, basis_functions_derivatives, k, x) | ||
| # Normalization | ||
| r = norm_H1(basis_functions_orthonormalized[k], basis_functions_orthonormalized_derivatives[k], nodes) | ||
| A[k, :] = A[k, :] / r | ||
| end | ||
| return basis_functions_orthonormalized, basis_functions_orthonormalized_derivatives | ||
| end | ||
|
|
||
| function vandermonde_matrix(functions, nodes) | ||
| N = length(nodes) | ||
| K = length(functions) | ||
| V = zeros(eltype(nodes), N, K) | ||
| for i in 1:N | ||
| for j in 1:K | ||
| V[i, j] = functions[j](nodes[i]) | ||
| end | ||
| end | ||
| return V | ||
| end | ||
|
|
||
| function create_S(sigma, N) | ||
| S = zeros(eltype(sigma), N, N) | ||
| set_S!(S, sigma, N) | ||
| return S | ||
| end | ||
|
|
||
| function set_S!(S, sigma, N) | ||
| k = 1 | ||
| for i in 1:N | ||
| for j in (i + 1):N | ||
| S[i, j] = sigma[k] | ||
| S[j, i] = -sigma[k] | ||
| k += 1 | ||
| end | ||
| end | ||
| end | ||
|
|
||
| sig(x) = 1 / (1 + exp(-x)) | ||
| sig_deriv(x) = sig(x) * (1 - sig(x)) | ||
|
|
||
| function create_P(rho, x_length) | ||
| P = Diagonal(sig.(rho)) | ||
| P *= x_length / sum(P) | ||
| return P | ||
| end | ||
|
|
||
| function construct_function_space_operator(basis_functions, nodes, | ||
| ::GlaubitzNordströmÖffner2023; | ||
| options = Options(g_tol = 1e-14, iterations = 10000)) | ||
| K = length(basis_functions) | ||
| N = length(nodes) | ||
| L = div(N * (N - 1), 2) | ||
| basis_functions_derivatives = [x -> ForwardDiff.derivative(basis_functions[i], x) for i in 1:K] | ||
| basis_functions_orthonormalized, basis_functions_orthonormalized_derivatives = orthonormalize_gram_schmidt(basis_functions, | ||
| basis_functions_derivatives, | ||
| nodes) | ||
| V = vandermonde_matrix(basis_functions_orthonormalized, nodes) | ||
| V_x = vandermonde_matrix(basis_functions_orthonormalized_derivatives, nodes) | ||
| # Here, W satisfies W'*W = I | ||
| # W = [V; -V_x] | ||
|
|
||
| B = spzeros(eltype(nodes), N, N) | ||
| B[1, 1] = -1 | ||
| B[N, N] = 1 | ||
|
|
||
| R = B * V / 2 | ||
| x_length = last(nodes) - first(nodes) | ||
| S = zeros(eltype(nodes), N, N) | ||
| A = zeros(eltype(nodes), N, K) | ||
| SV = zeros(eltype(nodes), N, K) | ||
| PV_x = zeros(eltype(nodes), N, K) | ||
| daij_dsigmak = zeros(eltype(nodes), N, K, L) | ||
| daij_drhok = zeros(eltype(nodes), N, K, N) | ||
| p = (V, V_x, R, x_length, S, A, SV, PV_x, daij_dsigmak, daij_drhok) | ||
|
|
||
| x0 = zeros(L + N) | ||
| fg!(F, G, x) = optimization_function_and_grad!(F, G, x, p) | ||
| result = optimize(Optim.only_fg!(fg!), x0, LBFGS(), options) | ||
|
|
||
| x = minimizer(result) | ||
| sigma = x[1:L] | ||
| rho = x[(L + 1):end] | ||
| S = create_S(sigma, N) | ||
| P = create_P(rho, x_length) | ||
| weights = diag(P) | ||
| Q = S + B/2 | ||
| D = inv(P) * Q | ||
| return weights, D | ||
| end | ||
|
|
||
| @views function optimization_function_and_grad!(F, G, x, p) | ||
| V, V_x, R, x_length, S, A, SV, PV_x, daij_dsigmak, daij_drhok = p | ||
| (N, _) = size(R) | ||
| L = div(N * (N - 1), 2) | ||
| sigma = x[1:L] | ||
| rho = x[(L + 1):end] | ||
| set_S!(S, sigma, N) | ||
| P = create_P(rho, x_length) | ||
| mul!(SV, S, V) | ||
| mul!(PV_x, P, V_x) | ||
| @. A = SV - PV_x + R | ||
| if !isnothing(G) | ||
| fill!(daij_dsigmak, zero(eltype(daij_dsigmak))) | ||
| for k in axes(daij_dsigmak, 3) | ||
| for j in axes(daij_dsigmak, 2) | ||
| for i in axes(daij_dsigmak, 1) | ||
| l_tilde = k + i - N * (i - 1) + div(i * (i - 1), 2) | ||
| # same as above, but needs more type conversions | ||
| # l_tilde = Int(k + i - (i - 1) * (N - i/2)) | ||
| if i + 1 <= l_tilde <= N | ||
| daij_dsigmak[i, j, k] += V[l_tilde, j] | ||
| else | ||
| C = N^2 - 3*N + 2*i - 2*k + 1/4 | ||
| if C >= 0 | ||
| D = sqrt(C) | ||
| D_plus_one_half = D + 0.5 | ||
| D_plus_one_half_trunc = trunc(D_plus_one_half) | ||
| if D_plus_one_half == D_plus_one_half_trunc | ||
| int_D_plus_one_half = trunc(Int, D_plus_one_half_trunc) | ||
| l_hat = N - int_D_plus_one_half | ||
| if 1 <= l_hat <= i - 1 | ||
| daij_dsigmak[i, j, k] -= V[l_hat, j] | ||
| end | ||
| end | ||
| end | ||
| end | ||
| end | ||
| end | ||
| end | ||
| sig_rho = sig.(rho) | ||
| sig_deriv_rho = sig_deriv.(rho) | ||
| sum_sig_rho = sum(sig_rho) | ||
| for k in axes(daij_drhok, 3) | ||
| for j in axes(daij_drhok, 2) | ||
| for i in axes(daij_drhok, 1) | ||
| factor1 = x_length * V_x[i, j] / sum_sig_rho^2 | ||
| factor = factor1 * sig_deriv_rho[k] | ||
| if k == i | ||
| daij_drhok[i, j, k] = -factor * (sum_sig_rho - sig_rho[k]) | ||
| else | ||
| daij_drhok[i, j, k] = factor * sig_rho[i] | ||
| end | ||
| end | ||
| end | ||
| end | ||
| for k in axes(daij_dsigmak, 3) | ||
| G[k] = 2 * dot(daij_dsigmak[:, :, k], A) | ||
| end | ||
| for k in axes(daij_drhok, 3) | ||
| G[L + k] = 2 * dot(daij_drhok[:, :, k], A) | ||
| end | ||
| end | ||
| if !isnothing(F) | ||
| return norm(A)^2 | ||
| end | ||
| end | ||
|
|
||
| end # module |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,63 @@ | ||
|
|
||
| """ | ||
| GlaubitzNordströmÖffner2023() | ||
| Function space SBP (FSBP) operators given in | ||
| - Glaubitz, Nordström, Öffner (2023) | ||
| Summation-by-parts operators for general function spaces. | ||
| SIAM Journal on Numerical Analysis 61, 2, pp. 733-754. | ||
| See also | ||
| - Glaubitz, Nordström, Öffner (2024) | ||
| An optimization-based construction procedure for function space based | ||
| summation-by-parts operators on arbitrary grids. | ||
| arXiv, arXiv:2405.08770v1. | ||
| See [`function_space_operator`](@ref). | ||
| """ | ||
| struct GlaubitzNordströmÖffner2023 <: SourceOfCoefficients end | ||
|
|
||
| function Base.show(io::IO, source::GlaubitzNordströmÖffner2023) | ||
| if get(io, :compact, false) | ||
| summary(io, source) | ||
| else | ||
| print(io, | ||
| "Glaubitz, Nordström, Öffner (2023) \n", | ||
| " Summation-by-parts operators for general function spaces \n", | ||
| " SIAM Journal on Numerical Analysis 61, 2, pp. 733-754. \n", | ||
| "See also \n", | ||
| " Glaubitz, Nordström, Öffner (2024) \n", | ||
| " An optimization-based construction procedure for function \n", | ||
| " space based summation-by-parts operators on arbitrary grids \n", | ||
| " arXiv, arXiv:2405.08770v1.") | ||
| end | ||
| end | ||
|
|
||
| # This function is extended in the package extension SummationByPartsOperatorsOptimExt | ||
| """ | ||
| function_space_operator(basis_functions, nodes, source; | ||
| derivative_order = 1, accuracy_order = 0, | ||
| options = Optim.Options(g_tol = 1e-14, iterations = 10000)) | ||
| Construct an operator that represents a first-derivative operator in a function space spanned by | ||
| the `basis_functions`, which is an iterable of functions. The operator is constructed on the | ||
| interval `[x_min, x_max]` with the nodes `nodes`, where `x_min` is taken as the minimal value in | ||
| `nodes` and `x_max` the maximal value. Note that the `nodes` will be sorted internally. The | ||
| `accuracy_order` is the order of the accuracy of the operator, which can optionally be passed, | ||
| but does not have any effect on the operator. The operator is constructed solving an optimization | ||
| problem with Optim.jl. You can specify the options for the optimization problem with the `options` | ||
| argument, see also the [documentation of Optim.jl](https://julianlsolvers.github.io/Optim.jl/stable/user/config/). | ||
| The operator that is returned follows the general interface. Currently, it is wrapped in a | ||
| [`MatrixDerivativeOperator`](@ref), but this might change in the future. | ||
| In order to use this function, the package `Optim` must be loaded. | ||
| See also [`GlaubitzNordströmÖffner2023`](@ref). | ||
| !!! compat "Julia 1.9" | ||
| This function requires at least Julia 1.9. | ||
| !!! warning "Experimental implementation" | ||
| This is an experimental feature and may change in future releases. | ||
| """ | ||
| function function_space_operator end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.