Merge branch 'main' into sc/glm_speed_coupled_semis

docs/literate/src/files/scalar_linear_advection_1d.jl

@@ -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

src/solvers/dgsem/basis_lobatto_legendre.jl

@@ -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)
@@ -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)
@@ -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)
@@ -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)
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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