Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Docstrings for some methods in basis Lobatto-Legendre #1874

Merged
merged 3 commits into from
Mar 15, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/literate/src/files/scalar_linear_advection_1d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ integral = sum(nodes.^3 .* weights)
# To approximate the solution, we need to get the polynomial coefficients $\{u_j^{Q_l}\}_{j=0}^N$
# for every element $Q_l$.

# After defining all nodes, we can implement the spatial coordinate $x$ and its initial value $u0$
# After defining all nodes, we can implement the spatial coordinate $x$ and its initial value $u0 = u(t_0)$
# for every node.
x = Matrix{Float64}(undef, length(nodes), n_elements)
for element in 1:n_elements
Expand Down
77 changes: 71 additions & 6 deletions src/solvers/dgsem/basis_lobatto_legendre.jl
Original file line number Diff line number Diff line change
Expand Up @@ -404,7 +404,8 @@ function calc_dsplit(nodes, weights)
return dsplit
end

# Calculate the polynomial derivative matrix D
# Calculate the polynomial derivative matrix D.
# This implements algorithm 37 "PolynomialDerivativeMatrix" from Kopriva's book.
function polynomial_derivative_matrix(nodes)
n_nodes = length(nodes)
d = zeros(n_nodes, n_nodes)
Expand All @@ -421,6 +422,7 @@ function polynomial_derivative_matrix(nodes)
end

# Calculate and interpolation matrix (Vandermonde matrix) between two given sets of nodes
# See algorithm 32 "PolynomialInterpolationMatrix" from Kopriva's book.
function polynomial_interpolation_matrix(nodes_in, nodes_out,
baryweights_in = barycentric_weights(nodes_in))
n_nodes_in = length(nodes_in)
Expand All @@ -433,6 +435,7 @@ function polynomial_interpolation_matrix(nodes_in, nodes_out,
return vandermonde
end

# This implements algorithm 32 "PolynomialInterpolationMatrix" from Kopriva's book.
function polynomial_interpolation_matrix!(vandermonde,
nodes_in, nodes_out,
baryweights_in)
Expand Down Expand Up @@ -463,7 +466,19 @@ function polynomial_interpolation_matrix!(vandermonde,
return vandermonde
end

# Calculate the barycentric weights for a given node distribution.
"""
barycentric_weights(nodes)

Calculate the barycentric weights for a given node distribution, i.e.,
```math
w_j = \\frac{1}{ \\prod_{k \\neq j} \\left( x_j - x_k \\right ) }
```

For details, see (especially Section 3)
- Jean-Paul Berrut and Lloyd N. Trefethen (2004).
Barycentric Lagrange Interpolation.
[DOI:10.1137/S0036144502417715](https://doi.org/10.1137/S0036144502417715)
"""
function barycentric_weights(nodes)
n_nodes = length(nodes)
weights = ones(n_nodes)
Expand Down Expand Up @@ -494,12 +509,31 @@ function calc_lhat(x, nodes, weights)
return lhat
end

# Calculate Lagrange polynomials for a given node distribution.
"""
lagrange_interpolating_polynomials(x, nodes, wbary)

Calculate Lagrange polynomials for a given node distribution with
associated barycentric weights `wbary` at a given point `x` on the
reference interval ``[-1, 1]``.

This returns all ``l_j(x)``, i.e., the Lagrange polynomials for each node ``x_j``.
Thus, to obtain the interpolating polynomial ``p(x)`` at ``x``, one has to
multiply the Lagrange polynomials with the nodal values ``u_j`` and sum them up:
``p(x) = \\sum_{j=1}^{n} u_j l_j(x)``.

For details, see e.g. Section 2 of
- Jean-Paul Berrut and Lloyd N. Trefethen (2004).
Barycentric Lagrange Interpolation.
[DOI:10.1137/S0036144502417715](https://doi.org/10.1137/S0036144502417715)
"""
function lagrange_interpolating_polynomials(x, nodes, wbary)
n_nodes = length(nodes)
polynomials = zeros(n_nodes)

for i in 1:n_nodes
# Avoid division by zero when `x` is close to node by using
# the Kronecker-delta property at nodes
# of the Lagrange interpolation polynomials.
if isapprox(x, nodes[i], rtol = eps(x))
polynomials[i] = 1
return polynomials
Expand All @@ -518,6 +552,17 @@ function lagrange_interpolating_polynomials(x, nodes, wbary)
return polynomials
end

"""
gauss_lobatto_nodes_weights(n_nodes::Integer)

Computes nodes ``x_j`` and weights ``w_j`` for the (Legendre-)Gauss-Lobatto quadrature.
This implements algorithm 25 "GaussLobattoNodesAndWeights" from the book

- David A. Kopriva, (2009).
Implementing spectral methods for partial differential equations:
Algorithms for scientists and engineers.
[DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
"""
# From FLUXO (but really from blue book by Kopriva)
function gauss_lobatto_nodes_weights(n_nodes::Integer)
# From Kopriva's book
Expand Down Expand Up @@ -585,7 +630,7 @@ function gauss_lobatto_nodes_weights(n_nodes::Integer)
return nodes, weights
end

# From FLUXO (but really from blue book by Kopriva)
# From FLUXO (but really from blue book by Kopriva, algorithm 24)
function calc_q_and_l(N::Integer, x::Float64)
L_Nm2 = 1.0
L_Nm1 = x
Expand All @@ -609,7 +654,17 @@ function calc_q_and_l(N::Integer, x::Float64)
end
calc_q_and_l(N::Integer, x::Real) = calc_q_and_l(N, convert(Float64, x))

# From FLUXO (but really from blue book by Kopriva)
"""
gauss_nodes_weights(n_nodes::Integer)

Computes nodes ``x_j`` and weights ``w_j`` for the Gauss-Legendre quadrature.
This implements algorithm 23 "LegendreGaussNodesAndWeights" from the book

- David A. Kopriva, (2009).
Implementing spectral methods for partial differential equations:
Algorithms for scientists and engineers.
[DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
"""
function gauss_nodes_weights(n_nodes::Integer)
# From Kopriva's book
n_iterations = 10
Expand Down Expand Up @@ -666,7 +721,17 @@ function gauss_nodes_weights(n_nodes::Integer)
end
end

# From FLUXO (but really from blue book by Kopriva)
"""
legendre_polynomial_and_derivative(N::Int, x::Real)

Computes the Legendre polynomial of degree `N` and its derivative at `x`.
This implements algorithm 22 "LegendrePolynomialAndDerivative" from the book

- David A. Kopriva, (2009).
Implementing spectral methods for partial differential equations:
Algorithms for scientists and engineers.
[DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
"""
function legendre_polynomial_and_derivative(N::Int, x::Real)
if N == 0
poly = 1.0
Expand Down
Loading