|
| 1 | +module SummationByPartsOperatorsOptimForwardDiffExt |
| 2 | + |
| 3 | +using Optim: Optim, Options, LBFGS, optimize, minimizer |
| 4 | +using ForwardDiff: ForwardDiff |
| 5 | + |
| 6 | +using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023, MatrixDerivativeOperator |
| 7 | +using LinearAlgebra: Diagonal, LowerTriangular, dot, diag, norm, mul! |
| 8 | +using SparseArrays: spzeros |
| 9 | + |
| 10 | +function SummationByPartsOperators.function_space_operator(basis_functions, nodes::Vector{T}, |
| 11 | + source::SourceOfCoefficients; |
| 12 | + derivative_order = 1, accuracy_order = 0, |
| 13 | + options = Options(g_tol = 1e-14, iterations = 10000)) where {T, SourceOfCoefficients} |
| 14 | + |
| 15 | + if derivative_order != 1 |
| 16 | + throw(ArgumentError("Derivative order $derivative_order not implemented.")) |
| 17 | + end |
| 18 | + sort!(nodes) |
| 19 | + weights, D = construct_function_space_operator(basis_functions, nodes, source; options = options) |
| 20 | + return MatrixDerivativeOperator(first(nodes), last(nodes), nodes, weights, D, accuracy_order, source) |
| 21 | +end |
| 22 | + |
| 23 | +function inner_H1(f, g, f_derivative, g_derivative, nodes) |
| 24 | + return sum(f.(nodes) .* g.(nodes) + f_derivative.(nodes) .* g_derivative.(nodes)) |
| 25 | +end |
| 26 | +norm_H1(f, f_derivative, nodes) = sqrt(inner_H1(f, f, f_derivative, f_derivative, nodes)) |
| 27 | + |
| 28 | +call_orthonormal_basis_function(A, basis_functions, k, x) = sum([basis_functions[i](x)*A[k, i] for i in 1:k]) |
| 29 | + |
| 30 | +# This will orthonormalize the basis functions using the Gram-Schmidt process to reduce the condition |
| 31 | +# number of the Vandermonde matrix. The matrix A transfers the old basis functions to the new orthonormalized by |
| 32 | +# 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 |
| 33 | +# basis functions. Analogously, we have g'(x) = A * f'(x). |
| 34 | +function orthonormalize_gram_schmidt(basis_functions, basis_functions_derivatives, nodes) |
| 35 | + K = length(basis_functions) |
| 36 | + |
| 37 | + A = LowerTriangular(zeros(eltype(nodes), K, K)) |
| 38 | + |
| 39 | + basis_functions_orthonormalized = Vector{Function}(undef, K) |
| 40 | + basis_functions_orthonormalized_derivatives = Vector{Function}(undef, K) |
| 41 | + |
| 42 | + for k = 1:K |
| 43 | + A[k, k] = 1 |
| 44 | + for j = 1:k-1 |
| 45 | + g(x) = call_orthonormal_basis_function(A, basis_functions, j, x) |
| 46 | + g_derivative(x) = call_orthonormal_basis_function(A, basis_functions_derivatives, j, x) |
| 47 | + inner_product = inner_H1(basis_functions[k], g, basis_functions_derivatives[k], g_derivative, nodes) |
| 48 | + norm_squared = inner_H1(g, g, g_derivative, g_derivative, nodes) |
| 49 | + A[k, :] = A[k, :] - inner_product/norm_squared * A[j, :] |
| 50 | + end |
| 51 | + |
| 52 | + basis_functions_orthonormalized[k] = x -> call_orthonormal_basis_function(A, basis_functions, k, x) |
| 53 | + basis_functions_orthonormalized_derivatives[k] = x -> call_orthonormal_basis_function(A, basis_functions_derivatives, k, x) |
| 54 | + # Normalization |
| 55 | + r = norm_H1(basis_functions_orthonormalized[k], basis_functions_orthonormalized_derivatives[k], nodes) |
| 56 | + A[k, :] = A[k, :] / r |
| 57 | + end |
| 58 | + return basis_functions_orthonormalized, basis_functions_orthonormalized_derivatives |
| 59 | +end |
| 60 | + |
| 61 | +function vandermonde_matrix(functions, nodes) |
| 62 | + N = length(nodes) |
| 63 | + K = length(functions) |
| 64 | + V = zeros(eltype(nodes), N, K) |
| 65 | + for i in 1:N |
| 66 | + for j in 1:K |
| 67 | + V[i, j] = functions[j](nodes[i]) |
| 68 | + end |
| 69 | + end |
| 70 | + return V |
| 71 | +end |
| 72 | + |
| 73 | +function create_S(sigma, N) |
| 74 | + S = zeros(eltype(sigma), N, N) |
| 75 | + set_S!(S, sigma, N) |
| 76 | + return S |
| 77 | +end |
| 78 | + |
| 79 | +function set_S!(S, sigma, N) |
| 80 | + k = 1 |
| 81 | + for i in 1:N |
| 82 | + for j in (i + 1):N |
| 83 | + S[i, j] = sigma[k] |
| 84 | + S[j, i] = -sigma[k] |
| 85 | + k += 1 |
| 86 | + end |
| 87 | + end |
| 88 | +end |
| 89 | + |
| 90 | +sig(x) = 1 / (1 + exp(-x)) |
| 91 | +sig_deriv(x) = sig(x) * (1 - sig(x)) |
| 92 | + |
| 93 | +function create_P(rho, x_length) |
| 94 | + P = Diagonal(sig.(rho)) |
| 95 | + P *= x_length / sum(P) |
| 96 | + return P |
| 97 | +end |
| 98 | + |
| 99 | +function construct_function_space_operator(basis_functions, nodes, |
| 100 | + ::GlaubitzNordströmÖffner2023; |
| 101 | + options = Options(g_tol = 1e-14, iterations = 10000)) |
| 102 | + K = length(basis_functions) |
| 103 | + N = length(nodes) |
| 104 | + L = div(N * (N - 1), 2) |
| 105 | + basis_functions_derivatives = [x -> ForwardDiff.derivative(basis_functions[i], x) for i in 1:K] |
| 106 | + basis_functions_orthonormalized, basis_functions_orthonormalized_derivatives = orthonormalize_gram_schmidt(basis_functions, |
| 107 | + basis_functions_derivatives, |
| 108 | + nodes) |
| 109 | + V = vandermonde_matrix(basis_functions_orthonormalized, nodes) |
| 110 | + V_x = vandermonde_matrix(basis_functions_orthonormalized_derivatives, nodes) |
| 111 | + # Here, W satisfies W'*W = I |
| 112 | + # W = [V; -V_x] |
| 113 | + |
| 114 | + B = spzeros(eltype(nodes), N, N) |
| 115 | + B[1, 1] = -1 |
| 116 | + B[N, N] = 1 |
| 117 | + |
| 118 | + R = B * V / 2 |
| 119 | + x_length = last(nodes) - first(nodes) |
| 120 | + S = zeros(eltype(nodes), N, N) |
| 121 | + A = zeros(eltype(nodes), N, K) |
| 122 | + SV = zeros(eltype(nodes), N, K) |
| 123 | + PV_x = zeros(eltype(nodes), N, K) |
| 124 | + daij_dsigmak = zeros(eltype(nodes), N, K, L) |
| 125 | + daij_drhok = zeros(eltype(nodes), N, K, N) |
| 126 | + p = (V, V_x, R, x_length, S, A, SV, PV_x, daij_dsigmak, daij_drhok) |
| 127 | + |
| 128 | + x0 = zeros(L + N) |
| 129 | + fg!(F, G, x) = optimization_function_and_grad!(F, G, x, p) |
| 130 | + result = optimize(Optim.only_fg!(fg!), x0, LBFGS(), options) |
| 131 | + |
| 132 | + x = minimizer(result) |
| 133 | + sigma = x[1:L] |
| 134 | + rho = x[(L + 1):end] |
| 135 | + S = create_S(sigma, N) |
| 136 | + P = create_P(rho, x_length) |
| 137 | + weights = diag(P) |
| 138 | + Q = S + B/2 |
| 139 | + D = inv(P) * Q |
| 140 | + return weights, D |
| 141 | +end |
| 142 | + |
| 143 | +@views function optimization_function_and_grad!(F, G, x, p) |
| 144 | + V, V_x, R, x_length, S, A, SV, PV_x, daij_dsigmak, daij_drhok = p |
| 145 | + (N, _) = size(R) |
| 146 | + L = div(N * (N - 1), 2) |
| 147 | + sigma = x[1:L] |
| 148 | + rho = x[(L + 1):end] |
| 149 | + set_S!(S, sigma, N) |
| 150 | + P = create_P(rho, x_length) |
| 151 | + mul!(SV, S, V) |
| 152 | + mul!(PV_x, P, V_x) |
| 153 | + @. A = SV - PV_x + R |
| 154 | + if !isnothing(G) |
| 155 | + fill!(daij_dsigmak, zero(eltype(daij_dsigmak))) |
| 156 | + for k in axes(daij_dsigmak, 3) |
| 157 | + for j in axes(daij_dsigmak, 2) |
| 158 | + for i in axes(daij_dsigmak, 1) |
| 159 | + l_tilde = k + i - N * (i - 1) + div(i * (i - 1), 2) |
| 160 | + # same as above, but needs more type conversions |
| 161 | + # l_tilde = Int(k + i - (i - 1) * (N - i/2)) |
| 162 | + if i + 1 <= l_tilde <= N |
| 163 | + daij_dsigmak[i, j, k] += V[l_tilde, j] |
| 164 | + else |
| 165 | + C = N^2 - 3*N + 2*i - 2*k + 1/4 |
| 166 | + if C >= 0 |
| 167 | + D = sqrt(C) |
| 168 | + D_plus_one_half = D + 0.5 |
| 169 | + D_plus_one_half_trunc = trunc(D_plus_one_half) |
| 170 | + if D_plus_one_half == D_plus_one_half_trunc |
| 171 | + int_D_plus_one_half = trunc(Int, D_plus_one_half_trunc) |
| 172 | + l_hat = N - int_D_plus_one_half |
| 173 | + if 1 <= l_hat <= i - 1 |
| 174 | + daij_dsigmak[i, j, k] -= V[l_hat, j] |
| 175 | + end |
| 176 | + end |
| 177 | + end |
| 178 | + end |
| 179 | + end |
| 180 | + end |
| 181 | + end |
| 182 | + sig_rho = sig.(rho) |
| 183 | + sig_deriv_rho = sig_deriv.(rho) |
| 184 | + sum_sig_rho = sum(sig_rho) |
| 185 | + for k in axes(daij_drhok, 3) |
| 186 | + for j in axes(daij_drhok, 2) |
| 187 | + for i in axes(daij_drhok, 1) |
| 188 | + factor1 = x_length * V_x[i, j] / sum_sig_rho^2 |
| 189 | + factor = factor1 * sig_deriv_rho[k] |
| 190 | + if k == i |
| 191 | + daij_drhok[i, j, k] = -factor * (sum_sig_rho - sig_rho[k]) |
| 192 | + else |
| 193 | + daij_drhok[i, j, k] = factor * sig_rho[i] |
| 194 | + end |
| 195 | + end |
| 196 | + end |
| 197 | + end |
| 198 | + for k in axes(daij_dsigmak, 3) |
| 199 | + G[k] = 2 * dot(daij_dsigmak[:, :, k], A) |
| 200 | + end |
| 201 | + for k in axes(daij_drhok, 3) |
| 202 | + G[L + k] = 2 * dot(daij_drhok[:, :, k], A) |
| 203 | + end |
| 204 | + end |
| 205 | + if !isnothing(F) |
| 206 | + return norm(A)^2 |
| 207 | + end |
| 208 | +end |
| 209 | + |
| 210 | +end # module |
0 commit comments