SimpleDiscontinuousGalerkin.jl API

SimpleDiscontinuousGalerkin

SimpleDiscontinuousGalerkin.jl is a Julia package that implements some basic discontinuous Galerkin (DG) methods for the solution of hyperbolic partial differential equations (PDEs). The package is designed to be simple and easy to use and understand. It is intended for educational purposes and to provide a starting point for more complex DG methods. For a more comprehensive and advanced implementation of DG methods, we recommend using the package Trixi.jl. This package can be understood as a minimalistic version of Trixi.jl, which is designed to be easy to understand and modify. Many design concepts are inspired by Trixi.jl, but the implementation is much simpler and more straightforward.

SimpleDiscontinuousGalerkin.jl builds on the foundations of the package SummationByPartsOperators.jl.

See also: SimpleDiscontinuousGalerkin.jl

BurgersEquation1D

The inviscid Burgers' equation

\[\partial_t u + \frac{1}{2} \partial_1 u^2 = 0\]

in one space dimension.

flux_ec(u_ll, u_rr, equations::BurgersEquation1D)

Entropy-conserving, symmetric flux for the inviscid Burgers' equation.

\[F(u_L, u_R) = \frac{u_L^2 + u_L u_R + u_R^2}{6}\]

initial_condition_convergence_test(x, t, equations::BurgersEquation1D)

A smooth initial condition used for convergence tests.

source_terms_convergence_test(u, x, t, equations::BurgersEquation1D)

Source terms used for convergence tests in combination with initial_condition_convergence_test.

CompressibleEulerEquations1D(gamma)

The compressible Euler equations

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho v_1 \\ \rho e_{\text{total}} \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + p \\ (\rho e_{\text{total}} +p) v_1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}\]

for an ideal gas with ratio of specific heats gamma in one space dimension. Here, $\rho$ is the density, $v_1$ the velocity, $e_{\text{total}}$ the specific total energy, and

\[p = (\gamma - 1) \left( \rho e_{\text{total}} - \frac{1}{2} \rho v_1^2 \right)\]

the pressure.

flux_kennedy_gruber(u_ll, u_rr, equations::CompressibleEulerEquations1D)

Kinetic energy preserving two-point flux by

• Kennedy and Gruber (2008) Reduced aliasing formulations of the convective terms within the Navier-Stokes equations for a compressible fluid DOI: 10.1016/j.jcp.2007.09.020

flux_ranocha(u_ll, u_rr, equations::CompressibleEulerEquations1D)

Entropy conserving and kinetic energy preserving two-point flux by

• Hendrik Ranocha (2018) Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods for Hyperbolic Balance Laws PhD thesis, TU Braunschweig

See also

• Hendrik Ranocha (2020) Entropy Conserving and Kinetic Energy Preserving Numerical Methods for the Euler Equations Using Summation-by-Parts Operators Proceedings of ICOSAHOM 2018

initial_condition_convergence_test(x, t, equations::CompressibleEulerEquations1D)

A smooth initial condition used for convergence tests in combination with source_terms_convergence_test (and BoundaryConditionDirichlet(initial_condition_convergence_test) in non-periodic domains).

initial_condition_density_wave(x, t, equations::CompressibleEulerEquations1D)

A sine wave in the density with constant velocity and pressure; reduces the compressible Euler equations to the linear advection equations. This setup is the test case for stability of EC fluxes from the paper

• Gregor J. Gassner, Magnus Svärd, Florian J. Hindenlang (2020) Stability issues of entropy-stable and/or split-form high-order schemes arXiv: 2007.09026

with the following parameters

• domain [-1, 1]
• mesh = 4x4
• polydeg = 5

initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations1D)

A weak blast wave taken from

• Sebastian Hennemann, Gregor J. Gassner (2020) A provably entropy stable subcell shock capturing approach for high order split form DG arXiv: 2008.12044

source_terms_convergence_test(u, x, t, equations::CompressibleEulerEquations1D)

Source terms used for convergence tests in combination with initial_condition_convergence_test (and BoundaryConditionDirichlet(initial_condition_convergence_test) in non-periodic domains).

AbstractEquations{NDIMS, NVARS}

An abstract supertype of specific equations such as the linear advection equation. The type parameters encode the number of spatial dimensions (NDIMS) and the number of primary variables (NVARS) of the physics model.

cons2cons(u, equations)

Return the conservative variables u. While this function is as trivial as identity, it is also as useful.

cons2entropy(u, equations)

Return the entropy variables from the conservative variables u for the given system of equations. The entropy variables are defined as the derivative of the entropy with respect to the conservative variables.

cons2prim(u, equations)

Return the primitive variables from the conservative variables u.

default_analysis_errors(equations)

Default analysis errors used by the AnalysisCallback.

default_analysis_integrals(equations)

Default analysis integrals used by the AnalysisCallback.

eachvariable(equations::AbstractEquations)

Return an iterator over the indices that specify the location in relevant data structures for the variables in equations. In particular, not the variables themselves are returned.

entropy(u, equations)

Return the entropy of the conservative variables u for the given system of equations.

entropy_timederivative

The semi-discrete time derivative of the entropy, which can be used in the AnalysisCallback to compute the time derivative of the entropy.

get_name(equations::AbstractEquations)

Return the canonical, human-readable name for the given system of equations.

Examples

julia> SimpleDiscontinuousGalerkin.get_name(LinearAdvectionEquation1D(1.0))
"LinearAdvectionEquation1D"

mass(u, equations)

Return the mass of the conservative variables u for the given system of equations.

prim2cons(q, equations)

Return the conservative variables from the primitive variables q.

LinearAdvectionEquation1D

The linear scalar advection equation

\[\partial_t u + a \partial_1 u = 0\]

in one space dimension with constant velocity a.

initial_condition_convergence_test(x, t, equations::LinearAdvectionEquation1D)

A smooth initial condition used for convergence tests.

MaxwellEquations1D(c = 299_792_458.0)

The Maxwell equations of electro dynamics

\[\frac{\partial}{\partial t} \begin{pmatrix} E \\ B \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} c^2 B \\ E \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\]

in one dimension with speed of light c = 299792458 m/s (in vacuum). In one dimension the Maxwell equations reduce to a wave equation. The orthogonal magnetic (e.g.B_y) and electric field (E_z) propagate as waves through the domain in x-direction. For reference, see

• e.g. p.15 of Numerical Methods for Conservation Laws: From Analysis to Algorithms https://doi.org/10.1137/1.9781611975109

• or equation (1) in https://inria.hal.science/hal-01720293/document

initial_condition_convergence_test(x, t, equations::MaxwellEquations1D)

A smooth initial condition used for convergence tests.

flux_hll

See FluxHLL.

flux_lax_friedrichs

See FluxLaxFriedrichs.

DissipationLocalLaxFriedrichs(max_abs_speed=max_abs_speed)

Create a local Lax-Friedrichs dissipation operator where the maximum absolute wave speed is estimated as max_abs_speed(u_ll, u_rr, equations), defaulting to max_abs_speed.

FluxHLL(min_max_speed=min_max_speed_davis)

Create an HLL (Harten, Lax, van Leer) numerical flux where the minimum and maximum wave speeds are estimated as λ_min, λ_max = min_max_speed(u_ll, u_rr, orientation_or_normal_direction, equations), defaulting to min_max_speed_davis. Original paper:

• Amiram Harten, Peter D. Lax, Bram van Leer (1983) On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws DOI: 10.1137/1025002

FluxLaxFriedrichs(max_abs_speed=max_abs_speed)

Local Lax-Friedrichs (Rusanov) flux with maximum wave speed estimate provided by max_abs_speed, cf. DissipationLocalLaxFriedrichs and max_abs_speed.

FluxPlusDissipation(numerical_flux, dissipation)

Combine a numerical_flux with a dissipation operator to create a new numerical flux.

RiemannProblem(u_ll, u_rr)

Create a Riemann problem with left and right states u_ll and u_rr, which can be solved with a RiemannSolver. This is used for the flux_godunov numerical flux.

RiemannSolver(prob, equations)

An exact Riemann solver of the RiemannProblem prob for equations.

RiemannSolverSolution

A solution of a Riemann problem, which is a vector of states at different times. Is returned when solveing a RiemannProblem with a RiemannSolver.

flux_central(u_ll, u_rr, equations::AbstractEquations)

The classical central numerical flux f((u_ll) + f(u_rr)) / 2. When this flux is used as volume flux, the discretization is equivalent to the classical weak form DG method (except floating point errors).

flux_godunov(u_ll, u_rr, equations)

Numerical flux using the RiemannSolver to solve Riemann problems with left and right states u_ll and u_rr for the given equations exactly.

max_abs_speed(u_ll, u_rr, equations)

Simple and fast estimate of the maximal wave speed of the Riemann problem with left and right states u_ll, u_rr, based only on the local wave speeds associated to u_ll and u_rr.

min_max_speed_davis(u_ll, u_rr, equations)

Simple and fast estimates of the minimal and maximal wave speed of the Riemann problem with left and right states u_ll, u_rr, usually based only on the local wave speeds associated to u_ll and u_rr.

• S.F. Davis (1988) Simplified Second-Order Godunov-Type Methods DOI: 10.1137/0909030

See eq. (10.38) from

• Eleuterio F. Toro (2009) Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction DOI: 10.1007/b79761

See also FluxHLL, min_max_speed_davis.

InhomogeneousMesh{NDIMS, RealT}
InhomogeneousMesh(coordinates)

Create an inhomogeneous one-dimensional mesh from the given coordinates.

Mesh

Struct that holds the information for a simple homogeneous one-dimensional mesh.

Mesh(xmin, xmax, N_elements)

Create a simple homogeneous one-dimensional mesh from xmin to xmax with N_elements elements.

OversetGridMesh{NDIMS, RealT, MeshLeft, MeshRight}
OversetGridMesh(mesh_left, mesh_right)

Create an overset grid (Chimera) mesh that combines two meshes, mesh_left and mesh_right, which have an overlap region.

element_spacing(mesh::AbstractMesh)

Return the spacing of the elements in the mesh. This is the length of each element, which is the same for all elements in a homogeneous Mesh and can be different for each element in an InhomogeneousMesh.

element_spacing(mesh::AbstractMesh, element)

Return the length of the element element in the mesh.

left_element_boundary(mesh::AbstractMesh, element)

Return the left boundary coordinate of the element element in the mesh.

total_volume(mesh::AbstractMesh)

Return the total volume of the mesh, which is the length of the domain in one dimension.

xmax(mesh::AbstractMesh)

Return the maximum coordinate of the mesh.

xmin(mesh::AbstractMesh)

Return the minimum coordinate of the mesh.

boundary_condition_do_nothing = BoundaryConditionDoNothing()

Imposing no boundary condition just evaluates the flux at the inner state.

boundary_condition_periodic = BoundaryConditionPeriodic()

A singleton struct indicating periodic boundary conditions.

BoundaryConditionDirichlet(boundary_value_function)

Create a Dirichlet boundary condition that uses the function boundary_value_function to specify the values at the boundary. This can be used to create a boundary condition that specifies exact boundary values by passing the exact solution of the equation. The passed boundary value function will be called with the same arguments as an initial condition function is called, i.e., as

boundary_value_function(x, t, equations)

where x specifies the coordinates, t is the current time, and equation is the corresponding system of equations.

Examples

julia> BoundaryConditionDirichlet(initial_condition_convergence_test)

DG(; basis, surface_integral, volume_integral)

Create a discontinuous Galerkin method. If basis isa LegendreDerivativeOperator, this creates a DGSEM.

eachnode(solver::DG, element)

Return an iterator over the indices that specify the location in relevant data structures for the nodes in a specific element in solver. In particular, not the nodes themselves are returned.

get_variable(u, v, ::DG)

Return the solution belonging to the variable v of the solution u at one time step as a vector at every node across all elements.

DGSEM(; RealT=Float64, polydeg::Integer,
        surface_flux=flux_central,
        surface_integral=SurfaceIntegralWeakForm(surface_flux),
        volume_integral=VolumeIntegralWeakForm())

Create a discontinuous Galerkin spectral element method (DGSEM) using a LegendreDerivativeOperator with polynomials of degree polydeg.

FDSBP(D; RealT=Float64,
         surface_flux=flux_central,
         surface_integral=SurfaceIntegralWeakForm(surface_flux),
         volume_integral=VolumeIntegralWeakForm())

Create a discontinuous Galerkin method using a summation-by-parts operator D from SummationByPartsOperators.jl. This is similar to the DGSEM, but uses a general derivative operator instead of a Legendre derivative operator.

PerElementBasis{BasisType}

Basis, which can hold a different SBP operator for each element. This is used in PerElementFDSBP to allow for different bases on each element.

PerElementFDSBP(bases::Vector{BasisType};
                surface_flux = flux_central,
                surface_integral = SurfaceIntegralWeakForm(surface_flux),
                volume_integral = VolumeIntegralWeakForm()) where BasisType

Create a discontinuous Galerkin method using different bases for each element. This is like FDSBP, but allows for a different SBP operator on each element. See also: PerElementBasis.

SurfaceIntegralStrongForm(surface_flux=flux_central, surface_flux_boundary=surface_flux)

The classical strong form surface integral type for FD/DG methods. It uses surface_flux for the interior fluxes and surface_flux_boundary for the boundary fluxes.

See also VolumeIntegralStrongForm.

SurfaceIntegralWeakForm(surface_flux=flux_central, surface_flux_boundary=surface_flux)

The classical weak form surface integral type for DG methods as explained in standard textbooks. It uses surface_flux for the interior fluxes and surface_flux_boundary for the boundary fluxes.

See also VolumeIntegralWeakForm.

References

• Kopriva (2009) Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engineers doi: 10.1007/978-90-481-2261-5
• Hesthaven, Warburton (2007) Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications doi: 10.1007/978-0-387-72067-8

VolumeIntegralFluxDifferencing(volume_flux=flux_central)

Volume integral type for DG methods based on SBP operators and flux differencing using a symmetric two-point volume_flux. This volume_flux needs to satisfy the interface of numerical fluxes.

To be used together with SurfaceIntegralWeakForm.

VolumeIntegralFluxDifferencingStrongForm(volume_flux=flux_central)

Volume integral type for DG methods based on SBP operators and flux differencing using a symmetric two-point volume_flux. This volume_flux needs to satisfy the interface of numerical fluxes.

This is the strong formulation, which means it should be used together with SurfaceIntegralStrongForm.

VolumeIntegralStrongForm()

The classical strong form volume integral type for FD/DG methods.

VolumeIntegralWeakForm()

The classical weak form volume integral type for DG methods as explained in standard textbooks.

References

• Kopriva (2009) Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engineers doi: 10.1007/978-90-481-2261-5
• Hesthaven, Warburton (2007) Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications doi: 10.1007/978-0-387-72067-8

Semidiscretization

A struct containing everything needed to describe a spatial semidiscretization of an equation.

Semidiscretization(mesh, equations, initial_condition, solver;
                   source_terms = nothing,
                   boundary_conditions)

Construct a semidiscretization of a PDE.

grid(semi)

Get the grid of a semidiscretization.

flat_grid(semi)

Return a vector of the coordinates of all nodes in semi, flattened across all elements. This is useful for plotting or other operations that require a single vector of coordinates.

get_variable(u, v, semi)

Return the solution belonging to the variable v of the solution u at one time step as a vector at every node across all elements.

semidiscretize(semi::Semidiscretization, tspan)

Wrap the semidiscretization semi as an ODE problem in the time interval tspan that can be passed to solve from the SciML ecosystem.

AnalysisCallback(semi; interval=0,
                       extra_analysis_errors=Symbol[],
                       extra_analysis_integrals=(),
                       io=stdout)

Analyze a numerical solution every interval time steps. The L2- and the L∞-norm for each component are computed by default. Additional errors can be computed, e.g. by passing extra_analysis_errors = (:conservation_error,).

Further scalar functions func in extra_analysis_integrals are applied to the numerical solution and integrated over the computational domain. Some examples for this are mass, and entropy. You can also write your own function with the same signature as the examples listed above and pass it via extra_analysis_integrals. The computed errors and integrals are saved for each timestep and can be obtained by calling errors and integrals.

During the Simulation, the AnalysisCallback will print information to io.

errors(analysis_callback)

Return the computed errors for each timestep as a named tuple. The shape of each entry is (nvariables, ntimesteps).

integrals(analysis_callback)

Return the computed integrals for each time step as a named tuple.

tstops(analysis_callback)

Return the time values that correspond to the saved values of the errors and integrals.

StepsizeCallback(; cfl=1.0, interval = 1)

Set the time step size according to a CFL condition with CFL number cfl if the time integration method isn't adaptive itself.

The supplied keyword argument cfl must be either a Real number or a function of time t returning a Real number. By default, the timestep will be adjusted at every step. For different values of interval, the timestep will be adjusted every interval steps.

SummaryCallback(io::IO = stdout)

Create and return a callback that resets the timer at the beginning of a simulation and prints the timer values at the end of the simulation.

calc_mean_convergence(eocs)

Calculate the mean convergence rates from the given experimental orders of convergence eocs. The eocs are expected to be in the format returned by convergence_test, i.e., a Dict where the keys are the error types (e.g., :l2, :linf) and the values are matrices with the EOCs for each variable in the columns and the iterations in the rows. Returns a Dict with the same keys as eocs and the mean convergence rates for all variables as values.

convergence_test([mod::Module=Main,] example, iterations; io::IO = stdout, kwargs...)
convergence_test([mod::Module=Main,] example, Ns::AbstractVector; io::IO = stdout, kwargs...)

Run multiple simulations using the setup given in example and compute the experimental order of convergence (EOC) in the $L^2$ and $L^\infty$ norm. If iterations is passed as integer, in each iteration, the resolution of the respective mesh will be doubled. If Ns is passed as vector, the simulations will be run for each value of Ns. Additional keyword arguments kwargs... and the optional module mod are passed directly to trixi_include.

Adjusted from Trixi.jl.

default_example()

Return the path to an example that can be used to quickly see SimpleDiscontinuousGalerkin.jl in action. See also examples_dir.

Copied from Trixi.jl.

examples_dir()

Return the directory where the example files provided by SimpleDiscontinuousGalerkin.jl are located. If SimpleDiscontinuousGalerkin.jl is installed as a regular package (with ]add SimpleDiscontinuousGalerkin), these files are read-only and should not be modified. To find out which files are available, use, e.g., readdir.

Copied from Trixi.jl.

Examples

readdir(examples_dir())