API Reference

HyQMOM

3D Hyperbolic Quadrature Method of Moments (HyQMOM) solver with MPI parallelization.

This package implements a moment-based kinetic solver for the Boltzmann equation with BGK collision operator, using domain decomposition for parallel execution.

run_simulation(; Np=nothing, Nx=20, Ny=20, tmax=0.1, num_workers=1, kwargs...)

High-level wrapper for running HyQMOM simulations with sensible defaults.

This function provides a convenient interface to the HyQMOM solver with automatic MPI initialization, parameter setup, and optional visualization. It handles both serial and parallel execution transparently.

Arguments

• Np: (Legacy) Global grid size for square grid (Np x Np), overrides Nx/Ny if provided
• Nx: Global grid size in x direction, default 20
• Ny: Global grid size in y direction, default 20
• Nz: Global grid size in z direction, default 1
• tmax: Maximum simulation time, default 0.1
• num_workers: Number of MPI ranks (must match mpiexec -n), default 1
• verbose: Print progress information, default true
• save_output: Save output to file, default false
• enable_plots: Generate PyPlot visualization after simulation, default false
• save_figures: Save figures to disk (requires enable_plots=true), default false
• output_dir: Directory for saved figures, default "."
• Kn: Knudsen number (controls collision frequency), default 1.0
• Ma: Mach number (controls jet velocity), default 0.0
• flag2D: 2D flag (1 for 2D, 0 for 3D), default 0
• CFL: CFL number for numerical stability, default 0.5
• homogeneous_z: Jets at all z levels (true) or only lower half (false), default true

Returns

Dictionary with simulation results:

• :M: Final moment field (Nx x Ny x Nz x 35) on rank 0, nothing on other ranks
• :final_time: Actual final time reached
• :time_steps: Number of time steps taken
• :xm, :ym, :zm: Grid coordinates (only on rank 0)
• :Nx, :Ny, :Nz: Grid sizes
• :tmax: Requested max time

Examples

# Single rank with square grid
results = run_simulation(Nx=20, Ny=20, tmax=0.1)

# Multi-rank (run with: mpiexec -n 4 julia script.jl)
results = run_simulation(Nx=40, Ny=40, tmax=0.2, num_workers=4)

# With visualization and non-square grid
results = run_simulation(Nx=60, Ny=40, tmax=0.1, enable_plots=true)

# Legacy: Use Np for square grids
results = run_simulation(Np=40, tmax=0.1)

See Also

• simulation_runner: Lower-level simulation function
• run_simulation_with_snapshots: Simulation with time-series output

simulation_runner(params)

Run MPI-parallel simulation loop.

This is the core time-stepping loop that orchestrates all components:

• Domain decomposition
• Initial conditions
• Time evolution with flux computation
• Halo exchange
• Realizability enforcement
• BGK collision
• Result gathering

Arguments

• params: Named tuple with simulation parameters:
  • Nx: Global grid size in x direction
  • Ny: Global grid size in y direction
  • Nz: Global grid size in z direction
  • tmax: Maximum simulation time
  • Kn: Knudsen number
  • Ma: Mach number
  • flag2D: 2D flag (1 for 2D, 0 for 3D)
  • CFL: CFL number
  • dx, dy, dz: Grid spacing
  • Nmom: Number of moments (35)
  • nnmax: Maximum number of time steps
  • dtmax: Maximum time step size
  • IC parameters: rhol, rhor, T, r110, r101, r011
  • symmetry_check_interval: How often to check symmetry
  • homogeneous_z: Whether jets exist at all z levels (validation mode)
  • enable_memory_tracking: Enable memory tracking (not implemented)
  • snapshot_interval: (optional) Save snapshots every N steps. If not provided or 0, no snapshots saved.
  • Standardized moments (S4) and central moments (C4) are automatically saved with each snapshot

Returns

If snapshot_interval is provided and > 0:

• Snapshots are streamed directly to a JLD2 file (one snapshot at a time, never all in memory)
• On rank 0: snapshot_filename (path to JLD2 file), grid_out
• On other ranks: nothing, nothing

Snapshot file structure:

• meta/params: Simulation parameters
• meta/snapshot_interval: Interval value
• meta/n_snapshots: Total number of snapshots saved
• grid: Grid structure
• snapshots/000001/{M, t, step[, S][, C]}: First snapshot
• snapshots/000002/...: Second snapshot, etc.

Otherwise (standard mode):

• M_final: Final moment array (global, gathered to rank 0), or nothing on other ranks
• final_time: Final simulation time
• time_steps: Number of time steps taken
• grid_out: Grid structure (only on rank 0)

Algorithm

See simulation_runner.m for detailed algorithm description.

run_simulation_with_snapshots(params; snapshot_interval=2)

Run simulation with time-series snapshot collection.

This function extends simulation_runner to automatically collect and save simulation snapshots at regular intervals. Snapshots are streamed to a JLD2 file for memory efficiency and can be visualized with interactive_3d_timeseries_streaming.

Arguments

• params: Simulation parameters (same as simulation_runner)
• snapshot_interval: Save snapshots every N time steps (default: 2)

Returns

• filename: Path to JLD2 snapshot file (rank 0 only)
• grid: Grid structure for visualization (rank 0 only)

On non-zero ranks, returns (nothing, nothing).

Examples

using HyQMOM, MPI

MPI.Init()

# Set up parameters
params = (
    Nx = 40, Ny = 40, Nz = 40,
    tmax = 0.1, Ma = 1.0, Kn = 1.0,
    CFL = 0.7, flag2D = 0,
    # ... other parameters
)

# Run with snapshot collection
filename, grid = run_simulation_with_snapshots(params; snapshot_interval=5)

# Visualize results (rank 0 only)
if MPI.Comm_rank(MPI.COMM_WORLD) == 0 && filename !== nothing
    interactive_3d_timeseries_streaming(filename, grid, params)
end

MPI.Finalize()

See Also

• simulation_runner: Core simulation function
• interactive_3d_timeseries_streaming: Visualization of snapshot files (requires GLMakie)

CubicRegion

Defines a cubic/box region with specified moment properties.

Fields

• center::NTuple{3,Float64}: Center position (x, y, z) in physical coordinates
• width::NTuple{3,Float64}: Width in each direction (dx, dy, dz)
• density::Float64: Density (rho) in the region
• velocity::NTuple{3,Float64}: Mean velocity (U, V, W)
• temperature::Float64: Temperature (for covariance matrix)

Example

# Cube centered at origin, 0.2 units wide, moving in +x direction
cube = CubicRegion(
    center = (0.0, 0.0, 0.0),
    width = (0.2, 0.2, 0.2),
    density = 1.0,
    velocity = (0.5, 0.0, 0.0),
    temperature = 1.0
)

initialize_moment_field(grid_params, background, regions; 
                       r110=0.0, r101=0.0, r011=0.0)

Initialize a 3D moment field with background state and cubic regions.

Arguments

• grid_params: NamedTuple with (Nx, Ny, Nz, xmin, xmax, ymin, ymax, zmin, zmax, xm, ym, zm)
• background: CubicRegion defining the background state (uniform)
• regions: Vector of CubicRegion objects to place in the domain
• r110, r101, r011: Correlation coefficients for covariance matrix

Returns

• M: 4D array (nx, ny, nz, Nmom) of initialized moments

Example

# Background at rest with low density
background = CubicRegion(
    center = (0.0, 0.0, 0.0),
    width = (Inf, Inf, Inf),  # Fills entire domain
    density = 0.01,
    velocity = (0.0, 0.0, 0.0),
    temperature = 1.0
)

# Two crossing jets
jet1 = CubicRegion(
    center = (-0.25, -0.25, 0.0),
    width = (0.1, 0.1, 0.1),
    density = 1.0,
    velocity = (0.707, 0.707, 0.0),  # Moving northeast
    temperature = 1.0
)

jet2 = CubicRegion(
    center = (0.25, 0.25, 0.0),
    width = (0.1, 0.1, 0.1),
    density = 1.0,
    velocity = (-0.707, -0.707, 0.0),  # Moving southwest
    temperature = 1.0
)

M = initialize_moment_field(grid_params, background, [jet1, jet2])

initialize_moment_field_mpi(decomp, grid_params, background, regions; kwargs...)

MPI-aware version that only initializes local subdomain.

Arguments

• decomp: MPI decomposition structure
• grid_params: Grid parameters
• background: Background CubicRegion
• regions: Vector of CubicRegion objects
• halo: Halo size (default: 2)

Returns

• M: 4D array with halos (nx+2halo, ny+2halo, nz, Nmom)

crossing_jets_ic(Nx, Ny, Nz, xmin, xmax, ymin, ymax, zmin, zmax;
                Ma=0.0, rhol=1.0, rhor=0.01, T=1.0, 
                jet_size=0.1, jet_offset=0.0)

Create standard crossing jets initial condition using the flexible system.

This is a convenience function that creates the classic two-jet configuration.

Arguments

• Grid parameters: Nx, Ny, Nz, domain bounds
• Ma: Mach number (determines jet velocity magnitude)
• rhol: Jet density
• rhor: Background density
• T: Temperature
• jet_size: Fraction of domain for jet width (default: 0.1 = 10%)
• jet_offset: Offset from center in units of jet_size (default: 0)

Returns

• background::CubicRegion
• regions::Vector{CubicRegion} with two jets

InitializeM4_35(M000, umean, vmean, wmean, C200, C110, C101, C020, C011, C002)

Compute 3D fourth-order joint Gaussian moments from physical parameters.

Arguments

• M000: Number density
• umean, vmean, wmean: Mean velocities (M100/M000, M010/M000, M001/M000)
• C200, C110, C101, C020, C011, C002: Covariance matrix elements

Returns

• M: 35-element vector of raw moments up to 4th order

M2CS4_35(M4)

Compute central and standardized moments from raw moments.

Converts 35 raw moments to central (C4) and standardized (S4) moments.

Arguments

• M4: 35-element vector of raw moments

Returns

• C4: 35-element vector of central moments
• S4: 35-element vector of standardized moments

Moments5_3D(M4)

Compute 3D HyQMOM closure for fluxes.

Computes 5th-order closure moments from 4th-order moments using the HyQMOM method.

Arguments

• M4: 35-element vector of 4th-order raw moments

Returns

• M5: 21-element vector of 5th-order raw moments
• C5: 21-element vector of 5th-order central moments
• S5: 21-element vector of 5th-order standardized moments

hyqmom_3D(S300, S400, S110, S210, S310, S120, S220, S030, S130, S040,
          S101, S201, S301, S102, S202, S003, S103, S004, S011, S111,
          S211, S021, S121, S031, S012, S112, S013, S022)

Find closures for 5th-order standardized moments in 3D using HyQMOM.

Returns 21 5th-order standardized moments computed from 4th-order moments.

realizability(operation::Symbol, args...)

Dispatcher for realizability operations.

Operations

• :2D - Find realizable moment set in 2D (5 outputs)
• :3D - Check and correct realizability of cross moments in 3D (29 outputs)
• :S111 - Check and correct realizability of S111 (1 output)
• :S2 - Check and correct realizability of 2nd-order moments (4 outputs)
• :S210 - Check and correct realizability of S210, S201 (2 outputs)
• :S211 - Check and correct realizability of S211 (1 output)
• :S310 - Check and correct realizability of S310 and S220 (2 outputs)
• :S310_220 - Check and correct realizability of S220 for S310 (1 output)
• :S220 - Check maximum bounds and correct S220 (1 output)

Examples

# 2D realizability
S21, S12, S31, S22, S13 = realizability(:2D, S30, S40, S11, S21, S31, S12, S22, S03, S13, S04)

# 3D realizability
S300, ..., flag220 = realizability(:3D, S300, ..., S022)

# S111 correction
S111r = realizability(:S111, S110, S101, S011, S210, S201, S120, S021, S102, S012, S111)

# S2 correction
S110r, S101r, S011r, S2r = realizability(:S2, S110, S101, S011)

realizable_2D(S30, S40, S11, S21, S31, S12, S22, S03, S13, S04)

Find realizable moment set in 2D.

Arguments

• S30, S40: Third and fourth-order moments along first axis
• S11: Second-order cross-moment
• S21, S31: Third-order cross-moments
• S12, S22: Cross-moments
• S03, S13, S04: Moments along second axis

Returns

• S21, S12, S31, S22, S13: Corrected realizable moments

Algorithm

Sequentially enforces realizability constraints:

1. Check realizability of S11 (first minor)
2. Check realizability of S12 and S21
3. Check realizability of S22 (second minor)
4. Check realizability of S13 and S31 given S22
5. Check third minor using Cholesky determinant and polynomial roots

realizable_3D(S300, S400, S110, S210, S310, S120, S220, S030, S130, S040,
              S101, S201, S301, S102, S202, S003, S103, S004, S011, S111,
              S211, S021, S121, S031, S012, S112, S013, S022)

Check and correct realizability of cross moments in 3D.

This is the most complex realizability function, handling 3D cross-moment constraints. It enforces positive-definiteness of the 6x6 realizability matrix.

Arguments

• 28 standardized moments from 3rd and 4th order

Returns

• 28 corrected moments + flag220 indicator

Algorithm

1. Compute H200, H020, H002 (variance-related quantities)
2. Check maximum bounds on S220, S202, S022
3. Check and correct S110, S101, S011 (2nd-order)
4. Handle degenerate cases (faces) when S2 ~= 0
5. Handle interior cases:
   • Check diagonal elements of delta2star3D matrix
   • Correct 3rd-order moments (S210, S201, S120, S021, S102, S012)
   • Correct S111
   • Apply S310_220 constraints
   • Compute lower bounds for S220, S202, S022
   • Recheck diagonal elements
   • Compute bound values for off-diagonal moments using bound_minor1
   • Apply corrections based on which diagonal elements are problematic
6. Iteratively refine until all constraints are satisfied

Notes

• This is a 830-line function in MATLAB - the full port requires careful validation against test cases
• The current implementation provides the structure; full details need to be filled in by carefully translating each section of the MATLAB code

Flux_closure35_and_realizable_3D(M4, flag2D, Ma)

Compute 3D fluxes for all moments and correct unrealizable moments.

This is the main closure function that orchestrates the entire moment pipeline:

1. Convert raw moments to central and standardized moments
2. Enforce univariate realizability
3. Check and correct 2D realizability in all planes
4. Handle edge/corner cases
5. Apply 3D realizability
6. Compute 5th-order closure via HyQMOM
7. Convert back to raw moments
8. Assemble flux vectors

Arguments

• M4: 35-element moment vector up to 4th order [M000,M100,M200,M300,M400,M010,M110,M210,M310,M020,M120,M220,M030,M130,M040, M001,M101,M201,M301,M002,M102,M202,M003,M103,M004,M011,M111,M211,M021,M121, M031,M012,M112,M013,M022]
• flag2D: 2D simulation flag (1 for 2D, 0 for 3D)
• Ma: Mach number

Returns

• Fx: X-direction flux moments (35 elements, 5th order)
• Fy: Y-direction flux moments (35 elements, 5th order)
• Fz: Z-direction flux moments (35 elements, 5th order)
• M4r: Realizable 4th-order moments (35 elements)

Algorithm

See Fluxclosure35andrealizable3D.m for detailed algorithm description.

flux_HLL(Wstar, W, l1, l2, F, N)

Compute HLL flux from intermediate states and wave speeds.

This is a simple 5-line function that computes the numerical flux using the HLL (Harten-Lax-van Leer) scheme.

Arguments

• Wstar: Intermediate state array
• W: State array
• l1: Left wave speeds
• l2: Right wave speeds
• F: Flux array
• N: Number of cells

Returns

• Flux: Numerical flux array

pas_HLL(M, F, dt, dx, vpmin, vpmax; apply_bc_left=true, apply_bc_right=true)

HLL flux update scheme.

Computes the updated moments using the HLL (Harten-Lax-van Leer) flux scheme.

Arguments

• M: Moment array (Np x Nmom)
• F: Flux array (Np x Nmom)
• dt: Time step
• dx: Spatial step
• vpmin: Minimum eigenvalues (Np)
• vpmax: Maximum eigenvalues (Np)
• apply_bc_left: Apply boundary conditions at left boundary (default: true)
• apply_bc_right: Apply boundary conditions at right boundary (default: true)

Returns

• Mp: Updated moment array

Notes

• Set apply_bc_left=false or apply_bc_right=false for processor boundaries in MPI

closure_and_eigenvalues(mom)

Compute moment closure and eigenvalue bounds using Chebyshev algorithm.

This function computes the (2N+1)-th order moment from moments of order 0 to 2N using the Chebyshev algorithm, and also computes the min/max eigenvalues.

Arguments

• mom: Vector of moments of order 0 to 2N (length 2N+1)

Returns

• Mp: Moment of order 2N+1
• vpmin: Minimum eigenvalue
• vpmax: Maximum eigenvalue

eigenvalues6_hyperbolic_3D(M, axis, flag2D, Ma)

Compute unified eigenvalues of 3D flux Jacobian.

Arguments

• M: 35-element moment vector
• axis: Direction (1=X direction, 2=Y direction)
• flag2D: 2D simulation flag
• Ma: Mach number

Returns

• v6min, v6max: Min/max eigenvalues for hyperbolicity
• Mr: Corrected moment vector

Algorithm

1. Compute eigenvalues from current moments
2. Check for complex eigenvalues
3. If complex, correct moments to ensure real eigenvalues
4. Recompute eigenvalues with corrected moments

eigenvalues6z_hyperbolic_3D(M::Vector{Float64}, flag2D::Int, Ma::Float64)

Eigenvalues of 3-D flux Jacobian in z direction.

Arguments

• M: 35-element moment vector
• flag2D: 2D simulation flag
• Ma: Mach number

Returns

• v6min: Minimum eigenvalue for hyperbolicity
• v6max: Maximum eigenvalue for hyperbolicity
• Mr: Corrected moment vector

Note

M = [M000,M100,M200,M300,M400,M010,M110,M210,M310,M020,M120,M220,M030,M130,M040, M001,M101,M201,M301,M002,M102,M202,M003,M103,M004,M011,M111,M211,M021,M121, M031,M012,M112,M013,M022]

collision35(M, dt, Kn)

Apply elastic BGK collision operator to moments.

Relaxes moments toward Maxwellian equilibrium using BGK model: dM/dt = (MG - M)/tc

Arguments

• M: 35-element moment vector
• dt: Time step
• Kn: Knudsen number

Returns

• Mout: Updated moments after collision

setup_mpi_cartesian_3d(Nx::Int, Ny::Int, Nz::Int, halo::Int, comm::MPI.Comm)

Create a 2D Cartesian domain decomposition over MPI ranks (no z decomposition).

Arguments

• Nx: Global grid size along x direction
• Ny: Global grid size along y direction
• Nz: Global grid size along z (no decomposition, all ranks have full z)
• halo: Halo width in cells (typically 2)
• comm: MPI communicator

Returns

A named tuple with fields:

• nx_global: Global grid size in x
• ny_global: Global grid size in y
• nz_global: Global grid size in z
• halo: Halo width
• dims: (Px, Py, 1) process grid dimensions (Pz=1, no z decomposition)
• coords: (rx, ry, 0) coordinates of this rank in process grid (0-based)
• rank: This rank index (0-based)
• neighbors: Named tuple with fields left, right, down, up (rank indices or -1)
• local_size: (nxlocal, nylocal, nz) interior size (without halos)
• istart_iend: (i0, i1) global inclusive interior index range for x
• jstart_jend: (j0, j1) global inclusive interior index range for y
• kstart_kend: (1, nz) global inclusive interior index range for z (always full)

Notes

• Uses an approximately square process grid (Px, Py, 1) with Px*Py=nprocs
• Uses block decomposition in x and y with remainder cells assigned to lower coords
• No decomposition in z: all ranks have the full z-dimension (Pz=1)

Example

decomp = setup_mpi_cartesian_3d(40, 40, 4, 2, MPI.COMM_WORLD)
# For 4 ranks: creates 2x2x1 process grid
# Each rank gets 20x20x4 interior cells

halo_exchange_3d!(A::Array{T,4}, decomp, bc::Symbol) where T

Exchange halos for a 3D subdomain with nv variables (no z decomposition).

Arguments

• A: Local array (nx+2h, ny+2h, nz, nv), interior: A[h+1:h+nx, h+1:h+ny, :, :]
• decomp: Struct from setupmpicartesian_3d
• bc: Boundary condition type (default: :copy)

Returns

• Modifies A in place with exchanged halos

Notes

• Performs left/right exchange first, then up/down exchange
• Applies physical BC at global boundaries before exchange
• No exchange in z-direction (all ranks have full z extent)
• Corners are filled implicitly after second phase

Example

A = zeros(nx+2*halo, ny+2*halo, nz, Nmom)
halo_exchange_3d!(A, decomp, :copy)

compute_halo_fluxes_and_wavespeeds_3d!(Fx, Fy, vpxmin_ext, vpxmax_ext, vpymin_ext, vpymax_ext,
                                       M, vpxmin, vpxmax, vpymin, vpymax, vpzmin, vpzmax,
                                       nx, ny, nz, halo, flag2D, Ma)

Compute fluxes and wave speeds in halo cells (3D physical space).

Description

After halo exchange, the moment data M in halo cells is available from neighbors. This function computes the corresponding fluxes and wave speeds in those halo cells, which are needed for the pas_HLL stencil at processor boundaries.

Arguments

• Fx: X-flux array with halos (nx+2halo, ny+2halo, nz, Nmom) [modified in place]
• Fy: Y-flux array with halos (nx+2halo, ny+2halo, nz, Nmom) [modified in place]
• vpxmin_ext: Extended min x-wave speed (nx+2*halo, ny, nz) [output]
• vpxmax_ext: Extended max x-wave speed (nx+2*halo, ny, nz) [output]
• vpymin_ext: Extended min y-wave speed (nx, ny+2*halo, nz) [output]
• vpymax_ext: Extended max y-wave speed (nx, ny+2*halo, nz) [output]
• M: Moment array with halos (nx+2halo, ny+2halo, nz, Nmom)
• vpxmin: Min x-wave speed, interior only (nx, ny, nz)
• vpxmax: Max x-wave speed, interior only (nx, ny, nz)
• vpymin: Min y-wave speed, interior only (nx, ny, nz)
• vpymax: Max y-wave speed, interior only (nx, ny, nz)
• vpzmin: Min z-wave speed, interior only (nx, ny, nz)
• vpzmax: Max z-wave speed, interior only (nx, ny, nz)
• nx: Interior size in x
• ny: Interior size in y
• nz: Interior size in z
• halo: Halo width
• flag2D: 2D flag for flux closure
• Ma: Mach number for flux closure

apply_flux_update_3d!(Mnp, M, F, vpmin, vpmax, vpmin_ext, vpmax_ext,
                      nx, ny, nz, halo, dt, ds, decomp, axis)

Apply flux update with processor boundary handling (3D physical space).

Arguments

• Mnp: Updated moment array (output, modified in place)
• M: Moment array with halos (nx+2halo, ny+2halo, nz, Nmom)
• F: Flux array (nx+2halo, ny+2halo, nz, Nmom)
• vpmin: Min wave speed, interior (nx, ny, nz) for all directions
• vpmax: Max wave speed, interior (nx, ny, nz) for all directions
• vpmin_ext: Min wave speed, extended (nx+2halo, ny, nz) for X, (nx, ny+2halo, nz) for Y, unused for Z
• vpmax_ext: Max wave speed, extended (nx+2halo, ny, nz) for X, (nx, ny+2halo, nz) for Y, unused for Z
• nx: Interior size in x
• ny: Interior size in y
• nz: Interior size in z
• halo: Halo width (for x-y only, no halos in z)
• dt: Time step size
• ds: Grid spacing (dx for X, dy for Y, dz for Z)
• decomp: Domain decomposition structure with neighbors field
• axis: 1 for X-direction, 2 for Y-direction, 3 for Z-direction

Description

Unified flux update for X, Y, and Z directions using pasHLL, with special handling for processor boundaries. At processor boundaries in x-y, includes one halo cell so pasHLL can see neighbor data. Z-direction is simpler since there's no MPI decomposition (Pz=1).

moment_idx(order)

Return linear indices for moment arrays.

Helper function for accessing moments in 3D arrays.

moment_idx(name)

Return the index of a moment by name (e.g., "M110" -> 7).

Arguments

• name: Moment name as a string (e.g., "M110", "M200")

Returns

• Integer index in the 35-element moment vector

M4_to_vars(M4)

Extract M4 array (35 moments) to individual variables.

Returns

35 individual moment variables in canonical order.

M4_to_vars(M4::AbstractArray{T,3})

Extract M4 3D array (5x5x5) to individual variables.

Returns

35 individual moment variables in canonical order.

M5_to_vars(M5)

Extract M5 array (56 moments) to individual variables.

Returns

56 individual moment variables in canonical order (up to 5th order).

M5_to_vars(M5::AbstractArray{T,3})

Extract M5 3D array (6x6x6) to individual variables.

Returns

56 individual moment variables in canonical order (up to 5th order).

axis_moment_slice(M, axis)

Extract 15-element moment vector for jacobian6 computation.

Arguments

• M: 35-element moment vector
• axis: Direction
  • 1 = X (UV plane)
  • 2 = Y (VU plane)
  • 3 = Z (UW plane)
  • 4 = VW plane

Returns

• M_slice: 15-element moment vector in jacobian6 order: [M000, M010, M020, M030, M040, M100, M110, M120, M130, M200, M210, M220, M300, M310, M400]

Examples

# X-direction eigenvalues (UV plane)
J6 = jacobian6(axis_moment_slice(M, 1)...)

# Y-direction eigenvalues (VU plane)
J6 = jacobian6(axis_moment_slice(M, 2)...)

Notes

This eliminates repeated index slicing logic in eigenvalue computations and plotting functions.

compute_standardized_field(M)

Compute standardized moments for entire field.

Arguments

• M: 4D array (Nx, Ny, Nz, Nmom) of raw moments

Returns

• S: 4D array (Nx, Ny, Nz, Nmom) of standardized moments

Notes

Applies M2CS4_35 to each grid point to extract standardized moments. This is useful for visualization and analysis of snapshot data.

compute_central_field(M)

Compute central moments for entire field.

Arguments

• M: 4D array (Nx, Ny, Nz, Nmom) of raw moments

Returns

• C: 4D array (Nx, Ny, Nz, Nmom) of central moments

Notes

Applies M2CS4_35 to each grid point to extract central moments.

get_standardized_moment(S, moment_name::String)

Extract a specific standardized moment from the field.

Arguments

• S: 4D array (Nx, Ny, Nz, Nmom) of standardized moments
• moment_name: Name of moment (e.g., "S110", "S101", "S022")

Returns

• 3D array (Nx, Ny, Nz) of the requested moment

Supported moments

Second order: S110, S101, S011 Third order: S300, S210, S201, S120, S111, S102, S030, S021, S012, S003 Fourth order: S400, S310, S301, S220, S211, S202, S130, S121, S112, S103, S040, S031, S022, S013, S004

get_central_moment(C, moment_name::String)

Extract a specific central moment from the field.

Arguments

• C: 4D array (Nx, Ny, Nz, Nmom) of central moments
• moment_name: Name of moment (e.g., "C110", "C101", "C022")

Returns

• 3D array (Nx, Ny, Nz) of the requested moment

Supported moments

Same indexing as standardized moments but with 'C' prefix.

delta2star3D(s300, s400, s110, s210, s310, s120, s220, s030, s130, s040,
             s101, s201, s301, s102, s202, s003, s103, s004, s011, s111,
             s211, s021, s121, s031, s012, s112, s013, s022)

Compute the 6x6 Delta2* matrix for 3D moment realizability checking.

This is autogenerated code translated from MATLAB. The computation involves ~200 lines of arithmetic operations to construct the realizability matrix.

delta2star3D_permutation(S300, S400, S110, S210, S310, S120, S220, S030, S130, S040,
                         S101, S201, S301, S102, S202, S003, S103, S004, S011, S111,
                         S211, S021, S121, S031, S012, S112, S013, S022)

Find 6 permutations of Delta2* (all must be positive for realizability).

This function computes the Delta2* matrix for all 6 permutations of the coordinate axes:

• E1: Delta2*_ijk (original)
• E2: Delta2*_ikj (swap j and k)
• E3: Delta2*_jik (swap i and j)
• E4: Delta2*_jki (cyclic permutation)
• E5: Delta2*_kij (cyclic permutation)
• E6: Delta2*_kji (swap i and k)

Arguments

All 28 standardized moments up to 4th order (excluding those determined by symmetry).

Returns

• E1, E2, E3, E4, E5, E6: Six 6x6 matrices, one for each permutation

Realizability

For moments to be realizable, all six matrices must have non-negative determinants for all principal minors.

jacobian6(m00, m01, m02, m03, m04, m10, m11, m12, m13, m20, m21, m22, m30, m31, m40)

Compute 6x6 Jacobian matrix for moment system.

This function was autogenerated by MATLAB Symbolic Math Toolbox.

Arguments

• 15 moment components (m00 through m40)

Returns

• 6x6 Jacobian matrix

M4toC4_3D(M000, M100, M200, M300, M400, M010, M110, M210, M310, M020, M120, M220, M030, M130, M040, M001, M101, M201, M301, M002, M102, M202, M003, M103, M004, M011, M111, M211, M021, M121, M031, M012, M112, M013, M022)

Convert 4th-order raw moments to central moments in 3D.

This function was autogenerated by MATLAB Symbolic Math Toolbox.

Arguments

• 35 raw moment components up to 4th order

Returns

• 5x5x5 array of central moments up to 4th order

C4toM4_3D(M000, umean, vmean, wmean, C200, C110, C101, C020, C011, C002,
          C300, C210, C201, C120, C111, C102, C030, C021, C012, C003,
          C400, C310, C301, C220, C211, C202, C130, C121, C112, C103, C040, C031, C022, C013, C004)

Convert central moments to raw moments (autogenerated from Symbolic Math Toolbox). Returns 5x5x5 array of raw moments.

S4toC4_3D_r(C200, C110, C101, C020, C011, C002, S300, S210, S201, S120, S111, S102, S030, S021, S012, S003, S400, S310, S301, S220, S211, S202, S130, S121, S112, S103, S040, S031, S022, S013, S004)

Convert standardized moments to central moments using matrix square root.

This function was autogenerated by MATLAB Symbolic Math Toolbox.

Arguments

• Second-order central moments (C200, C110, C101, C020, C011, C002)
• Standardized moments from 3rd to 4th order

Returns

• 5x5x5 array of central moments up to 4th order

Notes

• Uses matrix square root of covariance matrix
• Warns if covariance matrix is not positive definite

C5toM5_3D(M000, umean, vmean, wmean, C200, C110, C101, C020, C011, C002, C300, C210, C201, C120, C111, C102, C030, C021, C012, C003, C400, C310, C301, C220, C211, C202, C130, C121, C112, C103, C040, C031, C022, C013, C004, C500, C410, C320, C230, C140, C401, C302, C203, C104, C311, C221, C131, C212, C113, C122, C050, C041, C032, C023, C014, C005)

Convert 5th-order central moments to raw moments in 3D.

This function was autogenerated by MATLAB Symbolic Math Toolbox.

Arguments

• M000: Zeroth moment (mass)
• umean, vmean, wmean: Mean velocities
• Central moments up to 5th order (56 components)

Returns

• 6x6x6 array of raw moments up to 5th order