add more docs

tristanmontoya · Oct 6, 2024 · adb1905 · adb1905 · tristanmontoya · Oct 6, 2024
1 parent 55e4c8a
commit adb1905
Show file tree

Hide file tree

Showing 23 changed files with 1,525 additions and 1,479 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "StableSpectralElements"
 uuid = "fb992021-99c7-4c2d-a14b-5e48ac4045b2"
 authors = ["Tristan Montoya <montoya.tristan@gmail.com>"]
-version = "0.2.13-DEV"
+version = "0.2.13"
 
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"

diff --git a/docs/src/ConservationLaws.md b/docs/src/ConservationLaws.md
@@ -1,29 +1,25 @@
 # Module `ConservationLaws`
 
+The `ConservationLaws` module defines the systems of partial differential equations which are solved by StableSpectralElements.jl.
+
 ## Overview
 
-The equations to be solved are defined by subtypes of `AbstractConservationLaw` on which functions such as `physical_flux` and `numerical_flux` are dispatched. Objects of type `AbstractConservationLaw` contain two type parameters, `d` and `PDEType`, the former denoting the spatial dimension of the problem, which is inherited by all subtypes, and the latter being a subtype of `AbstractPDEType` denoting the particular type of PDE being solved, which is either `FirstOrder` or `SecondOrder`. Whereas first-order problems remove the dependence of the flux tensor on the solution gradient in order to obtain systems of the form
+The equations to be solved are defined by subtypes of [`AbstractConservationLaw`](@ref) on which functions such as `physical_flux` and `numerical_flux` are dispatched. Whereas first-order problems (i.e. subtypes of `AbstractConservationLaw{d, FirstOrder}`) remove the dependence of the flux tensor on the solution gradient in order to obtain systems of the form
 ```math
 \partial_t \underline{U}(\bm{x},t) + \bm{\nabla}_{\bm{x}} \cdot \underline{\bm{F}}(\underline{U}(\bm{x},t)) = \underline{0},
 ```
-second-order problems are treated by StableSpectralElements.jl as first-order systems of the form 
+second-order problems (i.e. subtypes of `AbstractConservationLaw{d, SecondOrder}`) are treated by StableSpectralElements.jl as first-order systems of the form 
 ```math
 \begin{aligned}
 \underline{\bm{Q}}(\bm{x},t) - \bm{\nabla}_{\bm{x}} \underline{U}(\bm{x},t) &= \underline{0},\\
 \partial_t \underline{U}(\bm{x},t) + \bm{\nabla}_{\bm{x}} \cdot \underline{\bm{F}}(\underline{U}(\bm{x},t), \underline{\bm{Q}}(\bm{x},t)) &= \underline{0}.
 \end{aligned}
 ```
-Currently, the linear advection and advection-diffusion equations, the inviscid and viscous Burgers' equations, and the compressible Euler equations are supported by StableSpectralElements.jl, but any system of the above form can in principle be implemented, provided that appropriate physical and numerical fluxes are defined.
+Currently, the linear advection and advection-diffusion equations ([`LinearAdvectionEquation`](@ref) and [`LinearAdvectionDiffusionEquation`](@ref)), the inviscid and viscous Burgers' equations ([`InviscidBurgersEquation`](@ref) and [`ViscousBurgersEquation`](@ref)), and the compressible Euler equations ([`EulerEquations`](@ref)) are supported by StableSpectralElements.jl, but any system of the above form can in principle be implemented, provided that appropriate physical and numerical fluxes are defined.
 
 ## Reference
 
-```@meta
-CurrentModule = ConservationLaws
-```
-```@docs
-    LinearAdvectionEquation
-    LinearAdvectionDiffusionEquation
-    InviscidBurgersEquation
-    ViscousBurgersEquation
-    EulerEquations
+```@autodocs
+Modules = [ConservationLaws]
+Order   = [:function, :type]
 ```
diff --git a/docs/src/Solvers.md b/docs/src/Solvers.md
@@ -1,19 +1,15 @@
 # Module `Solvers`
 
+The `Solvers` module implements the algorithms which evaluate the semi-discrete residual corresponding to the discretization of an [`AbstractConservationLaw`](@ref) in space using the operators and geometric information contained within the [`SpatialDiscretization`](@ref) to create an [`ODEProblem`](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/), which can be solved using your choice of time-marching method from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl).
+
 ## Overview
 
-StableSpectralElements.jl is based on a semi-discrete approach to the numerical solution of partial differential equations, in which the spatial discretization is performed first in order to obtain a system of ordinary differential equations of the form 
-```math
-\frac{\mathrm{d} }{\mathrm{d}t}\underline{u}(t) = \underline{R}(\underline{u}(t),t),
-```
-where $\underline{u}(t) \in \mathbb{R}^{N_p \cdot N_c \cdot N_e}$ is the global solution array containing the $N_p$ coefficients for each of the $N_c$ solution variables and each of the $N_e$ mesh elements. Such systems can be solved using standard time-marching methods from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) by computing the semi-discrete residual $\underline{R}(\underline{u}(t),t)$ using a Julia function with the following signature:
-```julia
-semi_discrete_residual(dudt::AbstractArray{Float64,3},
-                       u::AbstractArray{Float64, 3},
-                       solver::Solver,
-                       t::Float64)
-```
-The first parameter contains the time derivative to be computed in place, the second parameter is the current solution state, and the fourth parameter is the time $t$. The third parameter, which is of type `Solver`, contains all the information defining the spatial discretization as well as preallocated arrays used for temporary storage. The particular algorithm used for computing the semi-discrete residual is then dispatched based on the particular parametric subtype of `Solver` which is passed into the `semi_discrete_residual!` function.
+StableSpectralElements.jl is based on a semi-discrete or "method-of-lines" approach to the numerical solution of partial differential equations, in which the spatial discretization is performed first in order to obtain a system of ordinary differential equations of the form $\underline{u}'(t) = \underline{R}(\underline{u}(t),t)$,
+where $\underline{u}(t) \in \mathbb{R}^{N_p \cdot N_c \cdot N_e}$ is the global solution array containing the $N_p$ coefficients for each of the $N_c$ solution variables and each of the $N_e$ mesh elements. These systems can be solved using standard time-marching methods from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) by computing the semi-discrete residual $\underline{R}(\underline{u}(t),t)$ using the in-place function [`semi_discrete_residual!`](@ref).
 
 ## Reference
 
+```@autodocs
+Modules = [Solvers]
+Order   = [:function, :type]
+```
diff --git a/docs/src/SpatialDiscretizations.md b/docs/src/SpatialDiscretizations.md
@@ -1,5 +1,7 @@
 # Module `SpatialDiscretizations`
 
+The `SpatialDiscretizations` module defines the discretization on the reference element and provides the geometric data relating to the mesh and the mapping from reference to physical space.
+
 ## Overview
 Discretizations in StableSpectralElements.jl are constructed by first building a local approximation on a canonical reference element, denoted generically as $\hat{\Omega} \subset \mathbb{R}^d$, and using a bijective transformation $\bm{X}^{(\kappa)} : \hat{\Omega} \rightarrow \Omega^{(\kappa)}$ to construct the approximation on each physical element $\Omega^{(\kappa)} \subset \Omega$ of the mesh $\{ \Omega^{(\kappa)}\}_{\kappa \in \{1:N_e\}}$ in terms of the associated operators on the reference element. An example of such a mapping is shown below, where we also depict the collapsed coordinate transformation $\bm{\chi} : [-1,1]^d \to \hat{\Omega}$ which may be used to construct operators with a tensor-product structure on the reference simplex.
 
@@ -15,21 +17,12 @@ In order to define the different geometric reference elements, the subtypes `Lin
 \hat{\Omega}_{\mathrm{tet}} &= \big\{ \bm{\xi} \in [-1,1]^3 : \xi_1 + \xi_2 + \xi_3 \leq -1 \big\}.
 \end{aligned}
 ```
-These element types are used in the constructor for StableSpectralElements.jl's [`ReferenceApproximation`](@ref) type, along with a subtype of `AbstractApproximationType` ([`NodalTensor`](@ref), [`ModalTensor`](@ref), [`NodalMulti`] (@ref), [`ModalMulti`](@ref), [`NodalMultiDiagE`](@ref), or [`ModalMultiDiagE`](@ref)) specifying the nature of the local approximation. 
+These element types are used in the constructor for StableSpectralElements.jl's [`ReferenceApproximation`](@ref) type, along with a subtype of `AbstractApproximationType` ([`NodalTensor`](@ref), [`ModalTensor`](@ref), [`NodalMulti`](@ref), [`ModalMulti`](@ref), [`NodalMultiDiagE`](@ref), or [`ModalMultiDiagE`](@ref)) specifying the nature of the local approximation. 
 
-All the information used to define the spatial discretization on the physical domain $\Omega$ is contained within a `SpatialDiscretization` structure, which is constructed using a [`ReferenceApproximation`](@ref) and a `MeshData` from StartUpDG.jl, which are stored as the fields `reference_approximation` and `mesh`. When the constructor for a `SpatialDiscretization` is called, the grid metrics are computed and stored in the field `geometric_factors` of type `GeometricFactors`.
+All the information used to define the spatial discretization on the physical domain $\Omega$ is contained within a [`SpatialDiscretization`](@ref) structure, which is constructed using a [`ReferenceApproximation`](@ref) and a [`MeshData`](https://jlchan.github.io/StartUpDG.jl/stable/MeshData/) from StartUpDG.jl. When the constructor for a [`SpatialDiscretization`](@ref) is called, the grid metrics are computed and stored in a field of type [`GeometricFactors`](@ref).
 
 ## Reference
-
-```@meta
-CurrentModule = SpatialDiscretizations
-```
-```@docs
-    ReferenceApproximation
-    NodalTensor
-    ModalTensor
-    NodalMulti
-    ModalMulti
-    NodalMultiDiagE
-    ModalMultiDiagE
+```@autodocs
+Modules = [SpatialDiscretizations]
+Order   = [:function, :type]
 ```
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -5,7 +5,7 @@
 \partial_t \underline{U}(\bm{x},t) + \bm{\nabla}_{\bm{x}} \cdot \underline{\bm{F}}(\underline{U}(\bm{x},t), \bm{\nabla}_{\bm{x}}\underline{U}(\bm{x},t)) = \underline{0},
 ```
 for $t \in (0,T)$ with $T \in \mathbb{R}^+ $ and $\bm{x} \in \Omega \subset \mathbb{R}^d$, subject to appropriate initial and boundary conditions, where $\underline{U}(\bm{x},t)$ is the vector of solution variables and $\underline{\bm{F}}(\underline{U}(\bm{x},t),\bm{\nabla}_{\bm{x}}\underline{U}(\bm{x},t))$ is the flux tensor containing advective and/or diffusive contributions. 
-These equations are spatially discretized on curvilinear unstructured grids using energy-stable and entropy-stable discontinuous spectral element methods in order to generate `ODEProblem` objects suitable for time integration using [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) within the [SciML](https://sciml.ai/) ecosystem. StableSpectralElements.jl also includes postprocessing tools employing [WriteVTK.jl](https://github.com/jipolanco/WriteVTK.jl) for generating `.vtu` files, allowing for visualization of high-order numerical solutions on unstructured grids using [ParaView](https://www.paraview.org/) or other software. Shared-memory parallelization is supported through multithreading.
+These equations are spatially discretized on curvilinear unstructured grids using energy-stable and entropy-stable discontinuous spectral element methods in order to generate [`ODEProblem`](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/) objects suitable for time integration using [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) within the [SciML](https://sciml.ai/) ecosystem. StableSpectralElements.jl also includes postprocessing tools employing [WriteVTK.jl](https://github.com/jipolanco/WriteVTK.jl) for generating `.vtu` files, allowing for visualization of high-order numerical solutions on unstructured grids using [ParaView](https://www.paraview.org/) or other software. Shared-memory parallelization is supported through multithreading.
 
 The functionality provided by [StartUpDG.jl](https://github.com/jlchan/StartUpDG.jl) for the handling of mesh data structures, polynomial basis functions, and quadrature nodes is employed throughout this package. Moreover, StableSpectralElements.jl implements dispatched strategies for semi-discrete operator evaluation using [LinearMaps.jl](https://github.com/JuliaLinearAlgebra/LinearMaps.jl), allowing for the efficient matrix-free application of tensor-product operators.
 
@@ -68,5 +68,9 @@ If you use StableSpectralElements.jl in your research, please cite the following
 }
 ```
 
+The following repositories are associated with the above papers, and contain the code required to reproduce the results presented therein:
+- <https://github.com/tristanmontoya/ReproduceSBPSimplex>
+- <https://github.com/tristanmontoya/ReproduceEntropyStableDSEM>
+
 ## License
 This software is released under the [GPLv3 license](https://www.gnu.org/licenses/gpl-3.0.en.html).
diff --git a/examples/advection_2d.ipynb b/examples/advection_2d.ipynb
diff --git a/examples/burgers_1d.ipynb b/examples/burgers_1d.ipynb
diff --git a/examples/figures/burgers_solution.gif b/examples/figures/burgers_solution.gif
diff --git a/src/Analysis/benchmark.jl b/src/Analysis/benchmark.jl
@@ -1,5 +1,7 @@
 
-"""This removes the threading and considers just one element. Note that this probably will not work with the Euler equations as the facet states at adjacent elements are left undefined and thus may lead to non-physical states when used to compute the fluxes"""
+# This removes the threading and considers just one element. Note that this probably will 
+# not work with the Euler equations as the facet states at adjacent elements are left 
+# undefined and thus may lead to non-physical states when used to compute the fluxes.
 @views @timeit "semi-disc. residual" function rhs_benchmark!(dudt::AbstractArray{Float64,
             3},
         u::AbstractArray{Float64, 3},

diff --git a/src/Analysis/conservation.jl b/src/Analysis/conservation.jl
@@ -1,9 +1,7 @@
 abstract type ConservationAnalysis <: AbstractAnalysis end
 abstract type AbstractConservationAnalysisResults <: AbstractAnalysisResults end
 
-"""
-Evaluate change in ∫udx  
-"""
+# Evaluate change in ∫udx  
 struct PrimaryConservationAnalysis{V_type} <: ConservationAnalysis
     WJ::Vector{Diagonal{Float64, Vector{Float64}}}
     N_c::Int
@@ -13,9 +11,7 @@ struct PrimaryConservationAnalysis{V_type} <: ConservationAnalysis
     dict_name::String
 end
 
-"""
-Evaluate change in  ∫½u²dx  
-"""
+# Evaluate change in  ∫½u²dx  
 struct EnergyConservationAnalysis{V_type, MassSolver} <: ConservationAnalysis
     mass_solver::MassSolver
     N_c::Int

diff --git a/src/Analysis/refinement.jl b/src/Analysis/refinement.jl
@@ -1,4 +1,4 @@
-"""Analyze results from grid refinement studies"""
+# Analyze results from grid refinement studies
 struct RefinementAnalysis{ExactSolution} <: AbstractAnalysis
     exact_solution::ExactSolution
     sequence_path::String

diff --git a/src/ConservationLaws/ConservationLaws.jl b/src/ConservationLaws/ConservationLaws.jl
@@ -36,12 +36,19 @@ export physical_flux,
     NoTwoPointFlux,
     ExactSolution
 
+@doc raw"""
+    AbstractConservationLaw{d, PDEType, N_c}
+
+Abstract type for a conservation law, where `d` is the number of spatial dimensions,
+`PDEType <: AbstractPDEType` is either `FirstOrder` or `SecondOrder`, and `N_c` is the number of conservative
+variables.
+"""
 abstract type AbstractConservationLaw{d, PDEType, N_c} end
 abstract type AbstractPDEType end
-struct FirstOrder <: AbstractPDEType end
-struct SecondOrder <: AbstractPDEType end
+struct FirstOrder <: AbstractPDEType end # PDE with only first derivatives
+struct SecondOrder <: AbstractPDEType end # PDE with first and second derivatives
 
-"""First-order numerical fluxes"""
+# First-order numerical fluxes
 abstract type AbstractInviscidNumericalFlux end
 struct NoInviscidFlux <: AbstractInviscidNumericalFlux end
 struct LaxFriedrichsNumericalFlux <: AbstractInviscidNumericalFlux
@@ -54,12 +61,12 @@ struct EntropyConservativeNumericalFlux <: AbstractInviscidNumericalFlux end
 struct CentralNumericalFlux <: AbstractInviscidNumericalFlux end
 LaxFriedrichsNumericalFlux() = LaxFriedrichsNumericalFlux(1.0)
 
-"""Second-order numerical fluxes"""
+# Second-order numerical fluxes
 abstract type AbstractViscousNumericalFlux end
 struct BR1 <: AbstractViscousNumericalFlux end
 struct NoViscousFlux <: AbstractViscousNumericalFlux end
 
-"""Two-point fluxes (for split forms and entropy-stable schemes)"""
+# Two-point fluxes (for split forms and entropy-stable schemes)
 abstract type AbstractTwoPointFlux end
 struct ConservativeFlux <: AbstractTwoPointFlux end
 struct EntropyConservativeFlux <: AbstractTwoPointFlux end
@@ -120,13 +127,13 @@ end
     end
 end
 
-""" 
-Algorithm based on the Taylor series trick from Ismail and Roe (2009). There are further optimizations that could be made, but let's leave it like this for now.
-"""
+# Algorithm based on the Taylor series trick from Ismail and Roe (2009). There are further 
+# optimizations that could be made, but let's leave it like this for now.
 @inline function logmean(x::Float64, y::Float64)
     # f = (y/x - 1) / (y/x + 1)
     #   = (y - x) / (x + y)
-    # rearrange to avoid divisions using trick from https://trixi-framework.github.io/Trixi.jl/stable/reference-trixi/#Trixi.ln_mean
+    # rearrange to avoid divisions using trick from 
+    # https://trixi-framework.github.io/Trixi.jl/stable/reference-trixi/#Trixi.ln_mean
     f² = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y)
     if f² < 1.0e-4
         return (x + y) * 105 / (210 + f² * (70 + f² * (42 + f² * 30)))
@@ -148,7 +155,7 @@ end
     end
 end
 
-"""Generic structure for exact solution to PDE (may be deprecated in future versions)"""
+# Generic structure for exact solution to PDE (may be deprecated in future versions)
 struct ExactSolution{d, ConservationLaw, InitialData, SourceTerm} <: AbstractGridFunction{d}
     conservation_law::ConservationLaw
     initial_data::InitialData

diff --git a/src/ConservationLaws/burgers.jl b/src/ConservationLaws/burgers.jl
@@ -44,14 +44,11 @@ end
 
 const BurgersType{d} = Union{InviscidBurgersEquation{d}, ViscousBurgersEquation{d}}
 
+# Default constructors
 InviscidBurgersEquation() = InviscidBurgersEquation((1.0,))
 ViscousBurgersEquation(b::Float64) = ViscousBurgersEquation((1.0,), b)
 
-"""
-Evaluate the flux for the inviscid Burgers' equation
-
-`F(u) = a ½u^2`
-"""
+# Evaluate the flux F(u) = a u^2/2 for the inviscid Burgers' equation
 function physical_flux!(f::AbstractArray{Float64, 3},
         conservation_law::InviscidBurgersEquation{d},
         u::AbstractMatrix{Float64}) where {d}
@@ -60,11 +57,7 @@ function physical_flux!(f::AbstractArray{Float64, 3},
     end
 end
 
-"""
-Evaluate the flux for the viscous Burgers' equation
-
-`F(u,q) = a ½u^2 - bq`
-"""
+# Evaluate the flux F(u,q) = a u^2/2 - bq for the viscous Burgers' equation 
 function physical_flux!(f::AbstractArray{Float64, 3},
         conservation_law::ViscousBurgersEquation{d},
         u::AbstractMatrix{Float64},
@@ -75,11 +68,8 @@ function physical_flux!(f::AbstractArray{Float64, 3},
     end
 end
 
-"""
-Evaluate the interface normal solution for the viscous Burgers' equation using the BR1 approach
-
-`U*(u⁻, u⁺, n) = ½(u⁻ + u⁺)n`
-"""
+# Evaluate the interface normal solution U*(u⁻, u⁺, n) = ½(u⁻ + u⁺)n 
+# for the viscous Burgers' equation using the BR1 approach 
 @inline function numerical_flux!(u_nstar::AbstractArray{Float64, 3},
         ::ViscousBurgersEquation{d},
         ::BR1,
@@ -93,11 +83,8 @@ Evaluate the interface normal solution for the viscous Burgers' equation using t
     end
 end
 
-"""
-Evaluate the numerical flux for the viscous Burgers' equation using the BR1 approach
-
-F*(u⁻, u⁺, q⁻, q⁺, n) = ½(F²(u⁻,q⁻) + F²(u⁺, q⁺))⋅n
-"""
+# Evaluate the numerical flux for the viscous Burgers' equation using the BR1 approach
+# F*(u⁻, u⁺, q⁻, q⁺, n) = ½(F²(u⁻,q⁻) + F²(u⁺, q⁺))⋅n
 @inline function numerical_flux!(f_star::AbstractMatrix{Float64},
         conservation_law::ViscousBurgersEquation{d},
         ::BR1,