API Reference

Struct to hold precomputed constants for Gramian moment equations in 1D with 3D velocity.

# Fields
- `fp_func::F`: Compiled function for evaluating basis functions at given angles.
- `Angles::Vector{Tuple{Float64,Float64}}`: List of angles (θ, φ) for quadrature.
- `A_matrix::Matrix{Float64}`: Matrix A used in closure transformation (solves the linear system).
- `slab_indices::Vector{Int}`: Indices of valid moments in the slab geometry.
- `target_indices::Vector{Int}`: Indices of moments used for closure.
- `rotations::AbstractVector{<:AbstractMatrix{Float64}}`: Precomputed rotation matrices for each angle.

Discrete Velocity Method (DVM) Equations in 1D velocity space with N discrete velocity points.

∂_t f_k + c_k ∂_x f_k = - (1/Kn) (f_k - f̃_k),  k = 1,...,N

where f̃_k is the corrected Maxwellian distribution ensuring conservation of mass, momentum, and energy.
Discretizes the velocity space using N discrete velocity points between c_l and c_u.

# Fields
- `N::Int`: Number of equations to be solved.
- `c_l::Real`: Lower bound of the velocity domain.
- `c_u::Real`: Upper bound of the velocity domain.
- `c_vec::SVector{N, Float64}`: Velocity grid points.
- `Kn::Real`: Knudsen number.

Storage for the electric field values and related data for the Vlasov-Poisson system

# Fields:
- `E::Vector{Float64}`: Electric field values at spatial points
- `Lx::Float64`: Length of the domain (x_max - x_min) for Poisson solve
- `x_range::Vector{Float64}`: Spatial grid points
- `initialized::Bool`: Flag to check if initialized
- `ρ::Vector{Float64}`: Density values at spatial points (to be removed later)
- `MP1::Int`: Number of moments + 1
- `E_L2::Vector{Float64}`: Track L2 norm of electric field over time
- `times::Vector{Float64}`: Track actual times when E_L2 is recorded
- `variables`: Store solution variables at spatial points
- `coordinates::Vector{Float64}`: Store spatial coordinates

Abstract type for the Gramian Moment Equations. This type is a subtype of Trixi.AbstractEquations and represents the system of equations for the Gramian moments.

This class serves as a marker for the Gramian Moment Equations and does not contain any fields. The specific implementation of the equations is done in the subtypes, such as GramianMomentEquations and GramianMomentEquations3D. It allows to use different closures.

struct GramianMomentEquations1D{Mp1, N, RealT <: Real} <: Trixi.AbstractEquations{1, Mp1}

A struct representing one-dimensional Gramian moment equations for kinetic theory.

closure::Symbol

• inv_Kn::RealT: The inverse of the Knudsen number, used in the relaxation source term function GramianMomentEquations1D(M::Integer, Knudsen::Real, closure::Union{Symbol,String}=:ExtGram; χ_set="optimal")
• n::Int: Internal parameter derived from the number of moments
• extended::Bool: Flag indicating whether to use extended formulation

Constructor

GramianMomentEquations1D(M::Integer, Knudsen::Real, extended=true::Bool)

Constructs a GramianMomentEquations1D system with M moments and specified Knudsen number.

Arguments

• M::Integer: Number of moments (must be greater than 1)
• Knudsen::Real: Knudsen number for the system
• closure::Union{Symbol,String}: Closure type, can be "Gram", "ExtGram", or "Grad"

Notes

• The system has M+1 variables in total
• The parameter χ is calculated differently depending on whether M is even or odd

struct GramianMomentEquations1D3D{Mp1, N, RealT <: Real} <: Trixi.AbstractEquations{1, Mp1}

A struct representing one-dimensional Gramian moment equations for kinetic theory.

closure::Symbol

• inv_Kn::RealT: The inverse of the Knudsen number, used in the relaxation source term function GramianMomentEquations1D(M::Integer, Knudsen::Real, closure::Union{Symbol,String}=:ExtGram; χ_set="optimal")
• n::Int: Internal parameter derived from the number of moments
• extended::Bool: Flag indicating whether to use extended formulation

Constructor

GramianMomentEquations1D(M::Integer, Knudsen::Real, extended=true::Bool)

Constructs a GramianMomentEquations1D system with M moments and specified Knudsen number.

Arguments

• M::Integer: Number of moments (must be greater than 1) - for the moment hard-coded to M=4
• Knudsen::Real: Knudsen number for the system
• closure::Union{Symbol,String}: Closure type, can be "Gram", "ExtGram", or "Grad"

Notes

• The system has M+1 variables in total
• The parameter χ is calculated differently depending on whether M is even or odd

Initial conditions for the DVM equations in 1D velocity space.
Setups up a Riemann problem with left and right states.

Initial conditions for Landau damping problem in Vlasov-Poisson system

The probability distribution function is a small perturbtation of a Maxwellian:
f(x,c) = 1 / √(2π) * (1 + ϵ * cos(k * x)) * exp(- c^2 / 2)

Shock tube initial conditions with left and right distribution functions

Shock tube initial conditions with left and right distribution functions in 1D3D

Initial conditions for the two-stream instability problem in Vlasov-Poisson system

The probability distribution function is given by:
    f(x,c) = 1 / √(2π) * (1 + ϵ cos(k * x)) * (exp(- c^2 / 2) * c^2

Maxwellian distribution function in 1D velocity space

Maxwellian distribution function in 3D velocity space (1D spatial)

Callback struct that triangulates and writes the 2D tree mesh and the solution variables. Can also be used for 1D Tree meshes, in which case a continuous piecewise linear solution is constructed.

#Configuration:

-solution_variables: The set of solution variables to write (via Trixi solution variables function, such as 'cons2prim') -interval: Save every interval time-steps. An interval <= 0 means do not save interval-based.

-time_interval: Save at every simulation time multiple of the time interval. A time interval <= 0 disables this functionality. Cannot be used together with interval based saving. The callback guarantees that the interval multiples are hit by reducing the time-step size if needed.

-save_initial_solution, save_final_solution, output_directory

-file_format: supported file-formats are "dat", "szplt", "vtu" and "h5". '.dat' is the Tecplot Tabular data file written in ascii, readable by both Tecplot and Paraview. '.szplt' is the Tecplot proprietary format only readable by Tecplot. Use of szplt writing requires the LDLIBRARYPATH to contain the libtecio.so shared library file. Alternatively libtecio.so can also be present in the working directory. '.vtu' is the native Paraview format for an unstructured grid. '.h5' is a dump of the relevant data matrices into a hdf5 file, for more direct reading with e.g. Mathematica. '.tsv' 1D only file format containing the vertex data 'trixi' use Trixi.SaveSolutionCallback to generate the output, instead of this triangulation callback.

-append_solution: If true, write all solution data into one file. Only the dat, szplt, h5 and tsv file formats support the append option. vtu will simply ignore the append flag. Note that Tecplot only recognizes time evolutions when they are in one file, while Paraview is the exact opposite.

-verbose: print a short message to console every time output is written

-clear_out_dir: delete the contents of the output directory at initialization

computes c^n = v^n + 0 + (n over 2) v^(n-2)C^2 + ... + (n over 1) vC^(n-1) + C^n, where u[k] is the primitive moment arising from C^(k-1)

Creates set of callbacks

# Arguments
- `semi`: Semidiscretization (Trixi.jl)
- `tspan`: Time span of the simulation
- `basis`: DG basis
- `cfl=0.45`: Maximum CFL number
- `plot_interval=20`: Plot every `plot_interval` steps
- `time_interval=20`: Save solution `time_interval` often times
- `name="gram_solution"`: Name of the output files (*.tsv)

closure(u::AbstractVector, equations::GramianMomentEquations)

Given convective moments u[1..N] (N = M+1), compute the closure u_{N+1}
using Grad's closure.

Procedure:
    - Construct Grad's distribution function with unknowns α_k:
        f_G(c) = f_M(c; ρ,v,θ) * (1 + Σ_{k=0..M} α_k c^k)
    - Compute α by solving the linear system arising from moment matching:
        ∫ f_G(c) c^i dc = u_i, i=0..M
    - Compute closure moment:
        u_{M+1}^G = ∫ f_G(c) c^{M+1} dc

closure(u, equations::GramianMomentEquations, ::Val{:ExtGramEven})

Extended Gramian closure for the even case

closure(u, equations::GramianMomentEquations, ::Val{:GramOdd})

Extended Gramian closure for the odd case

closure(u, equations::GramianMomentEquations, ::Val{:GramEven})

Gramian closure for the even case

closure(u, equations::GramianMomentEquations, ::Val{:GramOdd})

Gramian closure for the odd case

closure(u, equations::GramianMomentEquations)

Decides which closure implementation is used based on the equations.closure value

Closure transformation for Gramian moment equations in 1D with 3D velocity.

ToDo: Make generic for arbitrary M

# Arguments
- `u`: Input moments in slab geometry.
- `equations::GramianMomentEquations1D3D`: The equations struct containing constants and closure type.

# Returns
- `SVector{4, Float64}`: Transformed moments after applying closure.

Compiles the function to get the transformed higher order moments for given degree.

Plots the total variation in space of density, momentum and energy over time.

# Arguments:
- `sol`: solution object from Trixi.jl
- `M::Int`: maximum moment degree
- `semi`: semi-discretization object from Trixi.jl

# Returns:
- `p`: Plots.jl plot object containing the total variation plots
- `mass`: array of total mass over time
- `momentum`: array of total momentum over time
- `energy`: array of total energy over time

Helper function to construct a block diagonal matrix from a list of matrices.

Compute the integral using the coordinate transform ξ = C/sqrt(2*f.θ) and divide out the exp(-c^2) contained in the gauss hermite quadrature rule (note: quite the inefficient implementation here)

convective_moments_1D3D(M; rho=1.0, u=(0,0,0), theta=1.0, q=8)

Compute the velocity moments up to `M` (all multi-indices with i+j+k <= M)
using tensor-product Gauss–Hermite with `q` nodes per dimension. Returns a Vector{Float64}
with the same ordering as `multi_index_list(M)`.

Derivative of the closure function w.r.t. the moments u

Derivative of the closure function w.r.t. the moments u

Helper for double factorial (n-1)!!

(n-1)!! = (n-1) * (n-3) * (n-5) * ... until 1 or 2

Jacobian of the flux function

Jacobian of the flux function

Compute the flux Jacobian for the DVM equations in 1D velocity space.

get_moment_val(u, i, j, k, indices_list)

Helper to retrieve U{ijk} from the state vector u based on the provided indiceslist. Returns 0.0 if the moment is not in the system.

Constructs the rotation matrix for given angles theta and phi.

Get the valid indices for slab geometry up to moment degree M (which are ≠ 0)

gramian(u, n::Int)

Gramian matrix
G_{ij} = u_{i+j}

Generates all sorted tuples of length 3 with elements in {1,2,3}

Gives the indices of the higher dimensional moments
u_xx, u_xy, u_xz, u_yy, u_yz, u_zz

index3D(total_degree)

Return a vector of multi-indices (i,j,k) with i+j+k == total_degree.
Ordering matches a lexicographic-style ordering and is intended for moment lists.

Positions of moments in slab geometry which are ≠ 0

Initialize constants for Gramian moment equations in 1D with 3D velocity.

# Arguments
- `M::Int`: Number of moments.
- `angles::Vector{Tuple{Float64,Float64}}`: List of angles (θ, φ) for quadrature.

# Returns
- `Constants`: Struct containing precomputed constants.

Generates the full list of multi-indices up to nmax

convert conservative to primitive variables

convert primitive [1, v, C^2, ...] variables to conservative [1, c, c^2, ....] variables

multi_index_list(M)

Return a flattened list of multi-indices for all shells 0..M.

Index of moments necessary for the closure (u, u_x, u_xx, u_xxx, …)

Plots the maximum eigenvalue (wave-speed) of the flux Jacobian over time.

Careful: Only works for polydeg=1 (linear basis functions)!

# Arguments:
- `semi`: semi-discretization object from Trixi.jl
- `sol`: solution object from Trixi.jl
- `M`: maximum moment degree
- `n_plots`: number of time levels to plot
- `x_lower`: lower bound of the spatial domain
- `x_upper`: upper bound of the spatial domain

# Returns:
- `p`: Plots.jl plot object containing the maximum eigenvalue plots

Calculates discretization, density, velocity and pressure

Careful: Only works for polydeg=1 (linear basis functions)!

# Arguments:
- `sol`: solution object from Trixi.jl
- `M`: maximum moment degree
- `x_lower`: lower bound of the spatial domain
- `x_upper`: upper bound of the spatial domain

# Returns:
- `x`: spatial discretization points
- `ρ`: density values at discretization points
- `v`: velocity values at discretization points
- `p`: pressure values at discretization points
- `p_plot`: Plots.jl plot object containing the plots of ρ, v, and p

Plot physical variables (density, velocity, pressure) from DVM solution.

# Arguments:
- `ρ::Vector{Float64}`: Density array.
- `v::Vector{Float64}`: Velocity array.
- `p::Vector{Float64}`: Pressure array.
- `x::Vector{Float64}`: Spatial coordinates array.
- `xlims::Tuple{Float64, Float64}=(-2.0, 2.0)`: Tuple specifying x-axis limits for the plot.

# Returns:
- `plt`: Plot object containing the plotted variables.

Compute the convective moments ∫ C^{n-1} f dc where C = c - v
and return the primitive moments [ρ, v, C^2, C^3, ..., C^{N-1}]

Reads the final time level solution from a .tvd file.

# Arguments:
- `filename`: path to the .tvd file

# Returns:
- `x`: Vector of spatial coordinates
- `u_solutions`: Matrix of solution variables (variables × npts)

Reads solution data from a .tvd file (output).

# Arguments:
- `filename`: path to the .tvd file

# Returns:
- `blocks`: Vector of matrices, each matrix corresponds to a time level

relaxation_source(u, x, t, equations::GramianMomentEquations1D3D)

RHS = -1/Kn (u - u_eq)

relaxation_source(u, x, t, equations::GramianMomentEquations1D{Mp1}) where {Mp1}

RHS = 1/Kn (u - u_eq)

Compute the DVM relaxation source term with conservation constraints.

Relaxation is: (-1/Kn) * (f - f̃)

where f̃ is the corrected Maxwellian distribution ensuring conservation of mass, momentum, and energy.
ρ = ∫ f dc
ρ v = ∫ c f dc
ρ θ = ∫ (c - v)² f dc

Steps:
1. Compute macroscopic moments ρ, v, θ from f.
2. Compute the Maxwellian distribution f_Maxwellian based on these moments.
3. Set up the Lagrange-multiplier system to enforce conservation of mass, momentum, and energy.
4. Solve for the correction Δf to f_Maxwellian.
5. Return the relaxation source term.

Constructs the full rotation matrix for moments up to given degree.

Helper function to set the 1D DVM equations with shock tube initial conditions up.

# Arguments
- `N`: Number of discrete velocity points.
- `Kn`: Knudsen number.
- `c_l`: Lower bound of the velocity domain.
- `c_u`: Upper bound of the velocity domain.
- `ic_left`: Function defining the left initial condition.
- `ic_right`: Function defining the right initial condition.
- `base_tree_level=8`: Base refinement level for the mesh.
- `surface_flux`=flux_lax_friedrichs: Numerical flux function for the surface integral.
- `volume_flux`=flux_central: Numerical flux function for the volume integral.

# Returns
- `basis`: DG basis functions.
- `mesh`: Computational mesh.
- `equations`: DVM equations object.
- `initial_condition`: Initial conditions object.
- `solver`: DGSEM solver object.
- `boundary_conditions`: Boundary conditions object.

Setup ready to create semidiscretizations of the Gramian moment equations in 1D

# Arguments
- `M`: Number of moments
- `Kn`: Knudsen number
- `closure`: Closure function
- `f_left`: Left distribution function for Riemann initial condition
- `f_right`: Right distribution function for Riemann initial condition
- `base_tree_level=8`: Base tree level for the mesh
- `surface_flux=flux_lax_friedrichs`: Surface flux function
- `volume_flux=flux_central`: Volume flux function
- `polydeg=1`: Polynomial degree for DG basis
- `domain=(-2.0, 2.0)`: Spatial domain
- `χ_set="optimal"`: Set of χ values for closure - only relevant for Extended Gramian Closure with even M
- `alpha_max=0.5`: Maximum shock capturing parameter
- `alpha_min=0.001`: Minimum shock capturing parameter
- `alpha_smooth=true`: Smooth shock capturing parameter

solve_alpha(u::AbstractVector, ρ::Real, v::Real, θ::Real; λ::Float64=0.0)

Given convective moments u[1..N], compute the Grad coefficients α_k by solving the linear system arising from moment matching:
    ∫ f_G(c) c^i dc = u_i, i=0..M
where f_G(c) = f_M(c; ρ,v,θ) * (1 + Σ_{k=0..M} α_k c^k)

Solve the Poisson equation ∂E/∂x = ρ - ⟨ρ⟩ with periodic boundary conditions using FFT

Constructs the tensor transformation matrix for rank-n tensors
Combines rotation matrices for each index of the tensor.

Callback function to solve the Poisson equation and store electric field data

Create a DiscreteCallback for the Vlasov-Poisson system to solve Poisson equation at each timestep

Source term that applies the source term from the electric field to the moment equations

Zero source, no RHS

Extract physical variables from the DVM solution.

# Arguments:
- `semi`: Semi-discretization object (Trixi.jl).
- `sol`: Solution object containing the solution at different time steps (Trixi.jl).
- `equations::DVMEquations1D`: DVM equations object (DVMEquations1D).

# Returns:
- `x::Vector{Float64}`: Spatial coordinates.
- `ρ::Vector{Float64}`: Density at each spatial coordinate.
- `v::Vector{Float64}`: Velocity at each spatial coordinate.
- `Θ::Vector{Float64}`: Temperature at each spatial coordinate.
- `p::Vector{Float64}`: Pressure at each spatial coordinate.

Trixi.flux(u, orientation::Integer, equations::GramianMomentEquations1D3D{Mp1}) where {Mp1}
Flux function: F(U) = u_{k+1}, k=0,…,M where u_{M+1} = C(u0, u1, …, uM)

Trixi.flux(u, orientation::Integer, equations::GramianMomentEquations1D{Mp1}) where {Mp1}

Flux function: F(U) = u_{k+1}, k=0,…,M where u_{M+1} = C(u0, u1, …, uM)

Compute the flux for the DVM equations in 1D velocity space.

f_i * c_i for i = 1,...,N

Calculate maximum wave speed for local Lax-Friedrichs-type dissipation

λ_max = v_x ± γ √(θ), with γ=5

Calculate maximum wave speed for local Lax-Friedrichs-type dissipation

Estimate the flux Jacobian eigenvalues by means of Gerschgorin

Estimate the flux Jacobian eigenvalues by means of Gerschgorin