Trixi.jl API

Trixi

Trixi.jl is a numerical simulation framework for hyperbolic conservation laws. A key objective for the framework is to be useful to both scientists and students. Therefore, next to having an extensible design with a fast implementation, Trixi.jl is focused on being easy to use for new or inexperienced users, including the installation and postprocessing procedures.

To get started, run your first simulation with Trixi.jl using

trixi_include(default_example())

See also: trixi-framework/Trixi.jl

boundary_condition_do_nothing = Trixi.BoundaryConditionDoNothing()

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

boundary_condition_periodic = Trixi.BoundaryConditionPeriodic()

A singleton struct indicating periodic boundary conditions.

flux_hll

See FluxHLL.

flux_hlle

See min_max_speed_einfeldt. This is a FluxHLL-type two-wave solver with special estimates of the wave speeds.

flux_lax_friedrichs

See FluxLaxFriedrichs.

AMRCallback(semi, controller [,adaptor=AdaptorAMR(semi)];
            interval,
            adapt_initial_condition=true,
            adapt_initial_condition_only_refine=true,
            dynamic_load_balancing=true)

Performs adaptive mesh refinement (AMR) every interval time steps for a given semidiscretization semi using the chosen controller.

AbstractEquations{NDIMS, NVARS}

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

AbstractMesh{NDIMS}

An abstract supertype of specific mesh types such as TreeMesh or StructuredMesh. The type parameters encode the number of spatial dimensions (NDIMS).

AcousticPerturbationEquations2D(v_mean_global, c_mean_global, rho_mean_global)

Acoustic perturbation equations (APE) in two space dimensions. The equations are given by

\[\begin{aligned} \frac{\partial\mathbf{v'}}{\partial t} + \nabla (\bar{\mathbf{v}}\cdot\mathbf{v'}) + \nabla\left( \frac{\bar{c}^2 \tilde{p}'}{\bar{\rho}} \right) &= 0 \\ \frac{\partial \tilde{p}'}{\partial t} + \nabla\cdot (\bar{\rho} \mathbf{v'} + \bar{\mathbf{v}} \tilde{p}') &= 0. \end{aligned}\]

The bar $\bar{(\cdot)}$ indicates time-averaged quantities. The unknowns of the APE are the perturbed velocities $\mathbf{v'} = (v_1', v_2')^T$ and the scaled perturbed pressure $\tilde{p}' = \frac{p'}{\bar{c}^2}$, where $p'$ denotes the perturbed pressure and the perturbed variables are defined by $\phi' = \phi - \bar{\phi}$.

In addition to the unknowns, Trixi.jl currently stores the mean values in the state vector, i.e. the state vector used internally is given by

\[\mathbf{u} = \begin{pmatrix} v_1' \\ v_2' \\ \tilde{p}' \\ \bar{v}_1 \\ \bar{v}_2 \\ \bar{c} \\ \bar{\rho} \end{pmatrix}.\]

This affects the implementation and use of these equations in various ways:

• The flux values corresponding to the mean values must be zero.
• The mean values have to be considered when defining initial conditions, boundary conditions or source terms.
• AnalysisCallback analyzes these variables too.
• Trixi.jl's visualization tools will visualize the mean values by default.

The constructor accepts a 2-tuple v_mean_global and scalars c_mean_global and rho_mean_global which can be used to make the definition of initial conditions for problems with constant mean flow more flexible. These values are ignored if the mean values are defined internally in an initial condition.

The equations are based on the APE-4 system introduced in the following paper:

• Roland Ewert and Wolfgang Schröder (2003) Acoustic perturbation equations based on flow decomposition via source filtering DOI: 10.1016/S0021-9991(03)00168-2

struct Adiabatic

Used to create a no-slip boundary condition with BoundaryConditionNavierStokesWall. The field boundary_value_normal_flux_function should be a function with signature boundary_value_normal_flux_function(x, t, equations) and return a scalar value for the normal heat flux at point x and time t.

AliveCallback(analysis_interval=0, alive_interval=analysis_interval÷10)

Inexpensive callback showing that a simulation is still running by printing some information such as the current time to the screen every alive_interval time steps. If analysis_interval ≂̸ 0, the output is omitted every analysis_interval time steps.

AnalysisCallback(semi; interval=0,
                        save_analysis=false,
                        output_directory="out",
                        analysis_filename="analysis.dat",
                        extra_analysis_errors=Symbol[],
                        extra_analysis_integrals=())

Analyze a numerical solution every interval time steps and print the results to the screen. If save_analysis, the results are also saved in joinpath(output_directory, analysis_filename).

Additional errors can be computed, e.g. by passing extra_analysis_errors = (:l2_error_primitive, :linf_error_primitive) or extra_analysis_errors = (:conservation_error,).

If you want to omit the computation (to safe compute-time) of the default_analysis_errors, specify analysis_errors = Symbol[]. Note: default_analysis_errors are :l2_error and :linf_error for all equations. If you want to compute extra_analysis_errors such as :conservation_error solely, i.e., without :l2_error, :linf_error you need to specify analysis_errors = [:conservation_error] instead of 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 entropy, energy_kinetic, energy_internal, and energy_total. You can also write your own function with the same signature as the examples listed above and pass it via extra_analysis_integrals. See the developer comments about Trixi.analyze, Trixi.pretty_form_utf, and Trixi.pretty_form_ascii for further information on how to create custom analysis quantities.

In addition, the analysis callback records and outputs a number of quantities that are useful for evaluating the computational performance, such as the total runtime, the performance index (time/DOF/rhs!), the time spent in garbage collection (GC), or the current memory usage (alloc'd memory).

AnalysisCallbackCoupled(semi, callbacks...)

Combine multiple analysis callbacks for coupled simulations with a SemidiscretizationCoupled. For each coupled system, an indididual AnalysisCallback must be created and passed to the AnalysisCallbackCoupled in order, i.e., in the same sequence as the indidvidual semidiscretizations are stored in the SemidiscretizationCoupled.

AnalysisSurfaceIntegral{Variable, NBoundaries}(semi,
                                               boundary_symbols::NTuple{NBoundaries, Symbol},
                                               variable)

This struct is used to compute the surface integral of a quantity of interest variable alongside the boundary/boundaries associated with particular name(s) given in boundary_symbol or boundary_symbols. For instance, this can be used to compute the lift LiftCoefficientPressure or drag coefficient DragCoefficientPressure of e.g. an airfoil with the boundary name :Airfoil in 2D.

• boundary_symbols::NTuple{NBoundaries, Symbol}: Name(s) of the boundary/boundaries where the quantity of interest is computed
• variable::Variable: Quantity of interest, like lift or drag

AveragingCallback(semi::SemidiscretizationHyperbolic, tspan; output_directory="out",
                  filename="averaging.h5")

A callback that averages the flow field described by semi which must be a semidiscretization of the compressible Euler equations in two dimensions. The callback records the mean velocity, mean speed of sound, mean density, and mean vorticity for each node over the time interval given by tspan and stores the results in an HDF5 file filename in the directory output_directory. Note that this callback does not support adaptive mesh refinement (AMRCallback).

BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter)

Boundary condition to glue two meshes together. Solution values at the boundary of another mesh will be used as boundary values. This requires the use of SemidiscretizationCoupled. The other mesh is specified by other_semi_index, which is the index of the mesh in the tuple of semidiscretizations.

Note that the elements and nodes of the two meshes at the coupled boundary must coincide. This is currently only implemented for StructuredMesh.

Arguments

• other_semi_index: the index in SemidiscretizationCoupled of the semidiscretization from which the values are copied
• indices::Tuple: node/cell indices at the boundary of the mesh in the other semidiscretization. See examples below.
• uEltype::Type: element type of solution
• coupling_converter::CouplingConverter: function to call for converting the solution state of one system to the other system

Examples

# Connect the left boundary of mesh 2 to our boundary such that our positive
# boundary direction will match the positive y direction of the other boundary
BoundaryConditionCoupled(2, (:begin, :i), Float64, fun)

# Connect the same two boundaries oppositely oriented
BoundaryConditionCoupled(2, (:begin, :i_backwards), Float64, fun)

# Using this as y_neg boundary will connect `our_cells[i, 1, j]` to `other_cells[j, end-i, end]`
BoundaryConditionCoupled(2, (:j, :i_backwards, :end), Float64, fun)

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)

struct BoundaryConditionNavierStokesWall

Creates a wall-type boundary conditions for the compressible Navier-Stokes equations. The fields boundary_condition_velocity and boundary_condition_heat_flux are intended to be boundary condition types such as the NoSlip velocity boundary condition and the Adiabatic or Isothermal heat boundary condition.

BoundaryConditionNeumann(boundary_normal_flux_function)

Similar to BoundaryConditionDirichlet, but creates a Neumann boundary condition for parabolic equations that uses the function boundary_normal_flux_function to specify the values of the normal flux at the boundary. The passed boundary value function will be called with the same arguments as an initial condition function is called, i.e., as

boundary_normal_flux_function(x, t, equations)

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

BoundsCheckCallback(; output_directory="out", save_errors=false, interval=1)

Subcell limiting techniques with SubcellLimiterIDP are constructed to adhere certain local or global bounds. To make sure that these bounds are actually met, this callback calculates the maximum deviation from the bounds. The maximum deviation per applied bound is printed to the screen at the end of the simulation. For more insights, when setting save_errors=true the occurring errors are exported every interval time steps during the simulation. Then, the maximum deviations since the last export are saved in "output_directory/deviations.txt". The BoundsCheckCallback has to be applied as a stage callback for the SSPRK time integration scheme.

  CarpenterKennedy2N43()

Carpenter, Kennedy (1994) Third order 2N storage RK schemes with error control

CarpenterKennedy2N54()

The following structures and methods provide a minimal implementation of the low-storage explicit Runge-Kutta method of

Carpenter, Kennedy (1994) Fourth order 2N storage RK schemes, Solution 3

using the same interface as OrdinaryDiffEq.jl.

CompressibleEulerEquations1D(gamma)

The compressible Euler equations

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho v_1 \\ \rho e \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + p \\ (\rho e +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$ the specific total energy rather than specific internal energy, and

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

the pressure.

CompressibleEulerEquations2D(gamma)

The compressible Euler equations

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho v_1 \\ \rho v_2 \\ \rho e \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ (\rho e +p) v_1 \end{pmatrix} + \frac{\partial}{\partial y} \begin{pmatrix} \rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ (\rho e +p) v_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}\]

for an ideal gas with ratio of specific heats gamma in two space dimensions. Here, $\rho$ is the density, $v_1$, $v_2$ the velocities, $e$ the specific total energy rather than specific internal energy, and

\[p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \right)\]

the pressure.

CompressibleEulerEquations3D(gamma)

The compressible Euler equations

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho v_1 \\ \rho v_2 \\ \rho v_3 \\ \rho e \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ \rho v_1 v_3 \\ ( \rho e +p) v_1 \end{pmatrix} + \frac{\partial}{\partial y} \begin{pmatrix} \rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ \rho v_1 v_3 \\ ( \rho e +p) v_2 \end{pmatrix} + \frac{\partial}{\partial z} \begin{pmatrix} \rho v_3 \\ \rho v_1 v_3 \\ \rho v_2 v_3 \\ \rho v_3^2 + p \\ ( \rho e +p) v_3 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}\]

for an ideal gas with ratio of specific heats gamma in three space dimensions. Here, $\rho$ is the density, $v_1$, $v_2$, $v_3$ the velocities, $e$ the specific total energy rather than specific internal energy, and

\[p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2+v_3^2) \right)\]

the pressure.

CompressibleEulerEquationsQuasi1D(gamma)

The quasi-1d compressible Euler equations (see Chan et al. DOI: 10.48550/arXiv.2307.12089 for details)

\[\frac{\partial}{\partial t} \begin{pmatrix} a \rho \\ a \rho v_1 \\ a e \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} a \rho v_1 \\ a \rho v_1^2 \\ a v_1 (e +p) \end{pmatrix} + a \frac{\partial}{\partial x} \begin{pmatrix} 0 \\ p \\ 0 \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$ the specific total energy rather than specific internal energy, $a$ the (possibly) variable nozzle width, and

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

the pressure.

The nozzle width function $a(x)$ is set inside the initial condition routine for a particular problem setup. To test the conservative form of the compressible Euler equations one can set the nozzle width variable $a$ to one.

In addition to the unknowns, Trixi.jl currently stores the nozzle width values at the approximation points despite being fixed in time. This affects the implementation and use of these equations in various ways:

• The flux values corresponding to the nozzle width must be zero.
• The nozzle width values must be included when defining initial conditions, boundary conditions or source terms.
• AnalysisCallback analyzes this variable.
• Trixi.jl's visualization tools will visualize the nozzle width by default.

CompressibleEulerMulticomponentEquations1D(; gammas, gas_constants)

Multicomponent version of the compressible Euler equations

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho v_1 \\ \rho e \\ \rho_1 \\ \rho_2 \\ \vdots \\ \rho_{n} \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho v_1^2 + p \\ (\rho e +p) v_1 \\ \rho_1 v_1 \\ \rho_2 v_1 \\ \vdots \\ \rho_{n} v_1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix}\]

for calorically perfect gas in one space dimension. Here, $\rho_i$ is the density of component $i$, $\rho=\sum_{i=1}^n\rho_i$ the sum of the individual $\rho_i$, $v_1$ the velocity, $e$ the specific total energy rather than specific internal energy, and

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

the pressure,

\[\gamma=\frac{\sum_{i=1}^n\rho_i C_{v,i}\gamma_i}{\sum_{i=1}^n\rho_i C_{v,i}}\]

total heat capacity ratio, $\gamma_i$ heat capacity ratio of component $i$,

\[C_{v,i}=\frac{R}{\gamma_i-1}\]

specific heat capacity at constant volume of component $i$.

In case of more than one component, the specific heat ratios gammas and the gas constants gas_constants should be passed as tuples, e.g., gammas=(1.4, 1.667).

The remaining variables like the specific heats at constant volume cv or the specific heats at constant pressure cp are then calculated considering a calorically perfect gas.

CompressibleEulerMulticomponentEquations2D(; gammas, gas_constants)

Multicomponent version of the compressible Euler equations

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho v_1 \\ \rho v_2 \\ \rho e \\ \rho_1 \\ \rho_2 \\ \vdots \\ \rho_{n} \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho v_1^2 + p \\ \rho v_1 v_2 \\ ( \rho e +p) v_1 \\ \rho_1 v_1 \\ \rho_2 v_1 \\ \vdots \\ \rho_{n} v_1 \end{pmatrix} + \frac{\partial}{\partial y} \begin{pmatrix} \rho v_1 v_2 \\ \rho v_2^2 + p \\ ( \rho e +p) v_2 \\ \rho_1 v_2 \\ \rho_2 v_2 \\ \vdots \\ \rho_{n} v_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix}\]

for calorically perfect gas in two space dimensions. Here, $\rho_i$ is the density of component $i$, $\rho=\sum_{i=1}^n\rho_i$ the sum of the individual $\rho_i$, $v_1$, $v_2$ the velocities, $e$ the specific total energy rather than specific internal energy, and

\[p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2 + v_2^2) \right)\]

the pressure,

\[\gamma=\frac{\sum_{i=1}^n\rho_i C_{v,i}\gamma_i}{\sum_{i=1}^n\rho_i C_{v,i}}\]

total heat capacity ratio, $\gamma_i$ heat capacity ratio of component $i$,

\[C_{v,i}=\frac{R}{\gamma_i-1}\]

specific heat capacity at constant volume of component $i$.

In case of more than one component, the specific heat ratios gammas and the gas constants gas_constants in [kJ/(kg*K)] should be passed as tuples, e.g., gammas=(1.4, 1.667).

The remaining variables like the specific heats at constant volume cv or the specific heats at constant pressure cp are then calculated considering a calorically perfect gas.

CompressibleNavierStokesDiffusion1D(equations; mu, Pr,
                                    gradient_variables=GradientVariablesPrimitive())

Contains the diffusion (i.e. parabolic) terms applied to mass, momenta, and total energy together with the advective terms from the CompressibleEulerEquations1D.

• equations: instance of the CompressibleEulerEquations1D
• mu: dynamic viscosity,
• Pr: Prandtl number,
• gradient_variables: which variables the gradients are taken with respect to. Defaults to GradientVariablesPrimitive().

Fluid properties such as the dynamic viscosity $\mu$ can be provided in any consistent unit system, e.g., [$\mu$] = kg m⁻¹ s⁻¹. The viscosity $\mu$ may be a constant or a function of the current state, e.g., depending on temperature (Sutherland's law): $\mu = \mu(T)$. In the latter case, the function mu needs to have the signature mu(u, equations).

The particular form of the compressible Navier-Stokes implemented is

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho v \\ \rho e \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho v \\ \rho v^2 + p \\ (\rho e + p) v \end{pmatrix} = \frac{\partial}{\partial x} \begin{pmatrix} 0 \\ \tau \\ \tau v - q \end{pmatrix}\]

where the system is closed with the ideal gas assumption giving

\[p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho v^2 \right)\]

as the pressure. The value of the adiabatic constant gamma is taken from the CompressibleEulerEquations1D. The terms on the right hand side of the system above are built from the viscous stress

\[\tau = \mu \frac{\partial}{\partial x} v\]

where the heat flux is

\[q = -\kappa \frac{\partial}{\partial x} \left(T\right),\quad T = \frac{p}{R\rho}\]

where $T$ is the temperature and $\kappa$ is the thermal conductivity for Fick's law. Under the assumption that the gas has a constant Prandtl number, the thermal conductivity is

\[\kappa = \frac{\gamma \mu R}{(\gamma - 1)\textrm{Pr}}.\]

From this combination of temperature $T$ and thermal conductivity $\kappa$ we see that the gas constant R cancels and the heat flux becomes

\[q = -\kappa \frac{\partial}{\partial x} \left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}} \frac{\partial}{\partial x} \left(\frac{p}{\rho}\right)\]

which is the form implemented below in the flux function.

In one spatial dimensions we require gradients for two quantities, e.g., primitive quantities

\[\frac{\partial}{\partial x} v,\, \frac{\partial}{\partial x} T\]

or the entropy variables

\[\frac{\partial}{\partial x} w_2,\, \frac{\partial}{\partial x} w_3\]

where

\[w_2 = \frac{\rho v1}{p},\, w_3 = -\frac{\rho}{p}\]

CompressibleNavierStokesDiffusion2D(equations; mu, Pr,
                                    gradient_variables=GradientVariablesPrimitive())

Contains the diffusion (i.e. parabolic) terms applied to mass, momenta, and total energy together with the advective terms from the CompressibleEulerEquations2D.

• equations: instance of the CompressibleEulerEquations2D
• mu: dynamic viscosity,
• Pr: Prandtl number,
• gradient_variables: which variables the gradients are taken with respect to. Defaults to GradientVariablesPrimitive().

Fluid properties such as the dynamic viscosity $\mu$ can be provided in any consistent unit system, e.g., [$\mu$] = kg m⁻¹ s⁻¹. The viscosity $\mu$ may be a constant or a function of the current state, e.g., depending on temperature (Sutherland's law): $\mu = \mu(T)$. In the latter case, the function mu needs to have the signature mu(u, equations).

The particular form of the compressible Navier-Stokes implemented is

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho \mathbf{v} \\ \rho e \end{pmatrix} + \nabla \cdot \begin{pmatrix} \rho \mathbf{v} \\ \rho \mathbf{v}\mathbf{v}^T + p \underline{I} \\ (\rho e + p) \mathbf{v} \end{pmatrix} = \nabla \cdot \begin{pmatrix} 0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \mathbf{q} \end{pmatrix}\]

where the system is closed with the ideal gas assumption giving

\[p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \right)\]

as the pressure. The value of the adiabatic constant gamma is taken from the CompressibleEulerEquations2D. The terms on the right hand side of the system above are built from the viscous stress tensor

\[\underline{\tau} = \mu \left(\nabla\mathbf{v} + \left(\nabla\mathbf{v}\right)^T\right) - \frac{2}{3} \mu \left(\nabla\cdot\mathbf{v}\right)\underline{I}\]

where $\underline{I}$ is the $2\times 2$ identity matrix and the heat flux is

\[\mathbf{q} = -\kappa\nabla\left(T\right),\quad T = \frac{p}{R\rho}\]

where $T$ is the temperature and $\kappa$ is the thermal conductivity for Fick's law. Under the assumption that the gas has a constant Prandtl number, the thermal conductivity is

\[\kappa = \frac{\gamma \mu R}{(\gamma - 1)\textrm{Pr}}.\]

From this combination of temperature $T$ and thermal conductivity $\kappa$ we see that the gas constant R cancels and the heat flux becomes

\[\mathbf{q} = -\kappa\nabla\left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}}\nabla\left(\frac{p}{\rho}\right)\]

which is the form implemented below in the flux function.

In two spatial dimensions we require gradients for three quantities, e.g., primitive quantities

\[\nabla v_1,\, \nabla v_2,\, \nabla T\]

or the entropy variables

\[\nabla w_2,\, \nabla w_3,\, \nabla w_4\]

where

\[w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p}\]

CompressibleNavierStokesDiffusion3D(equations; mu, Pr,
                                    gradient_variables=GradientVariablesPrimitive())

Contains the diffusion (i.e. parabolic) terms applied to mass, momenta, and total energy together with the advective terms from the CompressibleEulerEquations3D.

• equations: instance of the CompressibleEulerEquations3D
• mu: dynamic viscosity,
• Pr: Prandtl number,
• gradient_variables: which variables the gradients are taken with respect to. Defaults to GradientVariablesPrimitive().

Fluid properties such as the dynamic viscosity $\mu$ can be provided in any consistent unit system, e.g., [$\mu$] = kg m⁻¹ s⁻¹. The viscosity $\mu$ may be a constant or a function of the current state, e.g., depending on temperature (Sutherland's law): $\mu = \mu(T)$. In the latter case, the function mu needs to have the signature mu(u, equations).

The particular form of the compressible Navier-Stokes implemented is

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho \mathbf{v} \\ \rho e \end{pmatrix} + \nabla \cdot \begin{pmatrix} \rho \mathbf{v} \\ \rho \mathbf{v}\mathbf{v}^T + p \underline{I} \\ (\rho e + p) \mathbf{v} \end{pmatrix} = \nabla \cdot \begin{pmatrix} 0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \mathbf{q} \end{pmatrix}\]

where the system is closed with the ideal gas assumption giving

\[p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2+v_3^2) \right)\]

as the pressure. The value of the adiabatic constant gamma is taken from the CompressibleEulerEquations2D. The terms on the right hand side of the system above are built from the viscous stress tensor

\[\underline{\tau} = \mu \left(\nabla\mathbf{v} + \left(\nabla\mathbf{v}\right)^T\right) - \frac{2}{3} \mu \left(\nabla\cdot\mathbf{v}\right)\underline{I}\]

where $\underline{I}$ is the $3\times 3$ identity matrix and the heat flux is

\[\mathbf{q} = -\kappa\nabla\left(T\right),\quad T = \frac{p}{R\rho}\]

where $T$ is the temperature and $\kappa$ is the thermal conductivity for Fick's law. Under the assumption that the gas has a constant Prandtl number, the thermal conductivity is

\[\kappa = \frac{\gamma \mu R}{(\gamma - 1)\textrm{Pr}}.\]

From this combination of temperature $T$ and thermal conductivity $\kappa$ we see that the gas constant R cancels and the heat flux becomes

\[\mathbf{q} = -\kappa\nabla\left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}}\nabla\left(\frac{p}{\rho}\right)\]

which is the form implemented below in the flux function.

In two spatial dimensions we require gradients for three quantities, e.g., primitive quantities

\[\nabla v_1,\, \nabla v_2,\, \nabla v_3,\, \nabla T\]

or the entropy variables

\[\nabla w_2,\, \nabla w_3,\, \nabla w_4\, \nabla w_5\]

where

\[w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p},\, w_5 = -\frac{\rho}{p}\]

ControllerThreeLevel(semi, indicator; base_level=1,
                                      med_level=base_level, med_threshold=0.0,
                                      max_level=base_level, max_threshold=1.0)

An AMR controller based on three levels (in descending order of precedence):

• set the target level to max_level if indicator > max_threshold
• set the target level to med_level if indicator > med_threshold; if med_level < 0, set the target level to the current level
• set the target level to base_level otherwise

ControllerThreeLevelCombined(semi, indicator_primary, indicator_secondary;
                             base_level=1,
                             med_level=base_level, med_threshold=0.0,
                             max_level=base_level, max_threshold=1.0,
                             max_threshold_secondary=1.0)

An AMR controller based on three levels (in descending order of precedence):

• set the target level to max_level if indicator_primary > max_threshold
• set the target level to med_level if indicator_primary > med_threshold; if med_level < 0, set the target level to the current level
• set the target level to base_level otherwise

If indicator_secondary >= max_threshold_secondary, set the target level to max_level.

DG(; basis, mortar, surface_integral, volume_integral)

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

DGMulti(approximation_type::AbstractDerivativeOperator;
        element_type::AbstractElemShape,
        surface_flux=flux_central,
        surface_integral=SurfaceIntegralWeakForm(surface_flux),
        volume_integral=VolumeIntegralWeakForm(),
        kwargs...)

Create a summation by parts (SBP) discretization on the given element_type using a tensor product structure based on the 1D SBP derivative operator passed as approximation_type.

For more info, see the documentations of StartUpDG.jl and SummationByPartsOperators.jl.

DGMulti(; polydeg::Integer,
          element_type::AbstractElemShape,
          approximation_type=Polynomial(),
          surface_flux=flux_central,
          surface_integral=SurfaceIntegralWeakForm(surface_flux),
          volume_integral=VolumeIntegralWeakForm(),
          RefElemData_kwargs...)

Create a discontinuous Galerkin method which uses

• approximations of polynomial degree polydeg
• element type element_type (Tri(), Quad(), Tet(), and Hex() currently supported)

Optional:

• approximation_type (default is Polynomial(); SBP() also supported for Tri(), Quad(), and Hex() element types).
• RefElemData_kwargs are additional keyword arguments for RefElemData, such as quad_rule_vol. For more info, see the StartUpDG.jl docs.

DGMultiMesh{NDIMS, ...}

DGMultiMesh describes a mesh type which wraps StartUpDG.MeshData and boundary_faces in a dispatchable type. This is intended to store geometric data and connectivities for any type of mesh (Cartesian, affine, curved, structured/unstructured).

DGMultiMesh(dg::DGMulti{2, Tri}, triangulateIO, boundary_dict::Dict{Symbol, Int})

• dg::DGMulti contains information associated with to the reference element (e.g., quadrature, basis evaluation, differentiation, etc).
• triangulateIO is a TriangulateIO mesh representation
• boundary_dict is a Dict{Symbol, Int} which associates each integer TriangulateIO boundary tag with a Symbol.

DGMultiMesh(dg::DGMulti)

Constructs a single-element DGMultiMesh for a single periodic element given a DGMulti with approximation_type set to a periodic (finite difference) SBP operator from SummationByPartsOperators.jl.

DGMultiMesh(dg::DGMulti{NDIMS}, vertex_coordinates, EToV;
            is_on_boundary=nothing,
            periodicity=ntuple(_->false, NDIMS)) where {NDIMS}

• dg::DGMulti contains information associated with to the reference element (e.g., quadrature, basis evaluation, differentiation, etc).
• vertex_coordinates is a tuple of vectors containing x,y,... components of the vertex coordinates
• EToV is a 2D array containing element-to-vertex connectivities for each element
• is_on_boundary specifies boundary using a Dict{Symbol, <:Function}
• periodicity is a tuple of booleans specifying if the domain is periodic true/false in the (x,y,z) direction.

DGMultiMesh(dg::DGMulti{NDIMS}, cells_per_dimension, mapping;
            is_on_boundary=nothing,
            periodicity=ntuple(_ -> false, NDIMS), kwargs...) where {NDIMS}

Constructs a Curved() DGMultiMesh with element type dg.basis.element_type.

• mapping is a function which maps from a reference [-1, 1]^NDIMS domain to a mapped domain, e.g., xy = mapping(x, y) in 2D.
• is_on_boundary specifies boundary using a Dict{Symbol, <:Function}
• periodicity is a tuple of Bools specifying periodicity = true/false in the (x,y,z) direction.

DGMultiMesh(dg::DGMulti, cells_per_dimension;
            coordinates_min=(-1.0, -1.0), coordinates_max=(1.0, 1.0),
            is_on_boundary=nothing,
            periodicity=ntuple(_ -> false, NDIMS))

Constructs a Cartesian DGMultiMesh with element type dg.basis.element_type. The domain is the tensor product of the intervals [coordinates_min[i], coordinates_max[i]].

• is_on_boundary specifies boundary using a Dict{Symbol, <:Function}
• periodicity is a tuple of Bools specifying periodicity = true/false in the (x,y,z) direction.

DGMultiMesh(dg::DGMulti, filename::String)

• dg::DGMulti contains information associated with the reference element (e.g., quadrature, basis evaluation, differentiation, etc).
• filename is a path specifying a .mesh file generated by HOHQMesh.

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

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

DissipationGlobalLaxFriedrichs(λ)

Create a global Lax-Friedrichs dissipation operator with dissipation coefficient λ.

DissipationLocalLaxFriedrichs(max_abs_speed=max_abs_speed_naive)

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

DragCoefficientPressure(aoa, rhoinf, uinf, linf)

Compute the drag coefficient

\[C_{D,p} \coloneqq \frac{\oint_{\partial \Omega} p \boldsymbol n \cdot \psi_D \, \mathrm{d} S} {0.5 \rho_{\infty} U_{\infty}^2 L_{\infty}}\]

based on the pressure distribution along a boundary. Supposed to be used in conjunction with AnalysisSurfaceIntegral which stores the boundary information and semidiscretization.

• aoa::Real: Angle of attack in radians (for airfoils etc.)
• rhoinf::Real: Free-stream density
• uinf::Real: Free-stream velocity
• linf::Real: Reference length of geometry (e.g. airfoil chord length)

DragCoefficientShearStress(aoa, rhoinf, uinf, linf)

Compute the drag coefficient

\[C_{D,f} \coloneqq \frac{\oint_{\partial \Omega} \boldsymbol \tau_w \cdot \psi_D \, \mathrm{d} S} {0.5 \rho_{\infty} U_{\infty}^2 L_{\infty}}\]

based on the wall shear stress vector $\tau_w$ along a boundary. Supposed to be used in conjunction with AnalysisSurfaceIntegral which stores the boundary information and semidiscretization.

• aoa::Real: Angle of attack in radians (for airfoils etc.)
• rhoinf::Real: Free-stream density
• uinf::Real: Free-stream velocity
• linf::Real: Reference length of geometry (e.g. airfoil chord length)

EntropyBoundedLimiter{RealT <: Real}(; exp_entropy_decrease_max = 1f-13)

Entropy-bounded limiter by

• Lv, Ihme (2015) Entropy-bounded discontinuous Galerkin scheme for Euler equations doi: 10.1016/j.jcp.2015.04.026

This is an ideal-gas specific limiter that bounds the (unphysical) decrease of the thermodynamic entropy per element from one time step (or Runge-Kutta stage) to the next. The parameter exp_entropy_decrease_max is the maximum allowed exponentiated entropy decrease per element at each element's node.

In the original version of the paper, this value is set to zero to ensure that the entropy does not decrease, i.e., guarantee entropy stability in the sense of

• Tadmor (1986) A minimum entropy principle in the gas dynamics equations doi: 10.1016/0168-9274(86)90029-2

This, however, leads in general to very diffusive solutions for timesteps violating a CFL condition (Lemma 3 in Lv, Ihme (2015)) which is required for entropy stability in the mean values. Since most practical simulations will employ a significantly larger timestep, one can relax the strict entropy increase requirement by setting exp_entropy_decrease_max to a negative value. The limiter acts if the exponentiated entropy decrease on an element is larger than exp_entropy_decrease_max. This means that if the change in exponentiated entropy lies below exp_entropy_decrease_max (i.e., larger in absolute value) the limiter takes action. The choice of the tolerated exponentiated entropy decrease is a problem-specific parameter which balances the trade-off between accuracy and stability.

EulerAcousticsCouplingCallback

A callback that couples the acoustic perturbation equations and compressible Euler equations. Must be used in conjunction with SemidiscretizationEulerAcoustics. This callback manages the flow solver - which is always one time step ahead of the acoustics solver - and calculates the acoustic source term after each time step. The linearized Lamb vector is used as the source term, i.e.

\[\mathbf{s} = -(\mathbf{\omega'} \times \bar{\mathbf{v}} + \bar{\mathbf{\omega}} \times \mathbf{v'}),\]

where $\mathbf{v}$ denotes the velocity, $\mathbf{\omega}$ denotes the vorticity, the bar $\bar{(\cdot)}$ indicates time-averaged quantities (see AveragingCallback) and prime $(\cdot)'$ denotes perturbed quantities defined by $\phi' = \phi - \bar{\phi}$. Note that the perturbed quantities here are based entirely on the pure flow solution and should not be confused with the state variables of the acoustic perturbation equations.

In addition, this callback manages the time step size for both solvers and initializes the mean values of the acoustic perturbation equations using results obtained with the AveragingCallback.

• Michael Schlottke-Lakemper (2017) A direct-hybrid method for aeroacoustic analysis DOI: 10.18154/RWTH-2017-04082

EulerAcousticsCouplingCallback(ode_euler, averaging_file::AbstractString, alg,
                               cfl_acoustics::Real, cfl_euler::Real; kwargs...)

Creates an EulerAcousticsCouplingCallback based on the pure flow ODEProblem given by ode_euler. Creates an integrator using the time integration method alg and the keyword arguments to solve ode_euler (consult the OrdinaryDiffEq documentation for further information). Manages the step size for both solvers by using the minimum of the maximum step size obtained with CFL numbers cfl_acoustics for the acoustics solver and cfl_euler for and flow solver, respectively. The mean values for the acoustic perturbation equations are read from averaging_file (see AveragingCallback).

EulerAcousticsCouplingCallback(ode_euler,
                               averaging_callback::DiscreteCallback{<:Any, <:AveragingCallback},
                               alg, cfl_acoustics::Real, cfl_euler::Real; kwargs...)

Creates an EulerAcousticsCouplingCallback based on the pure flow ODEProblem given by ode_euler. Creates an integrator using the time integration method alg and the keyword arguments to solve ode_euler (consult the OrdinaryDiffEq documentation for further information). Manages the step size for both solvers by using the minimum of the maximum step size obtained with CFL numbers cfl_acoustics for the acoustics solver and cfl_euler for and flow solver, respectively. The mean values for the acoustic perturbation equations are read from averaging_callback (see AveragingCallback).

FDSBP(D_SBP; surface_integral, volume_integral)

Specialization of DG methods that uses general summation by parts (SBP) operators from SummationByPartsOperators.jl. In particular, this includes classical finite difference (FD) SBP methods. These methods have the same structure as classical DG methods - local operations on elements with connectivity through interfaces without imposing any continuity constraints.

D_SBP is an SBP derivative operator from SummationByPartsOperators.jl. The other arguments have the same meaning as in DG or DGSEM.

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

FluxHydrostaticReconstruction(numerical_flux, hydrostatic_reconstruction)

Allow for some kind of hydrostatic reconstruction of the solution state prior to the surface flux computation. This is a particular strategy to ensure that the method remains well-balanced for the shallow water equations, see ShallowWaterEquations1D or ShallowWaterEquations2D.

For example, the hydrostatic reconstruction from Audusse et al. is implemented in one and two spatial dimensions, see hydrostatic_reconstruction_audusse_etal or the original paper

• Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoit Perthame (2004) A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows DOI: 10.1137/S1064827503431090

Other hydrostatic reconstruction techniques are available, particularly to handle wet / dry fronts. A good overview of the development and application of hydrostatic reconstruction can be found in

• Guoxian Chen and Sebastian Noelle A unified surface-gradient and hydrostatic reconstruction scheme for the shallow water equations (2021) RWTH Aachen preprint
• Andreas Buttinger-Kreuzhuber, Zsolt Horváth, Sebastian Noelle, Günter Blöschl and Jürgen Waser (2019) A fast second-order shallow water scheme on two-dimensional structured grids over abrupt topography DOI: 10.1016/j.advwatres.2019.03.010

FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction,
             equations::CompressibleEulerEquations2D)

Low Mach number approximate Riemann solver (LMARS) for atmospheric flows using an estimate c of the speed of sound.

References:

• Xi Chen et al. (2013) A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian Coordinate DOI: 10.1175/MWR-D-12-00129.1

FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction,
             equations::CompressibleEulerEquations3D)

Low Mach number approximate Riemann solver (LMARS) for atmospheric flows using an estimate c of the speed of sound.

References:

• Xi Chen et al. (2013) A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian Coordinate DOI: 10.1175/MWR-D-12-00129.1

FluxLaxFriedrichs(max_abs_speed=max_abs_speed_naive)

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

FluxPlusDissipation(numerical_flux, dissipation)

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

FluxRotated(numerical_flux)

Compute a numerical_flux flux in direction of a normal vector by rotating the solution, computing the numerical flux in x-direction, and rotating the calculated flux back.

Requires a rotationally invariant equation with equation-specific functions rotate_to_x and rotate_from_x.

FluxUpwind(splitting)

A numerical flux f(u_left, u_right) = f⁺(u_left) + f⁻(u_right) based on flux vector splitting.

The SurfaceIntegralUpwind with a given splitting is equivalent to the SurfaceIntegralStrongForm with FluxUpwind(splitting) as numerical flux (up to floating point differences). Note, that SurfaceIntegralUpwind is only available on TreeMesh.

GlmSpeedCallback(; glm_scale=0.5, cfl, semi_indices=())

Update the divergence cleaning wave speed c_h according to the time step computed in StepsizeCallback for the ideal GLM-MHD equations. The cfl number should be set to the same value as for the time step size calculation. The glm_scale ensures that the GLM wave speed is lower than the fastest physical waves in the MHD solution and should thus be set to a value within the interval [0,1]. Note that glm_scale = 0 deactivates the divergence cleaning.

In case of coupled semidiscretizations, specify for which semi_index, i.e. index of the semidiscretization, the divergence cleaning should be applied. See also SemidiscretizationCoupled. Note: SemidiscretizationCoupled and all related features are considered experimental and may change at any time.

GradientVariablesPrimitive and GradientVariablesEntropy are gradient variable type parameters for CompressibleNavierStokesDiffusion1D. By default, the gradient variables are set to be GradientVariablesPrimitive. Specifying GradientVariablesEntropy instead uses the entropy variable formulation from

• Hughes, Mallet, Franca (1986) A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the compressible Euler and Navier-Stokes equations and the second law of thermodynamics. https://doi.org/10.1016/0045-7825(86)90127-1

Under GradientVariablesEntropy, the Navier-Stokes discretization is provably entropy stable.

HypDiffN3Erk3Sstar52()

Five stage, second-order accurate explicit Runge-Kutta scheme with stability region optimized for the hyperbolic diffusion equation with LLF flux and polynomials of degree polydeg=3.

HyperbolicDiffusionEquations1D

The linear hyperbolic diffusion equations in one space dimension. A description of this system can be found in Sec. 2.5 of the book

• Masatsuka (2013) I Do Like CFD, Too: Vol 1. Freely available at http://www.cfdbooks.com/

Further analysis can be found in the paper

• Nishikawa (2007) A first-order system approach for diffusion equation. I: Second-order residual-distribution schemes DOI: 10.1016/j.jcp.2007.07.029

HyperbolicDiffusionEquations2D

The linear hyperbolic diffusion equations in two space dimensions. A description of this system can be found in Sec. 2.5 of the book "I Do Like CFD, Too: Vol 1". The book is freely available at http://www.cfdbooks.com/ and further analysis can be found in the paper by Nishikawa DOI: 10.1016/j.jcp.2007.07.029

HyperbolicDiffusionEquations3D

The linear hyperbolic diffusion equations in three space dimensions. A description of this system can be found in Sec. 2.5 of the book "I Do Like CFD, Too: Vol 1". The book is freely available at http://www.cfdbooks.com/ and further analysis can be found in the paper by Nishikawa DOI: 10.1016/j.jcp.2007.07.029

IdealGlmMhdEquations1D(gamma)

The ideal compressible GLM-MHD equations for an ideal gas with ratio of specific heats gamma in one space dimension.

IdealGlmMhdEquations2D(gamma)

The ideal compressible GLM-MHD equations for an ideal gas with ratio of specific heats gamma in two space dimensions.

IdealGlmMhdEquations3D(gamma)

The ideal compressible GLM-MHD equations for an ideal gas with ratio of specific heats gamma in three space dimensions.

IdealGlmMhdMulticomponentEquations1D

The ideal compressible multicomponent GLM-MHD equations in one space dimension.

IdealGlmMhdMulticomponentEquations2D

The ideal compressible multicomponent GLM-MHD equations in two space dimensions.

IndicatorHennemannGassner(equations::AbstractEquations, basis;
                          alpha_max=0.5,
                          alpha_min=0.001,
                          alpha_smooth=true,
                          variable)
IndicatorHennemannGassner(semi::AbstractSemidiscretization;
                          alpha_max=0.5,
                          alpha_min=0.001,
                          alpha_smooth=true,
                          variable)

Indicator used for shock-capturing (when passing the equations and the basis) or adaptive mesh refinement (AMR, when passing the semi).

See also VolumeIntegralShockCapturingHG.

References

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

IndicatorLöhner (equivalent to IndicatorLoehner)

IndicatorLöhner(equations::AbstractEquations, basis;
                f_wave=0.2, variable)
IndicatorLöhner(semi::AbstractSemidiscretization;
                f_wave=0.2, variable)

AMR indicator adapted from a FEM indicator by Löhner (1987), also used in the FLASH code as standard AMR indicator. The indicator estimates a weighted second derivative of a specified variable locally.

When constructed to be used for AMR, pass the semi. Pass the equations, and basis if this indicator should be used for shock capturing.

References

• Löhner (1987) "An adaptive finite element scheme for transient problems in CFD" doi: 10.1016/0045-7825(87)90098-3
• https://flash.rochester.edu/site/flashcode/usersupport/flash4ug_4p62/node59.html#SECTION05163100000000000000

IndicatorMax(equations::AbstractEquations, basis; variable)
IndicatorMax(semi::AbstractSemidiscretization; variable)

A simple indicator returning the maximum of variable in an element. When constructed to be used for AMR, pass the semi. Pass the equations, and basis if this indicator should be used for shock capturing.

InviscidBurgersEquation1D

The inviscid Burgers' equation

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

in one space dimension.

struct Isothermal

Used to create a no-slip boundary condition with BoundaryConditionNavierStokesWall. The field boundary_value_function should be a function with signature boundary_value_function(x, t, equations) and return a scalar value for the temperature at point x and time t.

LaplaceDiffusion1D(diffusivity, equations)

LaplaceDiffusion1D represents a scalar diffusion term $\nabla \cdot (\kappa\nabla u))$ with diffusivity $\kappa$ applied to each solution component defined by equations.

LaplaceDiffusion2D(diffusivity, equations)

LaplaceDiffusion2D represents a scalar diffusion term $\nabla \cdot (\kappa\nabla u))$ with diffusivity $\kappa$ applied to each solution component defined by equations.

LaplaceDiffusion3D(diffusivity, equations)

LaplaceDiffusion3D represents a scalar diffusion term $\nabla \cdot (\kappa\nabla u))$ with diffusivity $\kappa$ applied to each solution component defined by equations.

LatticeBoltzmannEquations2D(; Ma, Re, collision_op=collision_bgk,
                           c=1, L=1, rho0=1, u0=nothing, nu=nothing)

The Lattice-Boltzmann equations

\[\partial_t u_\alpha + v_{\alpha,1} \partial_1 u_\alpha + v_{\alpha,2} \partial_2 u_\alpha = 0\]

in two space dimensions for the D2Q9 scheme.

The characteristic Mach number and Reynolds numbers are specified as Ma and Re. By the default, the collision operator collision_op is set to the BGK model. c, L, and rho0 specify the mean thermal molecular velocity, the characteristic length, and the reference density, respectively. They can usually be left to the default values. If desired, instead of the Mach number, one can set the macroscopic reference velocity u0 directly (Ma needs to be set to nothing in this case). Likewise, instead of the Reynolds number one can specify the kinematic viscosity nu directly (in this case, Re needs to be set to nothing).

The nine discrete velocity directions of the D2Q9 scheme are sorted as follows [4]:

  6     2     5       y
    ┌───┼───┐         │
    │       │         │
  3 ┼   9   ┼ 1        ──── x
    │       │        ╱
    └───┼───┘       ╱
  7     4     8    z

Note that usually the velocities are numbered from 0 to 8, where 0 corresponds to the zero velocity. Due to Julia using 1-based indexing, here we use indices from 1 to 9, where 1 through 8 correspond to the velocity directions in [4] and 9 is the zero velocity.

The corresponding opposite directions are:

• 1 ←→ 3
• 2 ←→ 4
• 3 ←→ 1
• 4 ←→ 2
• 5 ←→ 7
• 6 ←→ 8
• 7 ←→ 5
• 8 ←→ 6
• 9 ←→ 9

The main sources for the base implementation were

1. Misun Min, Taehun Lee, A spectral-element discontinuous Galerkin lattice Boltzmann method for nearly incompressible flows, J Comput Phys 230(1), 2011 doi:10.1016/j.jcp.2010.09.024
2. Karsten Golly, Anwendung der Lattice-Boltzmann Discontinuous Galerkin Spectral Element Method (LB-DGSEM) auf laminare und turbulente nahezu inkompressible Strömungen im dreidimensionalen Raum, Master Thesis, University of Cologne, 2018.
3. Dieter Hänel, Molekulare Gasdynamik, Springer-Verlag Berlin Heidelberg, 2004 doi:10.1007/3-540-35047-0
4. Dieter Krüger et al., The Lattice Boltzmann Method, Springer International Publishing, 2017 doi:10.1007/978-3-319-44649-3

LatticeBoltzmannEquations3D(; Ma, Re, collision_op=collision_bgk,
                           c=1, L=1, rho0=1, u0=nothing, nu=nothing)

The Lattice-Boltzmann equations

\[\partial_t u_\alpha + v_{\alpha,1} \partial_1 u_\alpha + v_{\alpha,2} \partial_2 u_\alpha + v_{\alpha,3} \partial_3 u_\alpha = 0\]

in three space dimensions for the D3Q27 scheme.

The characteristic Mach number and Reynolds numbers are specified as Ma and Re. By the default, the collision operator collision_op is set to the BGK model. c, L, and rho0 specify the mean thermal molecular velocity, the characteristic length, and the reference density, respectively. They can usually be left to the default values. If desired, instead of the Mach number, one can set the macroscopic reference velocity u0 directly (Ma needs to be set to nothing in this case). Likewise, instead of the Reynolds number one can specify the kinematic viscosity nu directly (in this case, Re needs to be set to nothing).

The twenty-seven discrete velocity directions of the D3Q27 scheme are sorted as follows [4]:

• plane at z = -1:

    24    17     21       y
       ┌───┼───┐          │
       │       │          │
    10 ┼   6   ┼ 15        ──── x
       │       │         ╱
       └───┼───┘        ╱
    20    12     26    z

• plane at z = 0:

    14     3     7        y
       ┌───┼───┐          │
       │       │          │
     2 ┼  27   ┼ 1         ──── x
       │       │         ╱
       └───┼───┘        ╱
     8     4     13    z

• plane at z = +1:

    25    11     19       y
       ┌───┼───┐          │
       │       │          │
    16 ┼   5   ┼ 9         ──── x
       │       │         ╱
       └───┼───┘        ╱
    22    18     23    z

Note that usually the velocities are numbered from 0 to 26, where 0 corresponds to the zero velocity. Due to Julia using 1-based indexing, here we use indices from 1 to 27, where 1 through 26 correspond to the velocity directions in [4] and 27 is the zero velocity.

The corresponding opposite directions are:

• 1 ←→ 2
• 2 ←→ 1
• 3 ←→ 4
• 4 ←→ 3
• 5 ←→ 6
• 6 ←→ 5
• 7 ←→ 8
• 8 ←→ 7
• 9 ←→ 10
• 10 ←→ 9
• 11 ←→ 12
• 12 ←→ 11
• 13 ←→ 14
• 14 ←→ 13
• 15 ←→ 16
• 16 ←→ 15
• 17 ←→ 18
• 18 ←→ 17
• 19 ←→ 20
• 20 ←→ 19
• 21 ←→ 22
• 22 ←→ 21
• 23 ←→ 24
• 24 ←→ 23
• 25 ←→ 26
• 26 ←→ 25
• 27 ←→ 27

The main sources for the base implementation were

1. Misun Min, Taehun Lee, A spectral-element discontinuous Galerkin lattice Boltzmann method for nearly incompressible flows, J Comput Phys 230(1), 2011 doi:10.1016/j.jcp.2010.09.024
2. Karsten Golly, Anwendung der Lattice-Boltzmann Discontinuous Galerkin Spectral Element Method (LB-DGSEM) auf laminare und turbulente nahezu inkompressible Strömungen im dreidimensionalen Raum, Master Thesis, University of Cologne, 2018.
3. Dieter Hänel, Molekulare Gasdynamik, Springer-Verlag Berlin Heidelberg, 2004 doi:10.1007/3-540-35047-0
4. Dieter Krüger et al., The Lattice Boltzmann Method, Springer International Publishing, 2017 doi:10.1007/978-3-319-44649-3

LiftCoefficientPressure(aoa, rhoinf, uinf, linf)

Compute the lift coefficient

\[C_{L,p} \coloneqq \frac{\oint_{\partial \Omega} p \boldsymbol n \cdot \psi_L \, \mathrm{d} S} {0.5 \rho_{\infty} U_{\infty}^2 L_{\infty}}\]

based on the pressure distribution along a boundary. Supposed to be used in conjunction with AnalysisSurfaceIntegral which stores the boundary information and semidiscretization.

• aoa::Real: Angle of attack in radians (for airfoils etc.)
• rhoinf::Real: Free-stream density
• uinf::Real: Free-stream velocity
• linf::Real: Reference length of geometry (e.g. airfoil chord length)

LiftCoefficientShearStress(aoa, rhoinf, uinf, linf)

Compute the lift coefficient

\[C_{L,f} \coloneqq \frac{\oint_{\partial \Omega} \boldsymbol \tau_w \cdot \psi_L \, \mathrm{d} S} {0.5 \rho_{\infty} U_{\infty}^2 L_{\infty}}\]

based on the wall shear stress vector $\tau_w$ along a boundary. Supposed to be used in conjunction with AnalysisSurfaceIntegral which stores the boundary information and semidiscretization.

• aoa::Real: Angle of attack in radians (for airfoils etc.)
• rhoinf::Real: Free-stream density
• uinf::Real: Free-stream velocity
• linf::Real: Reference length of geometry (e.g. airfoil chord length)

LinearScalarAdvectionEquation1D

The linear scalar advection equation

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

in one space dimension with constant velocity a.

LinearScalarAdvectionEquation2D

The linear scalar advection equation

\[\partial_t u + a_1 \partial_1 u + a_2 \partial_2 u = 0\]

in two space dimensions with constant velocity a.

LinearScalarAdvectionEquation3D

The linear scalar advection equation

\[\partial_t u + a_1 \partial_1 u + a_2 \partial_2 u + a_3 \partial_3 u = 0\]

in three space dimensions with constant velocity a.

LinearizedEulerEquations1D(v_mean_global, c_mean_global, rho_mean_global)

Linearized Euler equations in one space dimension. The equations are given by

\[\partial_t \begin{pmatrix} \rho' \\ v_1' \\ p' \end{pmatrix} + \partial_x \begin{pmatrix} \bar{\rho} v_1' + \bar{v_1} \rho ' \\ \bar{v_1} v_1' + \frac{p'}{\bar{\rho}} \\ \bar{v_1} p' + c^2 \bar{\rho} v_1' \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}\]

The bar $\bar{(\cdot)}$ indicates uniform mean flow variables and $c$ is the speed of sound. The unknowns are the perturbation quantities of the acoustic velocity $v_1'$, the pressure $p'$ and the density $\rho'$.

LinearizedEulerEquations2D(v_mean_global, c_mean_global, rho_mean_global)

Linearized Euler equations in two space dimensions. The equations are given by

\[\partial_t \begin{pmatrix} \rho' \\ v_1' \\ v_2' \\ p' \end{pmatrix} + \partial_x \begin{pmatrix} \bar{\rho} v_1' + \bar{v_1} \rho ' \\ \bar{v_1} v_1' + \frac{p'}{\bar{\rho}} \\ \bar{v_1} v_2' \\ \bar{v_1} p' + c^2 \bar{\rho} v_1' \end{pmatrix} + \partial_y \begin{pmatrix} \bar{\rho} v_2' + \bar{v_2} \rho ' \\ \bar{v_2} v_1' \\ \bar{v_2} v_2' + \frac{p'}{\bar{\rho}} \\ \bar{v_2} p' + c^2 \bar{\rho} v_2' \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}\]

The bar $\bar{(\cdot)}$ indicates uniform mean flow variables and $c$ is the speed of sound. The unknowns are the acoustic velocities $v' = (v_1', v_2')$, the pressure $p'$ and the density $\rho'$.

LinearizedEulerEquations3D(v_mean_global, c_mean_global, rho_mean_global)

Linearized Euler equations in three space dimensions. The equations are given by

\[\partial_t \begin{pmatrix} \rho' \\ v_1' \\ v_2' \\ v_3' \\ p' \end{pmatrix} + \partial_x \begin{pmatrix} \bar{\rho} v_1' + \bar{v_1} \rho ' \\ \bar{v_1} v_1' + \frac{p'}{\bar{\rho}} \\ \bar{v_1} v_2' \\ \bar{v_1} v_3' \\ \bar{v_1} p' + c^2 \bar{\rho} v_1' \end{pmatrix} + \partial_y \begin{pmatrix} \bar{\rho} v_2' + \bar{v_2} \rho ' \\ \bar{v_2} v_1' \\ \bar{v_2} v_2' + \frac{p'}{\bar{\rho}} \\ \bar{v_2} v_3' \\ \bar{v_2} p' + c^2 \bar{\rho} v_2' \end{pmatrix} + \partial_z \begin{pmatrix} \bar{\rho} v_3' + \bar{v_3} \rho ' \\ \bar{v_3} v_1' \\ \bar{v_3} v_2' \\ \bar{v_3} v_3' + \frac{p'}{\bar{\rho}} \\ \bar{v_3} p' + c^2 \bar{\rho} v_3' \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}\]

The bar $\bar{(\cdot)}$ indicates uniform mean flow variables and $c$ is the speed of sound. The unknowns are the acoustic velocities $v' = (v_1', v_2, v_3')$, the pressure $p'$ and the density $\rho'$.

LobattoLegendreBasis([RealT=Float64,] polydeg::Integer)

Create a nodal Lobatto-Legendre basis for polynomials of degree polydeg.

For the special case polydeg=0 the DG method reduces to a finite volume method. Therefore, this function sets the center point of the cell as single node. This exceptional case is currently only supported for TreeMesh!

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

struct NoSlip

Use to create a no-slip boundary condition with BoundaryConditionNavierStokesWall. The field boundary_value_function should be a function with signature boundary_value_function(x, t, equations) and should return a SVector{NDIMS} whose entries are the velocity vector at a point x and time t.

NonConservativeLocal()

Struct used for multiple dispatch on non-conservative flux functions in the format of "local * symmetric". When the argument nonconservative_type is of type NonConservativeLocal, the function returns the local part of the non-conservative term.

NonConservativeSymmetric()

Struct used for multiple dispatch on non-conservative flux functions in the format of "local * symmetric". When the argument nonconservative_type is of type NonConservativeSymmetric, the function returns the symmetric part of the non-conservative term.

P4estMesh{NDIMS, NDIMS_AMBIENT} <: AbstractMesh{NDIMS}

An unstructured curved mesh based on trees that uses the C library p4est to manage trees and mesh refinement.

The parameter NDIMS denotes the dimension of the spatial domain or manifold represented by the mesh itself, while NDIMS_AMBIENT denotes the dimension of the ambient space in which the mesh is embedded. For example, the type P4estMesh{3, 3} corresponds to a standard mesh for a three-dimensional volume, whereas P4estMesh{2, 3} corresponds to a mesh for a two-dimensional surface or shell in three-dimensional space.

P4estMesh(trees_per_dimension; polydeg,
          mapping=nothing, faces=nothing, coordinates_min=nothing, coordinates_max=nothing,
          RealT=Float64, initial_refinement_level=0, periodicity=true, unsaved_changes=true,
          p4est_partition_allow_for_coarsening=true)

Create a structured curved P4estMesh of the specified size.

There are three ways to map the mesh to the physical domain.

1. Define a mapping that maps the hypercube [-1, 1]^n.
2. Specify a Tuple faces of functions that parametrize each face.
3. Create a rectangular mesh by specifying coordinates_min and coordinates_max.

Non-periodic boundaries will be called :x_neg, :x_pos, :y_neg, :y_pos, :z_neg, :z_pos.

Arguments

• trees_per_dimension::NTupleE{NDIMS, Int}: the number of trees in each dimension.
• polydeg::Integer: polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree.
• mapping: a function of NDIMS variables to describe the mapping that transforms the reference mesh ([-1, 1]^n) to the physical domain. Use only one of mapping, faces and coordinates_min/coordinates_max.
• faces::NTuple{2*NDIMS}: a tuple of 2 * NDIMS functions that describe the faces of the domain. Each function must take NDIMS-1 arguments. faces[1] describes the face onto which the face in negative x-direction of the unit hypercube is mapped. The face in positive x-direction of the unit hypercube will be mapped onto the face described by faces[2]. faces[3:4] describe the faces in positive and negative y-direction respectively (in 2D and 3D). faces[5:6] describe the faces in positive and negative z-direction respectively (in 3D). Use only one of mapping, faces and coordinates_min/coordinates_max.
• coordinates_min: vector or tuple of the coordinates of the corner in the negative direction of each dimension to create a rectangular mesh. Use only one of mapping, faces and coordinates_min/coordinates_max.
• coordinates_max: vector or tuple of the coordinates of the corner in the positive direction of each dimension to create a rectangular mesh. Use only one of mapping, faces and coordinates_min/coordinates_max.
• RealT::Type: the type that should be used for coordinates.
• initial_refinement_level::Integer: refine the mesh uniformly to this level before the simulation starts.
• periodicity: either a Bool deciding if all of the boundaries are periodic or an NTuple{NDIMS, Bool} deciding for each dimension if the boundaries in this dimension are periodic.
• unsaved_changes::Bool: if set to true, the mesh will be saved to a mesh file.
• p4est_partition_allow_for_coarsening::Bool: Must be true when using AMR to make mesh adaptivity independent of domain partitioning. Should be false for static meshes to permit more fine-grained partitioning.

P4estMesh{NDIMS}(meshfile::String;
                 mapping=nothing, polydeg=1, RealT=Float64,
                 initial_refinement_level=0, unsaved_changes=true,
                 p4est_partition_allow_for_coarsening=true,
                 boundary_symbols = nothing)

Main mesh constructor for the P4estMesh that imports an unstructured, conforming mesh from an Abaqus mesh file (.inp). Each element of the conforming mesh parsed from the meshfile is created as a p4est tree datatype.

To create a curved unstructured mesh P4estMesh two strategies are available:

• p4est_mesh_from_hohqmesh_abaqus: High-order, curved boundary information created by HOHQMesh.jl is available in the meshfile. The mesh polynomial degree polydeg of the boundaries is provided from the meshfile. The computation of the mapped tree coordinates is done with transfinite interpolation with linear blending similar to UnstructuredMesh2D. Boundary name information is also parsed from the meshfile such that different boundary conditions can be set at each named boundary on a given tree.
• p4est_mesh_from_standard_abaqus: By default, with mapping=nothing and polydeg=1, this creates a straight-sided from the information parsed from the meshfile. If a mapping function is specified then it computes the mapped tree coordinates via polynomial interpolants with degree polydeg. The mesh created by this function will only have one boundary :all if boundary_symbols is not specified. If boundary_symbols is specified the mesh file will be parsed for nodesets defining the boundary nodes from which boundary edges (2D) and faces (3D) will be assigned.

Note that the mapping and polydeg keyword arguments are only used by the p4est_mesh_from_standard_abaqus function. The p4est_mesh_from_hohqmesh_abaqus function obtains the mesh polydeg directly from the meshfile and constructs the transfinite mapping internally.

The particular strategy is selected according to the header present in the meshfile where the constructor checks whether or not the meshfile was created with HOHQMesh.jl. If the Abaqus file header is not present in the meshfile then the P4estMesh is created with the function p4est_mesh_from_standard_abaqus.

The default keyword argument initial_refinement_level=0 corresponds to a forest where the number of trees is the same as the number of elements in the original meshfile. Increasing the initial_refinement_level allows one to uniformly refine the base mesh given in the meshfile to create a forest with more trees before the simulation begins. For example, if a two-dimensional base mesh contains 25 elements then setting initial_refinement_level=1 creates an initial forest of 2^2 * 25 = 100 trees.

Arguments

• meshfile::String: an uncurved Abaqus mesh file that can be imported by p4est.
• mapping: a function of NDIMS variables to describe the mapping that transforms the imported mesh to the physical domain. Use nothing for the identity map.
• polydeg::Integer: polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree. The default of 1 creates an uncurved geometry. Use a higher value if the mapping will curve the imported uncurved mesh.
• RealT::Type: the type that should be used for coordinates.
• initial_refinement_level::Integer: refine the mesh uniformly to this level before the simulation starts.
• unsaved_changes::Bool: if set to true, the mesh will be saved to a mesh file.
• p4est_partition_allow_for_coarsening::Bool: Must be true when using AMR to make mesh adaptivity independent of domain partitioning. Should be false for static meshes to permit more fine-grained partitioning.
• boundary_symbols::Vector{Symbol}: A vector of symbols that correspond to the boundary names in the meshfile. If nothing is passed then all boundaries are named :all.

PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt = nothing,
                  bS = 1.0, cS = 0.5)
PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization;
                  verbose = false, bS = 1.0, cS = 0.5)
PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64};
                  verbose = false, bS = 1.0, cS = 0.5)
Parameters:
- `num_stages` (`Int`): Number of stages in the PERK method.
- `base_path_monomial_coeffs` (`AbstractString`): Path to a file containing 
  monomial coefficients of the stability polynomial of PERK method.
  The coefficients should be stored in a text file at `joinpath(base_path_monomial_coeffs, "gamma_$(num_stages).txt")` and separated by line breaks.
- `dt_opt` (`Float64`, optional): Optimal time step size for the simulation setup. Can be `nothing` if it is unknown. 
   In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used.
- `tspan`: Time span of the simulation.
- `semi` (`AbstractSemidiscretization`): Semidiscretization setup.
-  `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the
  equation has been semidiscretized.
- `verbose` (`Bool`, optional): Verbosity flag, default is false.
- `bS` (`Float64`, optional): Value of b in the Butcher tableau at b_s, when 
  s is the number of stages, default is 1.0.
- `cS` (`Float64`, optional): Value of c in the Butcher tableau at c_s, when
  s is the number of stages, default is 0.5.

The following structures and methods provide a minimal implementation of the second-order paired explicit Runge-Kutta (PERK) method optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver).

• Brian Vermeire (2019). Paired explicit Runge-Kutta schemes for stiff systems of equations DOI: 10.1016/j.jcp.2019.05.014

Note: To use this integrator, the user must import the Convex and ECOS packages.

PairedExplicitRK3(num_stages, base_path_a_coeffs::AbstractString, dt_opt = nothing;
                  cS2 = 1.0f0)
PairedExplicitRK3(num_stages, tspan, semi::AbstractSemidiscretization;
                  verbose = false, cS2 = 1.0f0)
PairedExplicitRK3(num_stages, tspan, eig_vals::Vector{ComplexF64};
                  verbose = false, cS2 = 1.0f0)

Parameters:
- `num_stages` (`Int`): Number of stages in the paired explicit Runge-Kutta (P-ERK) method.
- `base_path_a_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix in 
  the Butcher tableau of the Runge Kutta method.
  The matrix should be stored in a text file at `joinpath(base_path_a_coeffs, "a_$(num_stages).txt")` and separated by line breaks.
- `dt_opt` (`Float64`, optional): Optimal time step size for the simulation setup. Can be `nothing` if it is unknown. 
   In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used.
- `tspan`: Time span of the simulation.
- `semi` (`AbstractSemidiscretization`): Semidiscretization setup.
-  `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the
  equation has been semidiscretized.
- `verbose` (`Bool`, optional): Verbosity flag, default is false.
- `cS2` (`Float64`, optional): Value of c in the Butcher tableau at c_{s-2}, when
  s is the number of stages, default is 1.0f0.

The following structures and methods provide an implementation of the third-order paired explicit Runge-Kutta (P-ERK) method optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver). The original paper is

• Nasab, Vermeire (2022)

Third-order Paired Explicit Runge-Kutta schemes for stiff systems of equations DOI: 10.1016/j.jcp.2022.111470 While the changes to SSPRK33 base-scheme are described in

• Doehring, Schlottke-Lakemper, Gassner, Torrilhon (2024)

Multirate Time-Integration based on Dynamic ODE Partitioning through Adaptively Refined Meshes for Compressible Fluid Dynamics DOI: 10.1016/j.jcp.2024.113223

Note: To use this integrator, the user must import the Convex, ECOS, and NLsolve packages.

ParametersEulerGravity(; background_density=0.0,
                         gravitational_constant=1.0,
                         cfl=1.0,
                         resid_tol=1.0e-4,
                         n_iterations_max=10^4,
                         timestep_gravity=timestep_gravity_erk52_3Sstar!)

Set up parameters for the gravitational part of a SemidiscretizationEulerGravity.

ParsaniKetchesonDeconinck3Sstar32()

Parsani, Ketcheson, Deconinck (2013) Optimized explicit RK schemes for the spectral difference method applied to wave propagation problems DOI: 10.1137/120885899

ParsaniKetchesonDeconinck3Sstar94()

Parsani, Ketcheson, Deconinck (2013) Optimized explicit RK schemes for the spectral difference method applied to wave propagation problems DOI: 10.1137/120885899

PerformanceCounter()

A PerformanceCounter can be used to track the runtime performance of some calls. Add a new runtime measurement via put!(counter, runtime) and get the averaged runtime of all measurements added so far via take!(counter), resetting the counter.

PerformanceCounterList{N}()

A PerformanceCounterList{N} can be used to track the runtime performance of calls to multiple functions, adding them up. Add a new runtime measurement via put!(counter.counters[i], runtime) and get the averaged runtime of all measurements added so far via take!(counter), resetting the counter.

PlotData1D

Holds all relevant data for creating 1D plots of multiple solution variables and to visualize the mesh.

PlotData1D(u, semi [or mesh, equations, solver, cache];
           solution_variables=nothing, nvisnodes=nothing)

Create a new PlotData1D object that can be used for visualizing 1D DGSEM solution data array u with Plots.jl. All relevant geometrical information is extracted from the semidiscretization semi. By default, the primitive variables (if existent) or the conservative variables (otherwise) from the solution are used for plotting. This can be changed by passing an appropriate conversion function to solution_variables.

nvisnodes specifies the number of visualization nodes to be used. If it is nothing, twice the number of solution DG nodes are used for visualization, and if set to 0, exactly the number of nodes in the DG elements are used.

When visualizing data from a two-dimensional simulation, a 1D slice is extracted for plotting. slice specifies the axis along which the slice is extracted and may be :x, or :y. The slice position is specified by a point that lies on it, which defaults to (0.0, 0.0). Both of these values are ignored when visualizing 1D data. This applies analogously to three-dimensional simulations, where slice may be :xy, :xz, or :yz.

Another way to visualize 2D/3D data is by creating a plot along a given curve. This is done with the keyword argument curve. It can be set to a list of 2D/3D points which define the curve. When using curve any other input from slice or point will be ignored.

PlotData1D(sol; kwargs...)

Create a PlotData1D object from a solution object created by either OrdinaryDiffEq.solve! (which returns a SciMLBase.ODESolution) or Trixi.jl's own solve! (which returns a TimeIntegratorSolution).

PlotData2D

Holds all relevant data for creating 2D plots of multiple solution variables and to visualize the mesh.

PolytropicEulerEquations2D(gamma, kappa)

The polytropic Euler equations

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho v_1 \\ \rho v_2 \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + \kappa\rho^\gamma \\ \rho v_1 v_2 \end{pmatrix} + \frac{\partial}{\partial y} \begin{pmatrix} \rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + \kappa\rho^\gamma \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}\]

for an ideal gas with ratio of specific heats gamma in two space dimensions. Here, $\rho$ is the density and $v_1$ andv_2 the velocities and

\[p = \kappa\rho^\gamma\]

the pressure, which we replaced using this relation.

PositivityPreservingLimiterZhangShu(; threshold, variables)

The fully-discrete positivity-preserving limiter of

• Zhang, Shu (2011) Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments doi: 10.1098/rspa.2011.0153

The limiter is applied to all scalar variables in their given order using the associated thresholds to determine the minimal acceptable values. The order of the variables is important and might have a strong influence on the robustness.

SaveRestartCallback(; interval=0,
                      save_final_restart=true,
                      output_directory="out")

Save the current numerical solution in a restart file every interval time steps.

SaveSolutionCallback(; interval::Integer=0,
                       dt=nothing,
                       save_initial_solution=true,
                       save_final_solution=true,
                       output_directory="out",
                       solution_variables=cons2prim)

Save the current numerical solution in regular intervals. Either pass interval to save every interval time steps or pass dt to save in intervals of dt in terms of integration time by adding additional (shortened) time steps where necessary (note that this may change the solution). solution_variables can be any callable that converts the conservative variables at a single point to a set of solution variables. The first parameter passed to solution_variables will be the set of conservative variables and the second parameter is the equation struct.

SemidiscretizationCoupled

A struct used to bundle multiple semidiscretizations. semidiscretize will return an ODEProblem that synchronizes time steps between the semidiscretizations. Each call of rhs! will call rhs! for each semidiscretization individually. The semidiscretizations can be coupled by gluing meshes together using BoundaryConditionCoupled.

SemidiscretizationCoupled(semis...)

Create a coupled semidiscretization that consists of the semidiscretizations passed as arguments.

SemidiscretizationEulerAcoustics(semi_acoustics::SemiAcoustics, semi_euler::SemiEuler;
                                 source_region=x->true, weights=x->1.0)

Construct a semidiscretization of the acoustic perturbation equations that is coupled with the compressible Euler equations via source terms for the perturbed velocity. Both semidiscretizations have to use the same mesh and solvers with a shared basis. The coupling region is described by a function source_region that maps the coordinates of a single node to true or false depending on whether the point lies within the coupling region or not. A weighting function weights that maps coordinates to weights is applied to the acoustic source terms. Note that this semidiscretization should be used in conjunction with EulerAcousticsCouplingCallback and only works in two dimensions.

SemidiscretizationEulerGravity

A struct containing everything needed to describe a spatial semidiscretization of a the compressible Euler equations with self-gravity, reformulating the Poisson equation for the gravitational potential as steady-state problem of the hyperblic diffusion equations.

• Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020) "A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics" arXiv: 2008.10593

SemidiscretizationEulerGravity(semi_euler::SemiEuler, semi_gravity::SemiGravity, parameters)

Construct a semidiscretization of the compressible Euler equations with self-gravity. parameters should be given as ParametersEulerGravity.

SemidiscretizationHyperbolic

A struct containing everything needed to describe a spatial semidiscretization of a hyperbolic conservation law.

SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
                             source_terms=nothing,
                             boundary_conditions=boundary_condition_periodic,
                             RealT=real(solver),
                             uEltype=RealT,
                             initial_cache=NamedTuple())

Construct a semidiscretization of a hyperbolic PDE.

SemidiscretizationHyperbolicParabolic

A struct containing everything needed to describe a spatial semidiscretization of a mixed hyperbolic-parabolic conservation law.

SemidiscretizationHyperbolicParabolic(mesh, both_equations, initial_condition, solver;
                                      solver_parabolic=default_parabolic_solver(),
                                      source_terms=nothing,
                                      both_boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic),
                                      RealT=real(solver),
                                      uEltype=RealT,
                                      both_initial_caches=(NamedTuple(), NamedTuple()))

Construct a semidiscretization of a hyperbolic-parabolic PDE.

ShallowWaterEquations1D(; gravity, H0 = 0)

Shallow water equations (SWE) in one space dimension. The equations are given by

\[\begin{aligned} \frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h v) &= 0 \\ \frac{\partial}{\partial t}(h v) + \frac{\partial}{\partial x}\left(h v^2 + \frac{g}{2}h^2\right) + g h \frac{\partial b}{\partial x} &= 0 \end{aligned}\]

The unknown quantities of the SWE are the water height $h$ and the velocity $v$. The gravitational constant is denoted by g and the (possibly) variable bottom topography function $b(x)$. Conservative variable water height $h$ is measured from the bottom topography $b$, therefore one also defines the total water height as $H = h + b$.

The additional quantity $H_0$ is also available to store a reference value for the total water height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness.

The bottom topography function $b(x)$ is set inside the initial condition routine for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography variable b to zero.

In addition to the unknowns, Trixi.jl currently stores the bottom topography values at the approximation points despite being fixed in time. This is done for convenience of computing the bottom topography gradients on the fly during the approximation as well as computing auxiliary quantities like the total water height $H$ or the entropy variables. This affects the implementation and use of these equations in various ways:

• The flux values corresponding to the bottom topography must be zero.
• The bottom topography values must be included when defining initial conditions, boundary conditions or source terms.
• AnalysisCallback analyzes this variable.
• Trixi.jl's visualization tools will visualize the bottom topography by default.

References for the SWE are many but a good introduction is available in Chapter 13 of the book:

• Randall J. LeVeque (2002) Finite Volume Methods for Hyperbolic Problems DOI: 10.1017/CBO9780511791253

ShallowWaterEquations2D(; gravity, H0 = 0)

Shallow water equations (SWE) in two space dimensions. The equations are given by

\[\begin{aligned} \frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h v_1) + \frac{\partial}{\partial y}(h v_2) &= 0 \\ \frac{\partial}{\partial t}(h v_1) + \frac{\partial}{\partial x}\left(h v_1^2 + \frac{g}{2}h^2\right) + \frac{\partial}{\partial y}(h v_1 v_2) + g h \frac{\partial b}{\partial x} &= 0 \\ \frac{\partial}{\partial t}(h v_2) + \frac{\partial}{\partial x}(h v_1 v_2) + \frac{\partial}{\partial y}\left(h v_2^2 + \frac{g}{2}h^2\right) + g h \frac{\partial b}{\partial y} &= 0. \end{aligned}\]

The unknown quantities of the SWE are the water height $h$ and the velocities $\mathbf{v} = (v_1, v_2)^T$. The gravitational constant is denoted by g and the (possibly) variable bottom topography function $b(x,y)$. Conservative variable water height $h$ is measured from the bottom topography $b$, therefore one also defines the total water height as $H = h + b$.

The additional quantity $H_0$ is also available to store a reference value for the total water height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness.

The bottom topography function $b(x,y)$ is set inside the initial condition routine for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography variable b to zero.

In addition to the unknowns, Trixi.jl currently stores the bottom topography values at the approximation points despite being fixed in time. This is done for convenience of computing the bottom topography gradients on the fly during the approximation as well as computing auxiliary quantities like the total water height $H$ or the entropy variables. This affects the implementation and use of these equations in various ways:

• The flux values corresponding to the bottom topography must be zero.
• The bottom topography values must be included when defining initial conditions, boundary conditions or source terms.
• AnalysisCallback analyzes this variable.
• Trixi.jl's visualization tools will visualize the bottom topography by default.

References for the SWE are many but a good introduction is available in Chapter 13 of the book:

• Randall J. LeVeque (2002) Finite Volume Methods for Hyperbolic Problems DOI: 10.1017/CBO9780511791253

ShallowWaterEquationsQuasi1D(; gravity, H0 = 0, threshold_limiter = nothing threshold_wet = nothing)

The quasi-1D shallow water equations (SWE). The equations are given by

\[\begin{aligned} \frac{\partial}{\partial t}(a h) + \frac{\partial}{\partial x}(a h v) &= 0 \\ \frac{\partial}{\partial t}(a h v) + \frac{\partial}{\partial x}(a h v^2) + g a h \frac{\partial}{\partial x}(h + b) &= 0 \end{aligned}\]

The unknown quantities of the Quasi-1D SWE are the water height $h$ and the scaled velocity $v$. The gravitational constant is denoted by g, the (possibly) variable bottom topography function $b(x)$, and (possibly) variable channel width $a(x)$. The water height $h$ is measured from the bottom topography $b$, therefore one also defines the total water height as $H = h + b$.

The additional quantity $H_0$ is also available to store a reference value for the total water height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness.

The bottom topography function $b(x)$ and channel width $a(x)$ are set inside the initial condition routine for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography variable b to zero and $a$ to one.

In addition to the unknowns, Trixi.jl currently stores the bottom topography and channel width values at the approximation points despite being fixed in time. This is done for convenience of computing the bottom topography gradients on the fly during the approximation as well as computing auxiliary quantities like the total water height $H$ or the entropy variables. This affects the implementation and use of these equations in various ways:

• The flux values corresponding to the bottom topography and channel width must be zero.
• The bottom topography and channel width values must be included when defining initial conditions, boundary conditions or source terms.
• AnalysisCallback analyzes this variable.
• Trixi.jl's visualization tools will visualize the bottom topography and channel width by default.

SimpleSSPRK33(; stage_callbacks=())

The third-order SSP Runge-Kutta method of Shu and Osher.

References

• Shu, Osher (1988) "Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes" (Eq. 2.18) DOI: 10.1016/0021-9991(88)90177-5

SteadyStateCallback(; abstol=1.0e-8, reltol=1.0e-6)

Terminates the integration when the residual_steady_state(du, equations) falls below the threshold specified by abstol, reltol.

StepsizeCallback(; cfl=1.0)

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

StructuredMesh{NDIMS} <: AbstractMesh{NDIMS}

A structured curved mesh.

Different numbers of cells per dimension are possible and arbitrary functions can be used as domain faces.

StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max; periodicity=true)

Create a StructuredMesh that represents a uncurved structured mesh with a rectangular domain.

Arguments

• cells_per_dimension::NTuple{NDIMS, Int}: the number of cells in each dimension.
• coordinates_min::NTuple{NDIMS, RealT}: coordinate of the corner in the negative direction of each dimension.
• coordinates_max::NTuple{NDIMS, RealT}: coordinate of the corner in the positive direction of each dimension.
• periodicity: either a Bool deciding if all of the boundaries are periodic or an NTuple{NDIMS, Bool} deciding for each dimension if the boundaries in this dimension are periodic.

StructuredMesh(cells_per_dimension, mapping;
               RealT=Float64,
               unsaved_changes=true, 
               mapping_as_string=mapping2string(mapping, length(cells_per_dimension), RealT=RealT))

Create a StructuredMesh of the given size and shape that uses RealT as coordinate type.

Arguments

• cells_per_dimension::NTupleE{NDIMS, Int}: the number of cells in each dimension.
• mapping: a function of NDIMS variables to describe the mapping, which transforms the reference mesh to the physical domain. If no mapping_as_string is defined, this function must be defined with the name mapping to allow for restarts. This will be changed in the future, see https://github.com/trixi-framework/Trixi.jl/issues/541.
• RealT::Type: the type that should be used for coordinates.
• periodicity: either a Bool deciding if all of the boundaries are periodic or an NTuple{NDIMS, Bool} deciding for each dimension if the boundaries in this dimension are periodic.
• unsaved_changes::Bool: if set to true, the mesh will be saved to a mesh file.
• mapping_as_string::String: the code that defines the mapping. If CodeTracking can't find the function definition, it can be passed directly here. The code string must define the mapping function with the name mapping. This will be changed in the future, see https://github.com/trixi-framework/Trixi.jl/issues/541.

StructuredMesh(cells_per_dimension, faces; RealT=Float64, unsaved_changes=true, faces_as_string=faces2string(faces))

Create a StructuredMesh of the given size and shape that uses RealT as coordinate type.

Arguments

• cells_per_dimension::NTupleE{NDIMS, Int}: the number of cells in each dimension.
• faces::NTuple{2*NDIMS}: a tuple of 2 * NDIMS functions that describe the faces of the domain. Each function must take NDIMS-1 arguments. faces[1] describes the face onto which the face in negative x-direction of the unit hypercube is mapped. The face in positive x-direction of the unit hypercube will be mapped onto the face described by faces[2]. faces[3:4] describe the faces in positive and negative y-direction respectively (in 2D and 3D). faces[5:6] describe the faces in positive and negative z-direction respectively (in 3D).
• RealT::Type: the type that should be used for coordinates.
• periodicity: either a Bool deciding if all of the boundaries are periodic or an NTuple{NDIMS, Bool} deciding for each dimension if the boundaries in this dimension are periodic.

StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS}

A view on a structured curved mesh.

StructuredMeshView(parent; indices_min, indices_max)

Create a StructuredMeshView on a StructuredMesh parent.

Arguments

• parent: the parent StructuredMesh.
• indices_min: starting indices of the parent mesh.
• indices_max: ending indices of the parent mesh.

SubcellLimiterIDP(equations::AbstractEquations, basis;
                  local_twosided_variables_cons = String[],
                  positivity_variables_cons = String[],
                  positivity_variables_nonlinear = [],
                  positivity_correction_factor = 0.1,
                  local_onesided_variables_nonlinear = [],
                  max_iterations_newton = 10,
                  newton_tolerances = (1.0e-12, 1.0e-14),
                  gamma_constant_newton = 2 * ndims(equations))

Subcell invariant domain preserving (IDP) limiting used with VolumeIntegralSubcellLimiting including:

• Local two-sided Zalesak-type limiting for conservative variables (local_twosided_variables_cons)
• Positivity limiting for conservative variables (positivity_variables_cons) and nonlinear variables

(positivity_variables_nonlinear)

• Local one-sided limiting for nonlinear variables, e.g. entropy_guermond_etal and entropy_math

with local_onesided_variables_nonlinear

To use these three limiting options use the following structure:

***Conservative variables*** to be limited are passed as a vector of strings, e.g. local_twosided_variables_cons = ["rho"] and positivity_variables_cons = ["rho"]. For ***nonlinear variables***, the wanted variable functions are passed within a vector: To ensure positivity use a plain vector including the desired variables, e.g. positivity_variables_nonlinear = [pressure]. For local one-sided limiting pass the variable function combined with the requested bound (min or max) as a tuple. For instance, to impose a lower local bound on the modified specific entropy by Guermond et al. use local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, min)].

The bounds are calculated using the low-order FV solution. The positivity limiter uses positivity_correction_factor such that u^new >= positivity_correction_factor * u^FV. Local and global limiting of nonlinear variables uses a Newton-bisection method with a maximum of max_iterations_newton iterations, relative and absolute tolerances of newton_tolerances and a provisional update constant gamma_constant_newton (gamma_constant_newton>=2*d, where d = #dimensions). See equation (20) of Pazner (2020) and equation (30) of Rueda-Ramírez et al. (2022).

References

• Rueda-Ramírez, Pazner, Gassner (2022) Subcell Limiting Strategies for Discontinuous Galerkin Spectral Element Methods DOI: 10.1016/j.compfluid.2022.105627
• Pazner (2020) Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting DOI: 10.1016/j.cma.2021.113876

SubcellLimiterIDPCorrection()

Perform antidiffusive correction stage for the a posteriori IDP limiter SubcellLimiterIDP called with VolumeIntegralSubcellLimiting.

References

• Rueda-Ramírez, Pazner, Gassner (2022) Subcell Limiting Strategies for Discontinuous Galerkin Spectral Element Methods DOI: 10.1016/j.compfluid.2022.105627
• Pazner (2020) Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting DOI: 10.1016/j.cma.2021.113876

SurfaceIntegralStrongForm(surface_flux=flux_central)

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

See also VolumeIntegralStrongForm.

SurfaceIntegralUpwind(splitting)

Couple elements with upwind simultaneous approximation terms (SATs) that use a particular flux splitting, e.g., splitting_steger_warming.

See also VolumeIntegralUpwind.

SurfaceIntegralWeakForm(surface_flux=flux_central)

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

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

T8codeMesh{NDIMS} <: AbstractMesh{NDIMS}

An unstructured curved mesh based on trees that uses the C library 't8code' to manage trees and mesh refinement.

T8codeMesh(trees_per_dimension; polydeg, mapping=identity,
           RealT=Float64, initial_refinement_level=0, periodicity=true)

Create a structured potentially curved 'T8codeMesh' of the specified size.

Non-periodic boundaries will be called ':xneg', ':xpos', ':yneg', ':ypos', ':zneg', ':zpos'.

Arguments

• 'treesperdimension::NTupleE{NDIMS, Int}': the number of trees in each dimension.
• 'polydeg::Integer': polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree.
• mapping: a function of NDIMS variables to describe the mapping that transforms the reference mesh ([-1, 1]^n) to the physical domain. Use only one of mapping, faces and coordinates_min/coordinates_max.
• faces::NTuple{2*NDIMS}: a tuple of 2 * NDIMS functions that describe the faces of the domain. Each function must take NDIMS-1 arguments. faces[1] describes the face onto which the face in negative x-direction of the unit hypercube is mapped. The face in positive x-direction of the unit hypercube will be mapped onto the face described by faces[2]. faces[3:4] describe the faces in positive and negative y-direction respectively (in 2D and 3D). faces[5:6] describe the faces in positive and negative z-direction respectively (in 3D). Use only one of mapping, faces and coordinates_min/coordinates_max.
• coordinates_min: vector or tuple of the coordinates of the corner in the negative direction of each dimension to create a rectangular mesh. Use only one of mapping, faces and coordinates_min/coordinates_max.
• coordinates_max: vector or tuple of the coordinates of the corner in the positive direction of each dimension to create a rectangular mesh. Use only one of mapping, faces and coordinates_min/coordinates_max.
• 'RealT::Type': the type that should be used for coordinates.
• 'initialrefinementlevel::Integer': refine the mesh uniformly to this level before the simulation starts.
• 'periodicity': either a 'Bool' deciding if all of the boundaries are periodic or an 'NTuple{NDIMS, Bool}' deciding for each dimension if the boundaries in this dimension are periodic.

T8codeMesh(conn::Ptr{p4est_connectivity}; kwargs...)

Main mesh constructor for the T8codeMesh that imports an unstructured, conforming mesh from a p4est_connectivity data structure.

Arguments

• conn::Ptr{p4est_connectivity}: Pointer to a P4est connectivity object.
• mapping: a function of NDIMS variables to describe the mapping that transforms the imported mesh to the physical domain. Use nothing for the identity map.
• polydeg::Integer: polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree. The default of 1 creates an uncurved geometry. Use a higher value if the mapping will curve the imported uncurved mesh.
• RealT::Type: the type that should be used for coordinates.
• initial_refinement_level::Integer: refine the mesh uniformly to this level before the simulation starts.

T8codeMesh(conn::Ptr{p8est_connectivity}; kwargs...)

Main mesh constructor for the T8codeMesh that imports an unstructured, conforming mesh from a p4est_connectivity data structure.

Arguments

• conn::Ptr{p4est_connectivity}: Pointer to a P4est connectivity object.
• mapping: a function of NDIMS variables to describe the mapping that transforms the imported mesh to the physical domain. Use nothing for the identity map.
• polydeg::Integer: polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree. The default of 1 creates an uncurved geometry. Use a higher value if the mapping will curve the imported uncurved mesh.
• RealT::Type: the type that should be used for coordinates.
• initial_refinement_level::Integer: refine the mesh uniformly to this level before the simulation starts.

T8codeMesh(cmesh::Ptr{t8_cmesh},
           mapping=nothing, polydeg=1, RealT=Float64,
           initial_refinement_level=0)

Main mesh constructor for the T8codeMesh that imports an unstructured, conforming mesh from a t8_cmesh data structure.

Arguments

• cmesh::Ptr{t8_cmesh}: Pointer to a cmesh object.
• mapping: a function of NDIMS variables to describe the mapping that transforms the imported mesh to the physical domain. Use nothing for the identity map.
• polydeg::Integer: polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree. The default of 1 creates an uncurved geometry. Use a higher value if the mapping will curve the imported uncurved mesh.
• RealT::Type: the type that should be used for coordinates.
• initial_refinement_level::Integer: refine the mesh uniformly to this level before the simulation starts.

T8codeMesh(meshfile::String, NDIMS; kwargs...)

Main mesh constructor for the T8codeMesh that imports an unstructured, conforming mesh from either a Gmsh mesh file (.msh) or Abaqus mesh file (.inp) which is determined by the file extension.

Arguments

• filepath::String: path to a Gmsh or Abaqus mesh file.
• ndims: Mesh file dimension: 2 or 3.

Optional Keyword Arguments

• mapping: A function of NDIMS variables to describe the mapping that transforms the imported mesh to the physical domain. Use nothing for the identity map.
• polydeg::Integer: Polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree. The default of 1 creates an uncurved geometry. Use a higher value if the mapping will curve the imported uncurved mesh.
• RealT::Type: The type that should be used for coordinates.
• initial_refinement_level::Integer: Refine the mesh uniformly to this level before the simulation starts.

T8codeMesh{NDIMS, RealT}(forest, boundary_names; polydeg = 1, mapping = nothing)

Main mesh constructor for the T8codeMesh wrapping around a given t8code forest object. This constructor is typically called by other T8codeMesh constructors.

Arguments

• forest: Pointer to a t8code forest.
• boundary_names: List of boundary names.
• polydeg::Integer: Polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree.
• mapping: A function of NDIMS variables to describe the mapping that transforms the imported mesh to the physical domain. Use nothing for the identity map.

T8codeMesh(meshfile::AbaqusFile{NDIMS};
           mapping=nothing, polydeg=1, RealT=Float64,
           initial_refinement_level=0, unsaved_changes=true,
           boundary_symbols = nothing)

Main mesh constructor for the T8codeMesh that imports an unstructured, conforming mesh from an Abaqus mesh file (.inp).

To create a curved unstructured T8codeMesh two strategies are available:

• HOHQMesh Abaqus: High-order, curved boundary information created by HOHQMesh.jl is available in the meshfile. The mesh polynomial degree polydeg of the boundaries is provided from the meshfile. The computation of the mapped tree coordinates is done with transfinite interpolation with linear blending similar to UnstructuredMesh2D. Boundary name information is also parsed from the meshfile such that different boundary conditions can be set at each named boundary on a given tree.

• Standard Abaqus: By default, with mapping=nothing and polydeg=1, this creates a straight-sided mesh from the information parsed from the meshfile. If a mapping function is specified then it computes the mapped tree coordinates via polynomial interpolants with degree polydeg. The mesh created by this function will only have one boundary :all if boundary_symbols is not specified. If boundary_symbols is specified the mesh file will be parsed for nodesets defining the boundary nodes from which boundary edges (2D) and faces (3D) will be assigned.

Note that the mapping and polydeg keyword arguments are only used by the HOHQMesh Abaqus option. The Standard Abaqus routine obtains the mesh polydeg directly from the meshfile and constructs the transfinite mapping internally.

The particular strategy is selected according to the header present in the meshfile where the constructor checks whether or not the meshfile was created with HOHQMesh.jl. If the Abaqus file header is not present in the meshfile then the T8codeMesh is created by Standard Abaqus.

The default keyword argument initial_refinement_level=0 corresponds to a forest where the number of trees is the same as the number of elements in the original meshfile. Increasing the initial_refinement_level allows one to uniformly refine the base mesh given in the meshfile to create a forest with more trees before the simulation begins. For example, if a two-dimensional base mesh contains 25 elements then setting initial_refinement_level=1 creates an initial forest of 2^2 * 25 = 100 trees.

Arguments

• meshfile::AbaqusFile{NDIMS}: Abaqus mesh file object of dimension NDIMS and given path to the file.

Optional Keyword Arguments

• mapping: A function of NDIMS variables to describe the mapping that transforms the imported mesh to the physical domain. Use nothing for the identity map.
• polydeg::Integer: Polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree. The default of 1 creates an uncurved geometry. Use a higher value if the mapping will curve the imported uncurved mesh.
• RealT::Type: The type that should be used for coordinates.
• initial_refinement_level::Integer: Refine the mesh uniformly to this level before the simulation starts.
• boundary_symbols::Vector{Symbol}: A vector of symbols that correspond to the boundary names in the meshfile. If nothing is passed then all boundaries are named :all.

T8codeMesh(meshfile::GmshFile{NDIMS}; kwargs...)

Main mesh constructor for the T8codeMesh that imports an unstructured, conforming mesh from a Gmsh mesh file (.msh).

Arguments

• meshfile::GmshFile{NDIMS}: Gmsh mesh file object of dimension NDIMS and give path to the file.

Optional Keyword Arguments

• mapping: A function of NDIMS variables to describe the mapping that transforms the imported mesh to the physical domain. Use nothing for the identity map.
• polydeg::Integer: Polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree. The default of 1 creates an uncurved geometry. Use a higher value if the mapping will curve the imported uncurved mesh.
• RealT::Type: The type that should be used for coordinates.
• initial_refinement_level::Integer: Refine the mesh uniformly to this level before the simulation starts.

TimeSeriesCallback(semi, point_coordinates;
                   interval=1, solution_variables=cons2cons,
                   output_directory="out", filename="time_series.h5",
                   RealT=real(solver), uEltype=eltype(cache.elements))

Create a callback that records point-wise data at points given in point_coordinates every interval time steps. The point coordinates are to be specified either as a vector of coordinate tuples or as a two-dimensional array where the first dimension is the point number and the second dimension is the coordinate dimension. By default, the conservative variables are recorded, but this can be controlled by passing a different conversion function to solution_variables.

After the last time step, the results are stored in an HDF5 file filename in directory output_directory.

The real data type RealT and data type for solution variables uEltype default to the respective types used in the solver and the cache.

Currently this callback is only implemented for TreeMesh and UnstructuredMesh2D.

TrafficFlowLWREquations1D

The classic Lighthill-Witham Richards (LWR) model for 1D traffic flow. The car density is denoted by $u \in [0, 1]$ and the maximum possible speed (e.g. due to speed limits) is $v_{\text{max}}$.

\[\partial_t u + v_{\text{max}} \partial_1 [u (1 - u)] = 0\]

For more details see e.g. Section 11.1 of

• Randall LeVeque (2002)

Finite Volume Methods for Hyperbolic Problems DOI: 10.1017/CBO9780511791253

TreeMesh{NDIMS} <: AbstractMesh{NDIMS}

A Cartesian mesh based on trees of hypercubes to support adaptive mesh refinement.

UnstructuredMesh2D <: AbstractMesh{2}

An unstructured (possibly curved) quadrilateral mesh.

UnstructuredMesh2D(filename; RealT=Float64, periodicity=false)

All mesh information, neighbour coupling, and boundary curve information is read in from a mesh file filename.

UnstructuredSortedBoundaryTypes

General container to sort the boundary conditions by type and name for some unstructured meshes/solvers. It stores a set of global indices for each boundary condition type and name to expedite computation during the call to calc_boundary_flux!. The original dictionary form of the boundary conditions set by the user in the elixir file is also stored for printing.

ViscousFormulationBassiRebay1()

The classical BR1 flux from

• F. Bassi, S. Rebay (1997) A High-Order Accurate Discontinuous Finite Element Method for the Numerical Solution of the Compressible Navier-Stokes Equations DOI: 10.1006/jcph.1996.5572

ViscousFormulationLocalDG(penalty_parameter)

The local DG (LDG) flux from "The Local Discontinuous Galerkin Method for Time-Dependent Convection-Diffusion Systems" by Cockburn and Shu (1998).

Note that, since this implementation does not involve the parabolic "upwinding" vector, the LDG solver is equivalent to ViscousFormulationBassiRebay1 with an LDG-type penalization.

• Cockburn and Shu (1998). The Local Discontinuous Galerkin Method for Time-Dependent Convection-Diffusion Systems DOI: 10.1137/S0036142997316712

VisualizationCallback(; interval=0,
                        solution_variables=cons2prim,
                        variable_names=[],
                        show_mesh=false,
                        plot_data_creator=PlotData2D,
                        plot_creator=show_plot,
                        plot_arguments...)

Create a callback that visualizes results during a simulation, also known as in-situ visualization.

The interval specifies the number of time step iterations after which a new plot is generated. The available variables to plot are configured with the solution_variables parameter, which acts the same way as for the SaveSolutionCallback. The variables to be actually plotted can be selected by providing a single string or a list of strings to variable_names, and if show_mesh is true, an additional plot with the mesh will be generated.

To customize the generated figure, plot_data_creator allows to use different plot data types. With plot_creator you can further specify an own function to visualize results, which must support the same interface as the default implementation show_plot. All remaining keyword arguments are collected and passed as additional arguments to the plotting command.

VolumeIntegralFluxDifferencing(volume_flux)

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 in Trixi.jl.

References

• LeFloch, Mercier, Rohde (2002) Fully Discrete, Entropy Conservative Schemes of Arbitrary Order doi: 10.1137/S003614290240069X
• Fisher, Carpenter (2013) High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains doi: 10.1016/j.jcp.2013.06.014
• Hendrik Ranocha (2017) Comparison of Some Entropy Conservative Numerical Fluxes for the Euler Equations arXiv: 1701.02264 doi: 10.1007/s10915-017-0618-1
• Chen, Shu (2017) Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws doi: 10.1016/j.jcp.2017.05.025

VolumeIntegralPureLGLFiniteVolume(volume_flux_fv)

A volume integral that only uses the subcell finite volume schemes of the VolumeIntegralShockCapturingHG.

This gives a formally O(1)-accurate finite volume scheme on an LGL-type subcell mesh (LGL = Legendre-Gauss-Lobatto).

References

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

VolumeIntegralShockCapturingHG(indicator; volume_flux_dg=flux_central,
                                          volume_flux_fv=flux_lax_friedrichs)

Shock-capturing volume integral type for DG methods using a convex blending of the finite volume method with numerical flux volume_flux_fv and the VolumeIntegralFluxDifferencing with volume flux volume_flux_dg. The amount of blending is determined by the indicator, e.g., IndicatorHennemannGassner.

References

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

VolumeIntegralStrongForm()

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

VolumeIntegralSubcellLimiting(limiter;
                              volume_flux_dg, volume_flux_fv)

A subcell limiting volume integral type for DG methods based on subcell blending approaches with a low-order FV method. Used with limiter SubcellLimiterIDP.

VolumeIntegralUpwind(splitting)

Specialized volume integral for finite difference summation-by-parts (FDSBP) solvers. Can be used together with the upwind SBP operators of Mattsson (2017) implemented in SummationByPartsOperators.jl. The splitting controls the discretization.

See also splitting_steger_warming, splitting_lax_friedrichs, splitting_vanleer_haenel.

References

• Mattsson (2017) Diagonal-norm upwind SBP operators doi: 10.1016/j.jcp.2017.01.042

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

VolumeIntegralWeakForm() is only implemented for conserved terms as non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, see VolumeIntegralFluxDifferencing. This treatment is required to achieve, e.g., entropy-stability or well-balancedness.

Base.getindex(pd::AbstractPlotData, variable_name)

Extract a single variable variable_name from pd for plotting with Plots.plot.

resize!(c::AbstractContainer, new_length) -> AbstractContainer

Resize c to contain new_length elements. If new_length is smaller than the current container length, the first new_length elements will be retained. If new_length is larger, the new elements are invalidated.

compute_coefficients!(u_ode, func, t, semi::AbstractSemidiscretization)

Same as compute_coefficients but stores the result in u_ode.

compute_coefficients(func, t, semi::AbstractSemidiscretization)

Compute the discrete coefficients of the continuous function func at time t associated with the semidiscretization semi. For example, the discrete coefficients of func for a discontinuous Galerkin spectral element method (DGSEM) are the values of func at the Lobatto-Legendre nodes. Similarly, a classical finite difference method will use the values of func at the nodes of the grid assoociated with the semidiscretization semi.

For semidiscretizations semi associated with an initial condition, func can be omitted to use the given initial condition at time t.

integrate(f, u, basis::LobattoLegendreBasis)

Map the function f to the coefficients u and integrate with respect to the quadrature rule given by basis.

integrate([func=(u_node,equations)->u_node,] u_ode, semi::AbstractSemidiscretization; normalize=true)

Call func(u_node, equations) for each vector of nodal variables u_node in u_ode and integrate the result using a quadrature associated with the semidiscretization semi.

If normalize is true, the result is divided by the total volume of the computational domain.

add_tstop!(integrator::AbstractPairedExplicitRKIntegrator, t)

Add a time stop during the time integration process. This function is called after the periodic SaveSolutionCallback to specify the next stop to save the solution.

add_tstop!(integrator::SimpleIntegratorSSP, t)

Add a time stop during the time integration process. This function is called after the periodic SaveSolutionCallback to specify the next stop to save the solution.

semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan)

Wrap the semidiscretization semi as a split ODE problem in the time interval tspan that can be passed to solve from the SciML ecosystem. The parabolic right-hand side is the first function of the split ODE problem and will be used by default by the implicit part of IMEX methods from the SciML ecosystem.

semidiscretize(semi::AbstractSemidiscretization, tspan, restart_file::AbstractString)

Wrap the semidiscretization semi as an ODE problem in the time interval tspan that can be passed to solve from the SciML ecosystem. The initial condition etc. is taken from the restart_file.

semidiscretize(semi::AbstractSemidiscretization, tspan)

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

DGMultiBasis(element_type, polydeg; approximation_type = Polynomial(), kwargs...)

Constructs a basis for DGMulti solvers. Returns a "StartUpDG.RefElemData" object. The kwargs arguments are additional keyword arguments for RefElemData, such as quad_rule_vol. These are the same as the RefElemData_kwargs used in DGMulti. For more info, see the StartUpDG.jl docs.

LBMCollisionCallback()

Apply the Lattice-Boltzmann method (LBM) collision operator before each time step. See LatticeBoltzmannEquations2D for further details.

P4estMeshCubedSphere(trees_per_face_dimension, layers, inner_radius, thickness;
                     polydeg, RealT=Float64,
                     initial_refinement_level=0, unsaved_changes=true,
                     p4est_partition_allow_for_coarsening=true)

Build a "Cubed Sphere" mesh as P4estMesh with 6 * trees_per_face_dimension^2 * layers trees.

The mesh will have two boundaries, :inside and :outside.

Arguments

• trees_per_face_dimension::Integer: the number of trees in the first two local dimensions of each face.
• layers::Integer: the number of trees in the third local dimension of each face, i.e., the number of layers of the sphere.
• inner_radius::Integer: the inner radius of the sphere.
• thickness::Integer: the thickness of the sphere. The outer radius will be inner_radius + thickness.
• polydeg::Integer: polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree.
• RealT::Type: the type that should be used for coordinates.
• initial_refinement_level::Integer: refine the mesh uniformly to this level before the simulation starts.
• unsaved_changes::Bool: if set to true, the mesh will be saved to a mesh file.
• p4est_partition_allow_for_coarsening::Bool: Must be true when using AMR to make mesh adaptivity independent of domain partitioning. Should be false for static meshes to permit more fine-grained partitioning.

PlotData2D(u, semi [or mesh, equations, solver, cache];
           solution_variables=nothing,
           grid_lines=true, max_supported_level=11, nvisnodes=nothing,
           slice=:xy, point=(0.0, 0.0, 0.0))

Create a new PlotData2D object that can be used for visualizing 2D/3D DGSEM solution data array u with Plots.jl. All relevant geometrical information is extracted from the semidiscretization semi. By default, the primitive variables (if existent) or the conservative variables (otherwise) from the solution are used for plotting. This can be changed by passing an appropriate conversion function to solution_variables.

If grid_lines is true, also extract grid vertices for visualizing the mesh. The output resolution is indirectly set via max_supported_level: all data is interpolated to 2^max_supported_level uniformly distributed points in each spatial direction, also setting the maximum allowed refinement level in the solution. nvisnodes specifies the number of visualization nodes to be used. If it is nothing, twice the number of solution DG nodes are used for visualization, and if set to 0, exactly the number of nodes in the DG elements are used.

When visualizing data from a three-dimensional simulation, a 2D slice is extracted for plotting. slice specifies the plane that is being sliced and may be :xy, :xz, or :yz. The slice position is specified by a point that lies on it, which defaults to (0.0, 0.0, 0.0). Both of these values are ignored when visualizing 2D data.

Examples

julia> using Trixi, Plots

julia> trixi_include(default_example())
[...]

julia> pd = PlotData2D(sol)
PlotData2D(...)

julia> plot(pd) # To plot all available variables

julia> plot(pd["scalar"]) # To plot only a single variable

julia> plot!(getmesh(pd)) # To add grid lines to the plot

PlotData2D(sol; kwargs...)

Create a PlotData2D object from a solution object created by either OrdinaryDiffEq.solve! (which returns a SciMLBase.ODESolution) or Trixi.jl's own solve! (which returns a TimeIntegratorSolution).

ScalarPlotData2D(u, semi::AbstractSemidiscretization; kwargs...)

Returns an PlotData2DTriangulated object which is used to visualize a single scalar field. u should be an array whose entries correspond to values of the scalar field at nodal points.

SummaryCallback()

Create and return a callback that prints a human-readable summary of the simulation setup at the beginning of a simulation and then resets the timer. When the returned callback is executed directly, the current timer values are shown.

TrivialCallback()

A callback that does nothing. This can be useful to disable some callbacks easily via trixi_include.

Trixi.adapt!(mesh::T8codeMesh, adapt_callback; kwargs...)

Adapt a T8codeMesh according to a user-defined adapt_callback.

Arguments

• mesh::T8codeMesh: Initialized mesh object.

• adapt_callback: A user-defined callback which tells the adaption routines if an element should be refined, coarsened or stay unchanged.

  The expected callback signature is as follows:

  `adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, user_data)`
    # Arguments
    - `forest`: Pointer to the analyzed forest.
    - `ltreeid`: Local index of the current tree where the analyzed elements are part of.
    - `eclass_scheme`: Element class of `elements`.
    - `lelemntid`: Local index of the first element in `elements`.
    - `elements`: Array of elements. If consecutive elements form a family
                  they are passed together, otherwise `elements` consists of just one element.
    - `is_family`: Boolean signifying if `elements` represents a family or not.
    - `user_data`: Void pointer to some arbitrary user data. Default value is `C_NULL`.
    # Returns
      -1 : Coarsen family of elements.
       0 : Stay unchanged.
       1 : Refine element.

• kwargs:

  • recursive = true: Adapt the forest recursively. If true the caller must ensure that the callback returns 0 for every analyzed element at some point to stop the recursion.
  • balance = true: Make sure the adapted forest is 2^(NDIMS-1):1 balanced.
  • partition = true: Partition the forest to redistribute elements evenly among MPI ranks.
  • ghost = true: Create a ghost layer for MPI data exchange.
  • user_data = C_NULL: Pointer to some arbitrary user-defined data.

adapt_to_mesh_level!(u_ode, semi, level)
adapt_to_mesh_level!(sol::Trixi.TrixiODESolution, level)

Like adapt_to_mesh_level, but modifies the solution and parts of the semidiscretization (mesh and caches) in place.

adapt_to_mesh_level(u_ode, semi, level)
adapt_to_mesh_level(sol::Trixi.TrixiODESolution, level)

Use the regular adaptive mesh refinement routines to adaptively refine/coarsen the solution u_ode with semidiscretization semi towards a uniformly refined grid with refinement level level. The solution and semidiscretization are copied such that the original objects remain unaltered.

A convenience method accepts an ODE solution object, from which solution and semidiscretization are extracted as needed.

See also: adapt_to_mesh_level!

Trixi.balance!(mesh::T8codeMesh)

Balance a T8codeMesh to ensure 2^(NDIMS-1):1 face neighbors.

barycentric_weights(nodes)

Calculate the barycentric weights for a given node distribution, i.e.,

\[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

boundary_condition_linear_x(u_inner, orientation, direction, x, t,
                            surface_flux_function,
                            equation::LinearScalarAdvectionEquation1D)

Boundary conditions for initial_condition_linear_x.

boundary_condition_linear_x(u_inner, orientation, direction, x, t,
                            surface_flux_function,
                            equation::LinearScalarAdvectionEquation2D)

Boundary conditions for initial_condition_linear_x.

boundary_condition_linear_x_y(u_inner, orientation, direction, x, t,
                              surface_flux_function,
                              equation::LinearScalarAdvectionEquation2D)

Boundary conditions for initial_condition_linear_x_y.

boundary_condition_linear_y(u_inner, orientation, direction, x, t,
                            surface_flux_function,
                            equation::LinearScalarAdvectionEquation2D)

Boundary conditions for initial_condition_linear_y.

boundary_condition_linear_z(u_inner, orientation, direction, x, t,
                            surface_flux_function,
                            equation::LinearScalarAdvectionEquation1D)

Boundary conditions for boundary_condition_linear_z.

boundary_condition_noslip_wall(u_inner, orientation, direction, x, t,
                               surface_flux_function,
                               equations::LatticeBoltzmannEquations2D)

No-slip wall boundary condition using the bounce-back approach.

boundary_condition_poisson_nonperiodic(u_inner, orientation, direction, x, t,
                                       surface_flux_function,
                                       equations::HyperbolicDiffusionEquations1D)

Boundary conditions used for convergence tests in combination with initial_condition_poisson_nonperiodic and source_terms_poisson_nonperiodic.

boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function,
                             equations::AcousticPerturbationEquations2D)

Use an orthogonal projection of the perturbed velocities to zero out the normal velocity while retaining the possibility of a tangential velocity in the boundary state. Further details are available in the paper:

• Marcus Bauer, Jürgen Dierke and Roland Ewert (2011) Application of a discontinuous Galerkin method to discretize acoustic perturbation equations DOI: 10.2514/1.J050333

boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t,
                             surface_flux_function, equations::CompressibleEulerEquations2D)

Should be used together with StructuredMesh.

boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t,
                             surface_flux_function, equations::CompressibleEulerEquations3D)

Should be used together with StructuredMesh.

boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function,
                             equations::CompressibleEulerEquations2D)

Determine the boundary numerical surface flux for a slip wall condition. Imposes a zero normal velocity at the wall. Density is taken from the internal solution state and pressure is computed as an exact solution of a 1D Riemann problem. Further details about this boundary state are available in the paper:

• J. J. W. van der Vegt and H. van der Ven (2002) Slip flow boundary conditions in discontinuous Galerkin discretizations of the Euler equations of gas dynamics PDF

Details about the 1D pressure Riemann solution can be found in Section 6.3.3 of the book

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

Should be used together with UnstructuredMesh2D.

boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function,
                             equations::CompressibleEulerEquations3D)

Determine the boundary numerical surface flux for a slip wall condition. Imposes a zero normal velocity at the wall. Density is taken from the internal solution state and pressure is computed as an exact solution of a 1D Riemann problem. Further details about this boundary state are available in the paper:

• J. J. W. van der Vegt and H. van der Ven (2002) Slip flow boundary conditions in discontinuous Galerkin discretizations of the Euler equations of gas dynamics PDF

Details about the 1D pressure Riemann solution can be found in Section 6.3.3 of the book

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

boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function,
                             equations::ShallowWaterEquations2D)

Create a boundary state by reflecting the normal velocity component and keep the tangential velocity component unchanged. The boundary water height is taken from the internal value. For details see Section 9.2.5 of the book:

• Eleuterio F. Toro (2001) Shock-Capturing Methods for Free-Surface Shallow Flows 1st edition ISBN 0471987662

boundary_condition_slip_wall(u_inner, orientation, direction, x, t,
                             surface_flux_function, equations::CompressibleEulerEquations1D)

Determine the boundary numerical surface flux for a slip wall condition. Imposes a zero normal velocity at the wall. Density is taken from the internal solution state and pressure is computed as an exact solution of a 1D Riemann problem. Further details about this boundary state are available in the paper:

• J. J. W. van der Vegt and H. van der Ven (2002) Slip flow boundary conditions in discontinuous Galerkin discretizations of the Euler equations of gas dynamics PDF

  Should be used together with TreeMesh.

boundary_condition_slip_wall(u_inner, orientation, direction, x, t,
                             surface_flux_function, equations::CompressibleEulerEquations2D)

Should be used together with TreeMesh.

boundary_condition_slip_wall(u_inner, orientation, direction, x, t,
                             surface_flux_function, equations::CompressibleEulerEquations3D)

Should be used together with TreeMesh.

boundary_condition_slip_wall(u_inner, orientation_or_normal, x, t, surface_flux_function,
                              equations::ShallowWaterEquations1D)

Create a boundary state by reflecting the normal velocity component and keep the tangential velocity component unchanged. The boundary water height is taken from the internal value.

For details see Section 9.2.5 of the book:

• Eleuterio F. Toro (2001) Shock-Capturing Methods for Free-Surface Shallow Flows 1st edition ISBN 0471987662

boundary_condition_slip_wall(u_inner, orientation, direction, x, t,
                             surface_flux_function, equations::ShallowWaterEquations2D)

Should be used together with TreeMesh.

boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function,
                        equations::AcousticPerturbationEquations2D)

Boundary conditions for a solid wall.

boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function,
                            equations::LinearizedEulerEquations1D)

Boundary conditions for a solid wall.

boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function,
                            equations::LinearizedEulerEquations2D)

Boundary conditions for a solid wall.

boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function,
                            equations::LinearizedEulerEquations3D)

Boundary conditions for a solid wall.

calc_error_norms([func=(u_node,equations)->u_node,] u_ode, t, analyzer, semi::AbstractSemidiscretization, cache_analysis)

Calculate discrete L2 and L∞ error norms of func applied to each nodal variable u_node in u_ode. If no exact solution is available, "errors" are calculated using some reference state and can be useful for regression tests.

calc_fast_wavespeed_roe(u_ll, u_rr, direction, equations::IdealGlmMhdEquations1D)

Compute the fast magnetoacoustic wave speed using Roe averages as given by

• Cargo and Gallice (1997) Roe Matrices for Ideal MHD and Systematic Construction of Roe Matrices for Systems of Conservation Laws DOI: 10.1006/jcph.1997.5773

calc_fast_wavespeed_roe(u_ll, u_rr, orientation_or_normal_direction, equations::IdealGlmMhdEquations2D)

Compute the fast magnetoacoustic wave speed using Roe averages as given by

• Cargo and Gallice (1997) Roe Matrices for Ideal MHD and Systematic Construction of Roe Matrices for Systems of Conservation Laws DOI: 10.1006/jcph.1997.5773

calc_fast_wavespeed_roe(u_ll, u_rr, orientation_or_normal_direction, equations::IdealGlmMhdEquations3D)

Compute the fast magnetoacoustic wave speed using Roe averages as given by

• Cargo and Gallice (1997) Roe Matrices for Ideal MHD and Systematic Construction of Roe Matrices for Systems of Conservation Laws DOI: 10.1006/jcph.1997.5773

calc_wavespeed_roe(u_ll, u_rr, direction::Integer,
                   equations::ShallowWaterEquations1D)

Calculate Roe-averaged velocity v_roe and wavespeed c_roe = sqrt{g * h_roe} See for instance equation (62) in

• Paul A. Ullrich, Christiane Jablonowski, and Bram van Leer (2010) High-order finite-volume methods for the shallow-water equations on the sphere DOI: 10.1016/j.jcp.2010.04.044

Or equation (9.17) in this lecture notes.

calc_wavespeed_roe(u_ll, u_rr, direction::Integer,
                   equations::ShallowWaterEquations2D)

Calculate Roe-averaged velocity v_roe and wavespeed c_roe = sqrt{g * h_roe} depending on direction. See for instance equation (62) in

• Paul A. Ullrich, Christiane Jablonowski, and Bram van Leer (2010) High-order finite-volume methods for the shallow-water equations on the sphere DOI: 10.1016/j.jcp.2010.04.044

Or this slides, slides 8 and 9.

calculate_cfl(ode_algorithm::AbstractPairedExplicitRK, ode)

This function computes the CFL number once using the initial condition of the problem and the optimal timestep (dt_opt) from the ODE algorithm.

collision_bgk(u, dt, equations::LatticeBoltzmannEquations2D)

Collision operator for the Bhatnagar, Gross, and Krook (BGK) model.

collision_bgk(u, dt, equations::LatticeBoltzmannEquations3D)

Collision operator for the Bhatnagar, Gross, and Krook (BGK) model.

cons2cons(u, equations)

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

cons2entropy(u, equations)

Convert the conserved variables u to the entropy variables for a given set of equations with chosen standard entropy.

u is a vector type of the correct length nvariables(equations). Notice the function doesn't include any error checks for the purpose of efficiency, so please make sure your input is correct. The inverse conversion is performed by entropy2cons.

cons2prim(u, equations)

Convert the conserved variables u to the primitive variables for a given set of equations. u is a vector type of the correct length nvariables(equations). Notice the function doesn't include any error checks for the purpose of efficiency, so please make sure your input is correct. The inverse conversion is performed by prim2cons.

convergence_test([mod::Module=Main,] elixir::AbstractString, iterations; kwargs...)

Run iterations Trixi.jl simulations using the setup given in elixir and compute the experimental order of convergence (EOC) in the $L^2$ and $L^\infty$ norm. In each iteration, the resolution of the respective mesh will be doubled. Additional keyword arguments kwargs... and the optional module mod are passed directly to trixi_include.

This function assumes that the spatial resolution is set via the keywords initial_refinement_level (an integer) or cells_per_dimension (a tuple of integers, one per spatial dimension).

default_analysis_errors(equations)

Default analysis errors (:l2_error and :linf_error) used by the AnalysisCallback.

default_analysis_integrals(equations)

Default analysis integrals used by the AnalysisCallback.

default_example()

Return the path to an example elixir that can be used to quickly see Trixi.jl in action on a TreeMesh. See also examples_dir and get_examples.

default_example_unstructured()

Return the path to an example elixir that can be used to quickly see Trixi.jl in action on an UnstructuredMesh2D. This simulation is run on the example curved, unstructured mesh given in the Trixi.jl documentation regarding unstructured meshes.

density(p::Real, equations::LatticeBoltzmannEquations2D)
density(u, equations::LatticeBoltzmannEquations2D)

Calculate the macroscopic density from the pressure p or the particle distribution functions u.

density(p::Real, equations::LatticeBoltzmannEquations3D)
density(u, equations::LatticeBoltzmannEquations3D)

Calculate the macroscopic density from the pressure p or the particle distribution functions u.

Trixi.download(src_url, file_path)

Download a file from given src_url to given file_path if file_path is not already a file. This function just returns file_path. This is a small wrapper of Downloads.download(src_url, file_path) that avoids race conditions when multiple MPI ranks are used.

dynamic_viscosity(u, equations)

Wrapper for the dynamic viscosity that calls dynamic_viscosity(u, equations.mu, equations), which dispatches on the type of equations.mu. For constant equations.mu, i.e., equations.mu is of Real-type it is returned directly. In all other cases, equations.mu is assumed to be a function with arguments u and equations and is called with these arguments.

each_dof_global(mesh::DGMultiMesh, dg::DGMulti, other_args...)

Return an iterator over the indices that specify the location in relevant data structures for the degrees of freedom (DOF) in dg. In particular, not the DOFs themselves are returned.

each_face_node(mesh::DGMultiMesh, dg::DGMulti, other_args...)

Return an iterator over the indices that specify the location in relevant data structures for the face nodes in dg. In particular, not the face_nodes themselves are returned.

each_face_node_global(mesh::DGMultiMesh, dg::DGMulti, other_args...)

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

each_quad_node(mesh::DGMultiMesh, dg::DGMulti, other_args...)

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

each_quad_node_global(mesh::DGMultiMesh, dg::DGMulti, other_args...)

Return an iterator over the indices that specify the location in relevant data structures for the global quadrature nodes in mesh. In particular, not the quadrature nodes themselves are returned.

eachboundary(dg::DG, cache)

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

eachcomponent(equations::AbstractCompressibleEulerMulticomponentEquations)

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

eachcomponent(equations::AbstractIdealGlmMhdMulticomponentEquations)

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

eachdim(mesh)

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

eachdirection(tree::AbstractTree)

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

eachelement(dg::DG, cache)

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

eachelement(mesh::DGMultiMesh, dg::DGMulti, other_args...)

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

eachelement(elements::ElementContainer1D)

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

eachelement(elements::ElementContainer2D)

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

eachelement(elements::ElementContainer3D)

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

eachelement(elements::UnstructuredElementContainer2D)

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

eachinterface(dg::DG, cache)

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

eachmortar(dg::DG, cache)

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

eachmpiinterface(dg::DG, cache)

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

eachmpimortar(dg::DG, cache)

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

eachnode(dg::DG)

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

eachnode(basis::LobattoLegendreBasis)

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

eachnode(analyzer::LobattoLegendreAnalyzer)

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

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.

energy_internal(u, equations)

Return the internal energy of the conserved variables u for a given set of equations, e.g., the CompressibleEulerEquations2D.

u is a vector of the conserved variables at a single node, i.e., a vector of the correct length nvariables(equations).

energy_kinetic(u, equations)

Return the kinetic energy of the conserved variables u for a given set of equations, e.g., the CompressibleEulerEquations2D.

u is a vector of the conserved variables at a single node, i.e., a vector of the correct length nvariables(equations).

energy_total(u, equations)

Return the total energy of the conserved variables u for a given set of equations, e.g., the CompressibleEulerEquations2D.

u is a vector of the conserved variables at a single node, i.e., a vector of the correct length nvariables(equations).

entropy(u, equations)

Return the chosen entropy of the conserved variables u for a given set of equations.

u is a vector of the conserved variables at a single node, i.e., a vector of the correct length nvariables(equations).

entropy2cons(w, equations)

Convert the entropy variables w based on a standard entropy to the conserved variables for a given set of equations. u is a vector type of the correct length nvariables(equations). Notice the function doesn't include any error checks for the purpose of efficiency, so please make sure your input is correct. The inverse conversion is performed by cons2entropy.

entropy_guermond_etal(u, equations::CompressibleEulerEquations2D)

Calculate the modified specific entropy of Guermond et al. (2019):

\[s_0 = p * \rho^{-\gamma} / (\gamma-1).\]

Note: This is not the "conventional" specific entropy $s = ln(p / \rho^\gamma)$.

• Guermond at al. (2019) Invariant domain preserving discretization-independent schemes and convex limiting for hyperbolic systems. DOI: 10.1016/j.cma.2018.11.036

equilibrium_distribution(alpha, rho, v1, v2, v3, equations::LatticeBoltzmannEquations3D)

Calculate the local equilibrium distribution for the distribution function with index alpha and given the macroscopic state defined by rho, v1, v2, v3.

equilibrium_distribution(alpha, rho, v1, v2, equations::LatticeBoltzmannEquations2D)

Calculate the local equilibrium distribution for the distribution function with index alpha and given the macroscopic state defined by rho, v1, v2.

examples_dir()

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

Examples

readdir(examples_dir())

flux(u, orientation_or_normal, equations)

Given the conservative variables u, calculate the (physical) flux in Cartesian direction orientation::Integer or in arbitrary direction normal::AbstractVector for the corresponding set of governing equations. orientation is 1, 2, and 3 for the x-, y-, and z-directions, respectively.

flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})

Enables calling flux with a non-integer argument normal_direction for one-dimensional equations. Returns the value of flux(u, 1, equations) scaled by normal_direction[1].

flux_central(u_ll, u_rr, orientation_or_normal_direction, 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).

@inline function fluxchanetal(ull, urr, orientation::Integer, equations::CompressibleEulerEquationsQuasi1D)

Conservative (symmetric) part of the entropy conservative flux for quasi 1D compressible Euler equations split form. This flux is a generalization of flux_ranocha for CompressibleEulerEquations1D. Further details are available in the paper:

• Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023) High order entropy stable schemes for the quasi-one-dimensional shallow water and compressible Euler equations DOI: 10.48550/arXiv.2307.12089

flux_chan_etal(u_ll, u_rr, orientation,
               equations::ShallowWaterEquationsQuasi1D)

Total energy conservative (mathematical entropy for quasi 1D shallow water equations) split form. When the bottom topography is nonzero this scheme will be well-balanced when used as a volume_flux. The surface_flux should still use, e.g., FluxPlusDissipation(flux_chan_etal, DissipationLocalLaxFriedrichs()).

Further details are available in the paper:

• Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023) High order entropy stable schemes for the quasi-one-dimensional shallow water and compressible Euler equations DOI: 10.48550/arXiv.2307.12089

flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)

Entropy conserving two-point flux by

• Chandrashekar (2013) Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes for Compressible Euler and Navier-Stokes Equations DOI: 10.4208/cicp.170712.010313a

flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations2D)

Entropy conserving two-point flux by

• Chandrashekar (2013) Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes for Compressible Euler and Navier-Stokes Equations DOI: 10.4208/cicp.170712.010313a

flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations3D)

Entropy conserving two-point flux by

• Chandrashekar (2013) Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes for Compressible Euler and Navier-Stokes Equations DOI: 10.4208/cicp.170712.010313a

flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations1D)

Entropy conserving two-point flux by

• Ayoub Gouasmi, Karthik Duraisamy (2020) "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations" arXiv:1904.00972v3 [math.NA] 4 Feb 2020

flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D)

Adaption of the entropy conserving two-point flux by

• Ayoub Gouasmi, Karthik Duraisamy (2020) "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations" arXiv:1904.00972v3 [math.NA] 4 Feb 2020

flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations1D)

Entropy conserving two-point flux by

• Derigs et al. (2018) Ideal GLM-MHD: About the entropy consistent nine-wave magnetic field divergence diminishing ideal magnetohydrodynamics equations DOI: 10.1016/j.jcp.2018.03.002

flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations2D)

Entropy conserving two-point flux by

• Derigs et al. (2018) Ideal GLM-MHD: About the entropy consistent nine-wave magnetic field divergence diminishing ideal magnetohydrodynamics equations DOI: 10.1016/j.jcp.2018.03.002

flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations3D)

Entropy conserving two-point flux by

• Derigs et al. (2018) Ideal GLM-MHD: About the entropy consistent nine-wave magnetic field divergence diminishing ideal magnetohydrodynamics equations DOI: 10.1016/j.jcp.2018.03.002

flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations1D)

Entropy conserving two-point flux adapted by

• Derigs et al. (2018) Ideal GLM-MHD: About the entropy consistent nine-wave magnetic field divergence diminishing ideal magnetohydrodynamics equations for multicomponent DOI: 10.1016/j.jcp.2018.03.002

flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdMulticomponentEquations2D)

Entropy conserving two-point flux adapted by

• Derigs et al. (2018) Ideal GLM-MHD: About the entropy consistent nine-wave magnetic field divergence diminishing ideal magnetohydrodynamics equations for multicomponent DOI: 10.1016/j.jcp.2018.03.002

flux_fjordholm_etal(u_ll, u_rr, orientation,
                    equations::ShallowWaterEquations1D)

Total energy conservative (mathematical entropy for shallow water equations). When the bottom topography is nonzero this should only be used as a surface flux otherwise the scheme will not be well-balanced. For well-balancedness in the volume flux use flux_wintermeyer_etal.

Details are available in Eq. (4.1) in the paper:

• Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography DOI: 10.1016/j.jcp.2011.03.042

flux_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction,
                    equations::ShallowWaterEquations2D)

Total energy conservative (mathematical entropy for shallow water equations). When the bottom topography is nonzero this should only be used as a surface flux otherwise the scheme will not be well-balanced. For well-balancedness in the volume flux use flux_wintermeyer_etal.

Details are available in Eq. (4.1) in the paper:

• Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography DOI: 10.1016/j.jcp.2011.03.042

flux_godunov(u_ll, u_rr, orientation_or_normal_direction,
             equations::LinearizedEulerEquations2D)

An upwind flux for the linearized Euler equations based on diagonalization of the physical flux matrix. Given the physical flux $Au$, $A=T \Lambda T^{-1}$ with $\Lambda$ being a diagonal matrix that holds the eigenvalues of $A$, decompose $\Lambda = \Lambda^+ + \Lambda^-$ where $\Lambda^+$ and $\Lambda^-$ are diagonal matrices holding the positive and negative eigenvalues of $A$, respectively. Then for left and right states $u_L, u_R$, the numerical flux calculated by this function is given by $A^+ u_L + A^- u_R$ where $A^{\pm} = T \Lambda^{\pm} T^{-1}$.

The diagonalization of the flux matrix can be found in

• R. F. Warming, Richard M. Beam and B. J. Hyett (1975) Diagonalization and simultaneous symmetrization of the gas-dynamic matrices DOI: 10.1090/S0025-5718-1975-0388967-5

flux_hindenlang_gassner(u_ll, u_rr, orientation_or_normal_direction,
                        equations::IdealGlmMhdEquations1D)

Entropy conserving and kinetic energy preserving two-point flux of Hindenlang and Gassner (2019), extending flux_ranocha to the MHD equations.

References

• Florian Hindenlang, Gregor Gassner (2019) A new entropy conservative two-point flux for ideal MHD equations derived from first principles. Presented at HONOM 2019: European workshop on high order numerical methods for evolutionary PDEs, theory and applications
• Hendrik Ranocha (2018) Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods for Hyperbolic Balance Laws PhD thesis, TU Braunschweig
• Hendrik Ranocha (2020) Entropy Conserving and Kinetic Energy Preserving Numerical Methods for the Euler Equations Using Summation-by-Parts Operators Proceedings of ICOSAHOM 2018

flux_hindenlang_gassner(u_ll, u_rr, orientation_or_normal_direction,
                        equations::IdealGlmMhdEquations2D)

Entropy conserving and kinetic energy preserving two-point flux of Hindenlang and Gassner (2019), extending flux_ranocha to the MHD equations.

References

• Florian Hindenlang, Gregor Gassner (2019) A new entropy conservative two-point flux for ideal MHD equations derived from first principles. Presented at HONOM 2019: European workshop on high order numerical methods for evolutionary PDEs, theory and applications
• Hendrik Ranocha (2018) Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods for Hyperbolic Balance Laws PhD thesis, TU Braunschweig
• Hendrik Ranocha (2020) Entropy Conserving and Kinetic Energy Preserving Numerical Methods for the Euler Equations Using Summation-by-Parts Operators Proceedings of ICOSAHOM 2018

flux_hindenlang_gassner(u_ll, u_rr, orientation_or_normal_direction,
                        equations::IdealGlmMhdEquations3D)

Entropy conserving and kinetic energy preserving two-point flux of Hindenlang and Gassner (2019), extending flux_ranocha to the MHD equations.

References

• Florian Hindenlang, Gregor Gassner (2019) A new entropy conservative two-point flux for ideal MHD equations derived from first principles. Presented at HONOM 2019: European workshop on high order numerical methods for evolutionary PDEs, theory and applications
• Hendrik Ranocha (2018) Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods for Hyperbolic Balance Laws PhD thesis, TU Braunschweig
• Hendrik Ranocha (2020) Entropy Conserving and Kinetic Energy Preserving Numerical Methods for the Euler Equations Using Summation-by-Parts Operators Proceedings of ICOSAHOM 2018

flux_hindenlang_gassner(u_ll, u_rr, orientation_or_normal_direction,
                        equations::IdealGlmMhdMulticomponentEquations1D)

Adaption of the entropy conserving and kinetic energy preserving two-point flux of Hindenlang (2019), extending flux_ranocha to the MHD equations.

References

• Florian Hindenlang, Gregor Gassner (2019) A new entropy conservative two-point flux for ideal MHD equations derived from first principles. Presented at HONOM 2019: European workshop on high order numerical methods for evolutionary PDEs, theory and applications
• Hendrik Ranocha (2018) Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods for Hyperbolic Balance Laws PhD thesis, TU Braunschweig
• Hendrik Ranocha (2020) Entropy Conserving and Kinetic Energy Preserving Numerical Methods for the Euler Equations Using Summation-by-Parts Operators Proceedings of ICOSAHOM 2018

flux_hindenlang_gassner(u_ll, u_rr, orientation_or_normal_direction,
                        equations::IdealGlmMhdMulticomponentEquations2D)

Adaption of the entropy conserving and kinetic energy preserving two-point flux of Hindenlang (2019), extending flux_ranocha to the MHD equations.

References

• Florian Hindenlang, Gregor Gassner (2019) A new entropy conservative two-point flux for ideal MHD equations derived from first principles. Presented at HONOM 2019: European workshop on high order numerical methods for evolutionary PDEs, theory and applications
• Hendrik Ranocha (2018) Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods for Hyperbolic Balance Laws PhD thesis, TU Braunschweig
• Hendrik Ranocha (2020) Entropy Conserving and Kinetic Energy Preserving Numerical Methods for the Euler Equations Using Summation-by-Parts Operators Proceedings of ICOSAHOM 2018

flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)

Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro Lecture slides Signal speeds: DOI: 10.1137/S1064827593260140

flux_hllc(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations2D)

Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro Lecture slides Signal speeds: DOI: 10.1137/S1064827593260140

flux_hllc(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations3D)

Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro Lecture slides Signal speeds: DOI: 10.1137/S1064827593260140

flux_hllc(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations1D)

• Li (2005)

An HLLC Riemann solver for magneto-hydrodynamics DOI: 10.1016/j.jcp.2004.08.020.

flux_kennedy_gruber(u_ll, u_rr, orientation, 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_kennedy_gruber(u_ll, u_rr, orientation_or_normal_direction,
                    equations::CompressibleEulerEquations2D)

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_kennedy_gruber(u_ll, u_rr, orientation_or_normal_direction,
                    equations::CompressibleEulerEquations3D)

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_nonconservative_audusse_etal(u_ll, u_rr, orientation::Integer,
                                  equations::ShallowWaterEquations1D)

Non-symmetric two-point surface flux that discretizes the nonconservative (source) term. The discretization uses the hydrostatic_reconstruction_audusse_etal on the conservative variables.

This hydrostatic reconstruction ensures that the finite volume numerical fluxes remain well-balanced for discontinuous bottom topographies ShallowWaterEquations1D. Should be used together with FluxHydrostaticReconstruction and hydrostatic_reconstruction_audusse_etal in the surface flux to ensure consistency.

Further details on the hydrostatic reconstruction and its motivation can be found in

• Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoit Perthame (2004) A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows DOI: 10.1137/S1064827503431090

flux_nonconservative_audusse_etal(u_ll, u_rr, orientation::Integer,
                                  equations::ShallowWaterEquations2D)
flux_nonconservative_audusse_etal(u_ll, u_rr,
                                  normal_direction::AbstractVector,
                                  equations::ShallowWaterEquations2D)

Non-symmetric two-point surface flux that discretizes the nonconservative (source) term. The discretization uses the hydrostatic_reconstruction_audusse_etal on the conservative variables.

This hydrostatic reconstruction ensures that the finite volume numerical fluxes remain well-balanced for discontinuous bottom topographies ShallowWaterEquations2D. Should be used together with FluxHydrostaticReconstruction and hydrostatic_reconstruction_audusse_etal in the surface flux to ensure consistency.

Further details for the hydrostatic reconstruction and its motivation can be found in

• Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoit Perthame (2004) A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows DOI: 10.1137/S1064827503431090

flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
                               equations::CompressibleEulerEquationsQuasi1D)
flux_nonconservative_chan_etal(u_ll, u_rr, normal_direction, 
                               equations::CompressibleEulerEquationsQuasi1D)
flux_nonconservative_chan_etal(u_ll, u_rr, normal_ll, normal_rr,
                               equations::CompressibleEulerEquationsQuasi1D)

Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the pressure CompressibleEulerEquationsQuasi1D and the nozzle width.

Further details are available in the paper:

• Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023) High order entropy stable schemes for the quasi-one-dimensional shallow water and compressible Euler equations DOI: 10.48550/arXiv.2307.12089

flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
                               equations::ShallowWaterEquationsQuasi1D)
flux_nonconservative_chan_etal(u_ll, u_rr, normal_direction::AbstractVector,
                               equations::ShallowWaterEquationsQuasi1D)    
flux_nonconservative_chan_etal(u_ll, u_rr, 
                               normal_ll::AbstractVector, normal_rr::AbstractVector,
                               equations::ShallowWaterEquationsQuasi1D)

Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography ShallowWaterEquationsQuasi1D and the channel width.

Further details are available in the paper:

• Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023) High order entropy stable schemes for the quasi-one-dimensional shallow water and compressible Euler equations DOI: 10.48550/arXiv.2307.12089

flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer,
                                    equations::ShallowWaterEquations1D)

Non-symmetric two-point surface flux discretizing the nonconservative (source) term of that contains the gradient of the bottom topography ShallowWaterEquations1D.

This flux can be used together with flux_fjordholm_etal at interfaces to ensure entropy conservation and well-balancedness.

Further details for the original finite volume formulation are available in

• Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography DOI: 10.1016/j.jcp.2011.03.042

and for curvilinear 2D case in the paper:

• Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry DOI: 10.1016/j.jcp.2017.03.036

flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer,
                                    equations::ShallowWaterEquations2D)
flux_nonconservative_fjordholm_etal(u_ll, u_rr,
                                    normal_direction::AbstractVector,
                                    equations::ShallowWaterEquations2D)

Non-symmetric two-point surface flux discretizing the nonconservative (source) term of that contains the gradient of the bottom topography ShallowWaterEquations2D.

This flux can be used together with flux_fjordholm_etal at interfaces to ensure entropy conservation and well-balancedness.

Further details for the original finite volume formulation are available in

• Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography DOI: 10.1016/j.jcp.2011.03.042

and for curvilinear 2D case in the paper:

• Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry DOI: 10.1016/j.jcp.2017.03.036

flux_nonconservative_powell(u_ll, u_rr, orientation::Integer,
                            equations::IdealGlmMhdEquations2D)
flux_nonconservative_powell(u_ll, u_rr,
                            normal_direction::AbstractVector,
                            equations::IdealGlmMhdEquations2D)

Non-symmetric two-point flux discretizing the nonconservative (source) term of Powell and the Galilean nonconservative term associated with the GLM multiplier of the IdealGlmMhdEquations2D.

On curvilinear meshes, the implementation differs from the reference since we use the averaged normal direction for the metrics dealiasing.

References

• Marvin Bohm, Andrew R.Winters, Gregor J. Gassner, Dominik Derigs, Florian Hindenlang, Joachim Saur An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. Part I: Theory and numerical verification DOI: 10.1016/j.jcp.2018.06.027

flux_nonconservative_powell(u_ll, u_rr, orientation::Integer,
                            equations::IdealGlmMhdEquations3D)
flux_nonconservative_powell(u_ll, u_rr,
                            normal_direction::AbstractVector,
                            equations::IdealGlmMhdEquations3D)

Non-symmetric two-point flux discretizing the nonconservative (source) term of Powell and the Galilean nonconservative term associated with the GLM multiplier of the IdealGlmMhdEquations3D.

On curvilinear meshes, the implementation differs from the reference since we use the averaged normal direction for the metrics dealiasing.

References

• Marvin Bohm, Andrew R.Winters, Gregor J. Gassner, Dominik Derigs, Florian Hindenlang, Joachim Saur An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. Part I: Theory and numerical verification DOI: 10.1016/j.jcp.2018.06.027

flux_nonconservative_powell(u_ll, u_rr, orientation::Integer,
                            equations::IdealGlmMhdMulticomponentEquations2D)

Non-symmetric two-point flux discretizing the nonconservative (source) term of Powell and the Galilean nonconservative term associated with the GLM multiplier of the IdealGlmMhdMulticomponentEquations2D.

References

• Marvin Bohm, Andrew R.Winters, Gregor J. Gassner, Dominik Derigs, Florian Hindenlang, Joachim Saur An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. Part I: Theory and numerical verification DOI: 10.1016/j.jcp.2018.06.027

flux_nonconservative_powell_local_symmetric(u_ll, orientation::Integer,
                                            equations::IdealGlmMhdEquations2D,
                                            nonconservative_type::NonConservativeSymmetric,
                                            nonconservative_term::Integer)

Symmetric part of the Powell and GLM non-conservative terms. Needed for the calculation of the non-conservative staggered "fluxes" for subcell limiting. See, e.g.,

• Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf.

This function is used to compute the subcell fluxes in dg2dsubcell_limiters.jl.

flux_nonconservative_powell_local_symmetric(u_ll, u_rr,
                                            orientation::Integer,
                                            equations::IdealGlmMhdEquations2D)

Non-symmetric two-point flux discretizing the nonconservative (source) term of Powell and the Galilean nonconservative term associated with the GLM multiplier of the IdealGlmMhdEquations2D.

This implementation uses a non-conservative term that can be written as the product of local and symmetric parts. It is equivalent to the non-conservative flux of Bohm et al. flux_nonconservative_powell for conforming meshes but it yields different results on non-conforming meshes(!). On curvilinear meshes this formulation applies the local normal direction compared to the averaged one used in flux_nonconservative_powell.

The two other flux functions with the same name return either the local or symmetric portion of the non-conservative flux based on the type of the nonconservativetype argument, employing multiple dispatch. They are used to compute the subcell fluxes in dg2dsubcelllimiters.jl.

References

• Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf.

flux_nonconservative_powell_local_symmetric(u_ll, orientation::Integer,
                                            equations::IdealGlmMhdEquations2D,
                                            nonconservative_type::NonConservativeLocal,
                                            nonconservative_term::Integer)

Local part of the Powell and GLM non-conservative terms. Needed for the calculation of the non-conservative staggered "fluxes" for subcell limiting. See, e.g.,

• Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf.

This function is used to compute the subcell fluxes in dg2dsubcell_limiters.jl.

flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
                                      equations::ShallowWaterEquations1D)

Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography ShallowWaterEquations1D.

Gives entropy conservation and well-balancedness on both the volume and surface when combined with flux_wintermeyer_etal.

Further details are available in the papers:

• Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry DOI: 10.1016/j.jcp.2017.03.036
• Patrick Ersing, Andrew R. Winters (2023) An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on curvilinear meshes DOI: 10.48550/arXiv.2306.12699

flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
                                      equations::ShallowWaterEquations2D)
flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
                                      normal_direction::AbstractVector,
                                      equations::ShallowWaterEquations2D)

Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography ShallowWaterEquations2D.

For the surface_flux either flux_wintermeyer_etal or flux_fjordholm_etal can be used to ensure well-balancedness and entropy conservation.

Further details are available in the papers:

• Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry DOI: 10.1016/j.jcp.2017.03.036
• Patrick Ersing, Andrew R. Winters (2023) An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on curvilinear meshes DOI: 10.48550/arXiv.2306.12699

flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, 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

flux_ranocha(u_ll, u_rr, orientation_or_normal_direction,
             equations::CompressibleEulerEquations2D)

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

flux_ranocha(u_ll, u_rr, orientation_or_normal_direction,
             equations::CompressibleEulerEquations3D)

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

flux_ranocha(u_ll, u_rr, orientation_or_normal_direction,
             equations::CompressibleEulerMulticomponentEquations1D)

Adaption of the 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

flux_ranocha(u_ll, u_rr, orientation_or_normal_direction,
             equations::CompressibleEulerMulticomponentEquations2D)

Adaption of the 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

flux_ranocha_turbo(u_ll, u_rr, orientation_or_normal_direction, equations)

Equivalent to flux_ranocha except that it may use specialized methods, e.g., when used with VolumeIntegralFluxDifferencing. These specialized methods may enable better use of SIMD instructions to increase runtime efficiency on modern hardware.

flux_shima_etal(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)

This flux is is a modification of the original kinetic energy preserving two-point flux by

• Yuichi Kuya, Kosuke Totani and Soshi Kawai (2018) Kinetic energy and entropy preserving schemes for compressible flows by split convective forms DOI: 10.1016/j.jcp.2018.08.058

The modification is in the energy flux to guarantee pressure equilibrium and was developed by

• Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020) Preventing spurious pressure oscillations in split convective form discretizations for compressible flows DOI: 10.1016/j.jcp.2020.110060

flux_shima_etal(u_ll, u_rr, orientation_or_normal_direction,
                equations::CompressibleEulerEquations2D)

This flux is is a modification of the original kinetic energy preserving two-point flux by

• Yuichi Kuya, Kosuke Totani and Soshi Kawai (2018) Kinetic energy and entropy preserving schemes for compressible flows by split convective forms DOI: 10.1016/j.jcp.2018.08.058

The modification is in the energy flux to guarantee pressure equilibrium and was developed by

• Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020) Preventing spurious pressure oscillations in split convective form discretizations for compressible flows DOI: 10.1016/j.jcp.2020.110060

flux_shima_etal(u_ll, u_rr, orientation_or_normal_direction,
                equations::CompressibleEulerEquations3D)

This flux is is a modification of the original kinetic energy preserving two-point flux by

• Yuichi Kuya, Kosuke Totani and Soshi Kawai (2018) Kinetic energy and entropy preserving schemes for compressible flows by split convective forms DOI: 10.1016/j.jcp.2018.08.058

The modification is in the energy flux to guarantee pressure equilibrium and was developed by

• Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020) Preventing spurious pressure oscillations in split convective form discretizations for compressible flows DOI: 10.1016/j.jcp.2020.110060

flux_shima_etal_turbo(u_ll, u_rr, orientation_or_normal_direction, equations)

Equivalent to flux_shima_etal except that it may use specialized methods, e.g., when used with VolumeIntegralFluxDifferencing. These specialized methods may enable better use of SIMD instructions to increase runtime efficiency on modern hardware.

flux_wintermeyer_etal(u_ll, u_rr, orientation,
                      equations::ShallowWaterEquations1D)

Total energy conservative (mathematical entropy for shallow water equations) split form. When the bottom topography is nonzero this scheme will be well-balanced when used as a volume_flux. For the surface_flux either flux_wintermeyer_etal or flux_fjordholm_etal can be used to ensure well-balancedness and entropy conservation.

Further details are available in Theorem 1 of the paper:

• Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry DOI: 10.1016/j.jcp.2017.03.036

flux_wintermeyer_etal(u_ll, u_rr, orientation_or_normal_direction,
                      equations::ShallowWaterEquations2D)

Total energy conservative (mathematical entropy for shallow water equations) split form. When the bottom topography is nonzero this scheme will be well-balanced when used as a volume_flux. For the surface_flux either flux_wintermeyer_etal or flux_fjordholm_etal can be used to ensure well-balancedness and entropy conservation.

Further details are available in Theorem 1 of the paper:

• Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry DOI: 10.1016/j.jcp.2017.03.036

flux_winters_etal(u_ll, u_rr, orientation_or_normal_direction,
                  equations::PolytropicEulerEquations2D)

Entropy conserving two-point flux for isothermal or polytropic gases. Requires a special weighted Stolarsky mean for the evaluation of the density denoted here as stolarsky_mean. Note, for isothermal gases where gamma = 1 this stolarsky_mean becomes the ln_mean.

For details see Section 3.2 of the following reference

• Andrew R. Winters, Christof Czernik, Moritz B. Schily & Gregor J. Gassner (2020) Entropy stable numerical approximations for the isothermal and polytropic Euler equations DOI: 10.1007/s10543-019-00789-w

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

get_boundary_outer_state(u_inner, t,
                         boundary_condition::BoundaryConditionDirichlet,
                         orientation_or_normal, direction,
                         equations, dg, cache, indices...)

For subcell limiting, the calculation of local bounds for non-periodic domains requires the boundary outer state. This function returns the boundary value for BoundaryConditionDirichlet at time t and for node with spatial indices indices at the boundary with orientation_or_normal and direction.

Should be used together with TreeMesh or StructuredMesh.

get_examples()

Return a list of all example elixirs that are provided by Trixi.jl. See also examples_dir and default_example.

get_name(x)

Returns a name of x ready for pretty printing. By default, return string(y) if x isa Val{y} and return string(x) otherwise.

Examples

julia> Trixi.get_name("test")
"test"

julia> Trixi.get_name(Val(:test))
"test"

get_name(equations::AbstractEquations)

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

Examples

julia> Trixi.get_name(CompressibleEulerEquations1D(1.4))
"CompressibleEulerEquations1D"

getmesh(pd::AbstractPlotData)

Extract grid lines from pd for plotting with Plots.plot.

global_mean_vars(equations::AcousticPerturbationEquations2D)

Returns the global mean variables stored in equations. This makes it easier to define flexible initial conditions for problems with constant mean flow.

have_nonconservative_terms(equations)

Trait function determining whether equations represent a conservation law with or without nonconservative terms. Classical conservation laws such as the CompressibleEulerEquations2D do not have nonconservative terms. The ShallowWaterEquations2D with non-constant bottom topography are an example of equations with nonconservative terms. The return value will be True() or False() to allow dispatching on the return type.

hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, orientation::Integer,
                                        equations::ShallowWaterEquations1D)

A particular type of hydrostatic reconstruction on the water height to guarantee well-balancedness for a general bottom topography ShallowWaterEquations1D. The reconstructed solution states u_ll_star and u_rr_star variables are then used to evaluate the surface numerical flux at the interface. Use in combination with the generic numerical flux routine FluxHydrostaticReconstruction.

Further details on this hydrostatic reconstruction and its motivation can be found in

• Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoit Perthame (2004) A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows DOI: 10.1137/S1064827503431090

hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, orientation_or_normal_direction,
                                        equations::ShallowWaterEquations2D)

A particular type of hydrostatic reconstruction on the water height to guarantee well-balancedness for a general bottom topography ShallowWaterEquations2D. The reconstructed solution states u_ll_star and u_rr_star variables are used to evaluate the surface numerical flux at the interface. Use in combination with the generic numerical flux routine FluxHydrostaticReconstruction.

Further details for the hydrostatic reconstruction and its motivation can be found in

• Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoit Perthame (2004) A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows DOI: 10.1137/S1064827503431090

init_mpi()

Initialize MPI by calling MPI.Initialized(). The function will check if MPI is already initialized and if yes, do nothing, thus it is safe to call it multiple times.

init_p4est()

Initialize p4est by calling p4est_init and setting the log level to SC_LP_ERROR. This function will check if p4est is already initialized and if yes, do nothing, thus it is safe to call it multiple times.

init_t8code()

Initialize t8code by calling sc_init, p4est_init, and t8_init while setting the log level to SC_LP_ERROR. This function will check if t8code is already initialized and if yes, do nothing, thus it is safe to call it multiple times.

initial_condition_constant(x, t, equations::AcousticPerturbationEquations2D)

A constant initial condition where the state variables are zero and the mean flow is constant. Uses the global mean values from equations.

initial_condition_constant(x, t, equations::CompressibleEulerEquations1D)

A constant initial condition to test free-stream preservation.

initial_condition_constant(x, t, equations::CompressibleEulerEquations2D)

A constant initial condition to test free-stream preservation.

initial_condition_constant(x, t, equations::CompressibleEulerEquations3D)

A constant initial condition to test free-stream preservation.

initial_condition_constant(x, t, equations::IdealGlmMhdEquations1D)

A constant initial condition to test free-stream preservation.

initial_condition_constant(x, t, equations::IdealGlmMhdEquations2D)

A constant initial condition to test free-stream preservation.

initial_condition_constant(x, t, equations::IdealGlmMhdEquations3D)

A constant initial condition to test free-stream preservation.

initial_condition_constant(x, t, equations::InviscidBurgersEquation1D)

A constant initial condition to test free-stream preservation.

initial_condition_constant(x, t, equations::LatticeBoltzmannEquations2D)

A constant initial condition to test free-stream preservation.

initial_condition_constant(x, t, equations::LatticeBoltzmannEquations3D)

A constant initial condition to test free-stream preservation.

initial_condition_constant(x, t, equations::LinearScalarAdvectionEquation1D)

A constant initial condition to test free-stream preservation.

initial_condition_constant(x, t, equations::LinearScalarAdvectionEquation2D)

A constant initial condition to test free-stream preservation.

initial_condition_constant(x, t, equations::LinearScalarAdvectionEquation1D)

A constant initial condition to test free-stream preservation.

initial_condition_convergence_test(x, t, equations::AcousticPerturbationEquations2D)

A smooth initial condition used for convergence tests in combination with source_terms_convergence_test. Uses the global mean values from equations.

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_convergence_test(x, t, equations::CompressibleEulerEquations2D)

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_convergence_test(x, t, equations::CompressibleEulerEquations3D)

A smooth initial condition used for convergence tests in combination with source_terms_convergence_test.

initial_condition_convergence_test(x, t, equations::CompressibleEulerEquationsQuasi1D)

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_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations1D)

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_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations2D)

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_convergence_test(x, t, equations::IdealGlmMhdEquations1D)

An Alfvén wave as smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquations2D)

An Alfvén wave as smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquations3D)

An Alfvén wave as smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::IdealGlmMhdMulticomponentEquations1D)

An Alfvén wave as smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::IdealGlmMhdMulticomponentEquations2D)

An Alfvén wave as smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::InviscidBurgersEquation1D)

A smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::LinearScalarAdvectionEquation1D)

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

initial_condition_convergence_test(x, t, equations::LinearScalarAdvectionEquation2D)

A smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::LinearScalarAdvectionEquation1D)

A smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::LinearizedEulerEquations1D)

A smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::LinearizedEulerEquations2D)

A smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::LinearizedEulerEquations3D)

A smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::MaxwellEquations1D)

A smooth initial condition used for convergence tests.

initial_condition_convergence_test(x, t, equations::PolytropicEulerEquations2D)

Manufactured smooth initial condition used for convergence tests in combination with source_terms_convergence_test.

initial_condition_convergence_test(x, t, equations::ShallowWaterEquations1D)

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_convergence_test(x, t, equations::ShallowWaterEquations2D)

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_convergence_test(x, t, equations::ShallowWaterEquationsQuasi1D)

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_convergence_test(x, t, equations::TrafficFlowLWREquations1D)

A smooth initial condition used for convergence tests.

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 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_density_wave(x, t, equations::CompressibleEulerEquations2D)

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 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_eoc_test_coupled_euler_gravity(x, t, equations::CompressibleEulerEquations1D)

One dimensional variant of the setup used for convergence tests of the Euler equations with self-gravity from

• Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020) A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics arXiv: 2008.10593

initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::CompressibleEulerEquations2D)

Setup used for convergence tests of the Euler equations with self-gravity used in

• Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020) A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics arXiv: 2008.10593

in combination with source_terms_eoc_test_coupled_euler_gravity or source_terms_eoc_test_euler.

initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::CompressibleEulerEquations3D)

Setup used for convergence tests of the Euler equations with self-gravity used in

• Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020) A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics arXiv: 2008.10593

in combination with source_terms_eoc_test_coupled_euler_gravity or source_terms_eoc_test_euler.

initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::HyperbolicDiffusionEquations1D)

Setup used for convergence tests of the Euler equations with self-gravity used in

• Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020) A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics arXiv: 2008.10593

in combination with source_terms_harmonic.

initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::HyperbolicDiffusionEquations2D)

Setup used for convergence tests of the Euler equations with self-gravity used in

• Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020) A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics arXiv: 2008.10593

in combination with source_terms_harmonic.

initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::HyperbolicDiffusionEquations3D)

Setup used for convergence tests of the Euler equations with self-gravity used in

• Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020) A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics arXiv: 2008.10593

in combination with source_terms_harmonic.

initial_condition_gauss(x, t, equations::AcousticPerturbationEquations2D)

A Gaussian pulse in a constant mean flow. Uses the global mean values from equations.

initial_condition_gauss(x, t, equations::LinearScalarAdvectionEquation1D)

A Gaussian pulse used together with BoundaryConditionDirichlet(initial_condition_gauss).

initial_condition_gauss(x, t, equation::LinearScalarAdvectionEquation2D)

A Gaussian pulse used together with BoundaryConditionDirichlet(initial_condition_gauss).

initial_condition_gauss(x, t, equations::LinearScalarAdvectionEquation1D)

A Gaussian pulse.

initial_condition_linear_x(x, t, equations::LinearScalarAdvectionEquation1D)

A linear function of x[1] used together with boundary_condition_linear_x.

initial_condition_linear_x(x, t, equations::LinearScalarAdvectionEquation2D)

A linear function of x[1] used together with boundary_condition_linear_x.

initial_condition_linear_x_y(x, t, equations::LinearScalarAdvectionEquation2D)

A linear function of x[1] + x[2] used together with boundary_condition_linear_x_y.

initial_condition_linear_y(x, t, equations::LinearScalarAdvectionEquation2D)

A linear function of x[1] used together with boundary_condition_linear_y.

initial_condition_linear_z(x, t, equations::LinearScalarAdvectionEquation1D)

A linear function of x[3] used together with boundary_condition_linear_z.

initial_condition_poisson_nonperiodic(x, t, equations::HyperbolicDiffusionEquations1D)

A non-priodic smooth initial condition. Can be used for convergence tests in combination with source_terms_poisson_nonperiodic and boundary_condition_poisson_nonperiodic.