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.

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
)

CubicRegion(; center, width, density, velocity, temperature=1.0)

Keyword constructor for CubicRegion.

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.

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

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.

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

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.

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

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).

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

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

S_to_C_batch(S110, S101, S011, S300, ..., sC200, sC020, sC002)

Batch convert standardized moments to central moments.

Algorithm

For each moment: C_ijk = S_ijk * sC200^i * sC020^j * sC002^k

Returns

All central moments in the same order as inputs.

apply_flux_update(M, F, vpmin, vpmax, vpmin_ext, vpmax_ext, nx, ny, halo, dt, ds, decomp, axis)

Apply flux update with processor boundary handling.

Description

Unified flux update for both X and Y directions using pas_HLL, with special handling for processor boundaries. At processor boundaries, includes one halo cell so pas_HLL can see neighbor data.

Arguments

• M: Moment array with halos (nx+2halo, ny+2halo, Nmom)
• F: Flux array with halos (nx+2halo, ny+2halo, Nmom)
• vpmin: Min wave speed, interior (nx, ny) for X or (nx, ny) for Y
• vpmax: Max wave speed, interior (nx, ny) for X or (nx, ny) for Y
• vpmin_ext: Min wave speed, extended (nx+2halo, ny) for X or (nx, ny+2halo) for Y
• vpmax_ext: Max wave speed, extended (nx+2halo, ny) for X or (nx, ny+2halo) for Y
• nx: Interior size in x
• ny: Interior size in y
• halo: Halo width
• dt: Time step size
• ds: Grid spacing (dx for X, dy for Y)
• decomp: Domain decomposition structure with neighbors field
• axis: 1 for X-direction, 2 for Y-direction

Returns

• Mnp: Updated moment array after flux update (nx+2halo, ny+2halo, Nmom)

Algorithm

For each line perpendicular to the flux direction:

1. Determine extent based on processor boundaries
2. Extract moments, fluxes, and wave speeds
3. Call pas_HLL with appropriate BC flags
4. Write back interior portion

Notes

The key complexity is handling processor boundaries correctly:

• Interior boundaries: Include one halo cell for pas_HLL
• Physical boundaries: Apply boundary conditions

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).

apply_physical_bc_2d!(A::Array{T,3}, decomp, bc) where T

Apply physical boundary conditions at global boundaries.

Arguments

• A: Local array with halos
• decomp: Decomposition structure
• bc: Boundary condition type (:copy, :periodic, etc.)

Algorithm

Only applies BCs at physical (global) boundaries, not at processor boundaries. Uses copy BC: halo cells copy from nearest interior cell.

Notes

This is called before halo exchange to set physical boundary values.

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

Fill halos at global boundaries based on bc type.

Arguments

• A: Local array (nx+2h, ny+2h, nz, nv)
• decomp: Domain decomposition struct
• bc: Boundary condition type

Supported bc types

• :copy - Neumann-like (copy nearest interior cell)

Notes

• Also applies copy BC in z-direction at global z boundaries (no decomposition in z)

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.

block_partition_1d(n::Int, P::Int, r::Int)

Compute 1D block decomposition for rank r (0-indexed).

Divides n points among P processes as evenly as possible. First (n mod P) ranks get (floor(n/P) + 1) points each.

Arguments

• n: Total number of points
• P: Number of processes
• r: Rank (0-indexed)

Returns

• n_local: Number of points for this rank
• i0: Starting global index (1-indexed)
• i1: Ending global index (1-indexed)

Examples

n_local, i0, i1 = block_partition_1d(100, 4, 0)  # First rank
# Returns (25, 1, 25)

bound_minor1(S003, S011, S012, S021, S030, S101, S102, S110, S111, S120, S201, S210, S300)

Compute bounds for minor moments based on standardized moments.

This function was autogenerated by MATLAB Symbolic Math Toolbox.

Arguments

• 13 standardized moment components

Returns

• 6x2 matrix of moment bounds

cell_overlaps_cube(cell_center::NTuple{3,Float64}, cell_size::NTuple{3,Float64}, cube::CubicRegion)

Check if a grid cell (defined by center and size) overlaps with a cubic region.

More robust than pointincube() for detecting when narrow jets intersect coarse grid cells. A cell and cube overlap if they intersect in all three dimensions.

Arguments

• cell_center: (x, y, z) coordinates of cell center
• cell_size: (dx, dy, dz) cell dimensions
• cube: CubicRegion to check overlap with

Returns

• true if cell and cube overlap in all three dimensions

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

Check and correct 2D moments in 3D code.

Algorithm

If |S11| >= 1, collapses to boundary case. Otherwise, applies realizability(:2D, ...) to correct cross moments.

Returns

10 corrected moments

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

Apply check2D realizability to all three coordinate planes.

Returns

30 corrected moments (10 per plane: XY, XZ, YZ)

Algorithm

Applies check2D to each of the three coordinate planes independently.

choose_process_grid(nl::Int)

Choose nearly-square factorization Px*Py=nl where Px and Py are as close as possible to sqrt(nl).

Arguments

• nl: Number of processes

Returns

• Px, Py: Process grid dimensions

Examples

Px, Py = choose_process_grid(6)  # Returns (2, 3) or (3, 2)

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

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

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.

compute_eigenvalues_for_axis(M, axis)

Helper: Compute eigenvalues for given axis.

compute_halo_fluxes_and_wavespeeds!(M, Fx, Fy, vpxmin, vpxmax, vpymin, vpymax,
                                    nx, ny, halo, flag2D, Ma)

Compute fluxes and wave speeds in halo cells.

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

• M: Moment array with halos (nx+2halo, ny+2halo, Nmom)
• Fx: X-flux array with halos (nx+2halo, ny+2halo, Nmom) - modified in-place
• Fy: Y-flux array with halos (nx+2halo, ny+2halo, Nmom) - modified in-place
• vpxmin: Min x-wave speed, interior only (nx, ny)
• vpxmax: Max x-wave speed, interior only (nx, ny)
• vpymin: Min y-wave speed, interior only (nx, ny)
• vpymax: Max y-wave speed, interior only (nx, ny)
• nx: Interior size in x
• ny: Interior size in y
• halo: Halo width
• flag2D: 2D flag for flux closure
• Ma: Mach number for flux closure

Returns

• vpxmin_ext: Extended min x-wave speed (nx+2*halo, ny)
• vpxmax_ext: Extended max x-wave speed (nx+2*halo, ny)
• vpymin_ext: Extended min y-wave speed (nx, ny+2*halo)
• vpymax_ext: Extended max y-wave speed (nx, ny+2*halo)

Algorithm

1. Create extended wave speed arrays
2. Copy interior wave speeds
3. Compute fluxes and wave speeds in left halo
4. Compute fluxes and wave speeds in right halo
5. Compute fluxes and wave speeds in bottom halo
6. Compute fluxes and wave speeds in top halo

Notes

• Modifies Fx and Fy in-place
• Returns extended wave speed arrays
• Uses closureandeigenvalues for 5-moment eigenvalues
• Uses eigenvalues6hyperbolic3D for 6-moment eigenvalues

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

compute_jacobian_eigenvalues(moments_a, moments_b)

Compute 6x6 Jacobian eigenvalues for two moment sets.

Arguments

• moments_a: 15-element vector of moments for first Jacobian
• moments_b: 15-element vector of moments for second Jacobian

Returns

• v6min, v6max: Min and max eigenvalues across both Jacobians
• lam6a, lam6b: Full eigenvalue vectors for both Jacobians

Algorithm

Constructs two 6x6 Jacobian matrices using jacobian6 and computes their eigenvalues. Returns the global min/max across both matrices.

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.

correct_moments_for_real_eigenvalues(M, axis, lam6a, lam6b, flag2D, Ma)

Helper: Correct moments to ensure real eigenvalues.

correct_uv_plane(C, S110, S300, S030, S400, S040, S220, C200, C020)

Helper: Correct UV plane moments.

correct_uw_plane(C, S101, S300, S003, S400, S004, S202, C200, C002)

Helper: Correct UW plane moments.

correct_vw_plane(C, S011, S030, S003, S040, S004, S022, C020, C002)

Helper: Correct VW plane moments.

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

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.

delta2starchol_L3(s03, s04, s11, s12, s13, s21, s22, s30, s31, s40)

Compute the L3 determinant check using Cholesky-based approach for realizability.

edge_corner_correction(R110, R101, R011, ...)

Correct moments at edges/corners of realizability domain.

Handles cases where one or more 2D correlations (R110, R101, R011) are non-realizable (<= 0), placing the state at edges or corners of the realizability domain.

Arguments

• R110, R101, R011: Realizability indicators (1 - S###^2)
• Various S###r corrected moments from check2Dallplanes

Returns

28 edge/corner-corrected standardized moments

Algorithm

Routes to appropriate edge or corner correction based on which R values are <= 0:

• Edge 1: S110 = +/-1 (xy-plane boundary)
• Edge 2: S101 = +/-1 (xz-plane boundary)
• Edge 3: S011 = +/-1 (yz-plane boundary)
• Corner: All three at boundaries

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]

enforce_univariate(S3, S4, h2min, s3max)

Enforce realizability bounds on univariate moments.

Ensures H = S4 - S3^2 - 1 >= h2min and |S3| <= s3max.

Arguments

• S3: Third standardized moment
• S4: Fourth standardized moment
• h2min: Minimum allowed H value
• s3max: Maximum allowed |S3| value

Returns

• S3: Corrected third moment
• S4: Corrected fourth moment
• H: Corrected H value

extract_all_C4_moments(C4)

Helper: Extract all C4 moments into a named tuple.

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

gather_M(M_interior::Array{T,3}, i0i1::Tuple{Int,Int}, j0j1::Tuple{Int,Int}, 
         Np::Int, Nmom::Int, comm::MPI.Comm) where T

Gather 2D moment arrays from all ranks to rank 0.

Arguments

• M_interior: Local interior moment array (nx, ny, Nmom)
• i0i1: Global index range for x (i0, i1)
• j0j1: Global index range for y (j0, j1)
• Np: Global grid size
• Nmom: Number of moments
• comm: MPI communicator

Returns

• M_full: Full moment array (only on rank 0, nothing on other ranks)

Algorithm

Rank 0 receives data from all other ranks and assembles the full array. Other ranks send their data to rank 0.

gather_M(M_interior::Array{T,4}, i0i1::Tuple{Int,Int}, j0j1::Tuple{Int,Int}, k0k1::Tuple{Int,Int},
         Np::Int, Nz::Int, Nmom::Int, comm::MPI.Comm) where T

Gather 3D moment arrays from all ranks to rank 0.

Arguments

• M_interior: Local interior moment array (nx, ny, nz, Nmom)
• i0i1: Global index range for x (i0, i1)
• j0j1: Global index range for y (j0, j1)
• k0k1: Global index range for z (k0, k1) - always (1, Nz) for 2D MPI decomposition
• Np: Global grid size in x and y
• Nz: Global grid size in z
• Nmom: Number of moments
• comm: MPI communicator

Returns

• M_full: Full moment array (only on rank 0, nothing on other ranks)

Algorithm

Rank 0 receives data from all other ranks and assembles the full array. Other ranks send their data to rank 0.

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.

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

halo_exchange_2d!(A::Array{T,3}, decomp, bc=nothing) where T

Exchange halos for a 2D subdomain with nv variables.

Arguments

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

Returns

• Updates A in-place with exchanged halos

Algorithm

1. Apply physical BC at global boundaries before exchange
2. Perform left/right exchange first
3. Then perform up/down exchange
4. Corners are filled implicitly after second phase

Notes

• Uses blocking send/receive (MPI.Send! and MPI.Recv!)
• For single-rank case, only applies physical BCs
• Physical boundaries use copy BC by default

Examples

A = zeros(nx+2*halo, ny+2*halo, Nmom)
# ... fill interior ...
halo_exchange_2d!(A, decomp)
# Now A has valid halo data

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)

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.

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)

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

lower_bound_S220(S011, S012, S021, S101, S102, S110, S111, S120, S201, S210)

Compute lower bound for S220 moment based on standardized moments.

This function was autogenerated by MATLAB Symbolic Math Toolbox.

Arguments

• 10 standardized moment components

Returns

• 3-element vector of S220 bounds

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

moment_names(order)

Return canonical moment names for given order.

Arguments

• order: Moment order (4 or 5)

Returns

• Vector of moment name strings

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

place_cubic_region!(M, region::CubicRegion, xm, ym, zm, r110, r101, r011; use_overlap=true)

Place a cubic region into an existing moment field (modifies M in-place).

Arguments

• M: Moment field array (Nx, Ny, Nz, Nmom)
• region: CubicRegion to place
• xm, ym, zm: Grid cell centers
• r110, r101, r011: Correlation coefficients
• use_overlap: If true, use cell overlap detection (more robust). If false, use point-in-cube (legacy)

point_in_cube(point::NTuple{3,Float64}, cube::CubicRegion)

Check if a point (x, y, z) is inside a cubic region.

Note: This checks if the point (typically a grid cell center) is inside the cube. For narrow jets on coarse grids, jets may be missed if their width is smaller than the grid spacing. To avoid this, ensure jet width >= grid spacing or use celloverlapscube() for more robust overlap detection.

rank_from_coords(rx::Int, ry::Int, Px::Int)

Convert (rx, ry) 0-based coordinates to 0-based rank index.

Arguments

• rx: X-coordinate in process grid (0-based)
• ry: Y-coordinate in process grid (0-based)
• Px: Number of processes in X-direction

Returns

• Rank index (0-based)

Algorithm

Uses row-major ordering: rank = ry * Px + rx

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)

realizability_S111(S110, S101, S011, S210, S201, S120, S021, S102, S012, S111)

Check and correct realizability of S111 moment.

Arguments

• Second-order standardized moments: S110, S101, S011
• Third-order standardized moments: S210, S201, S120, S021, S102, S012
• S111: Current value of S111

Returns

• S111r: Corrected value of S111

Algorithm

Computes bounds based on three different projections (A110, A101, A011) and clamps S111 to the feasible range [Rmin, Rmax].

realizability_S2(S110, S101, S011)

Check and correct realizability of 2nd-order moments.

Arguments

• S110, S101, S011: Second-order standardized cross-moments

Returns

• S110r, S101r, S011r: Corrected second-order moments
• S2r: Determinant-like quantity (should be non-negative)

Algorithm

Computes S2 = 1 + 2S110S101*S011 - (S110^2 + S101^2 + S011^2). If S2 < 0, scales all moments by a factor xr in (0,1) found via root-finding to make S2 = 0.

realizability_S210(S110, S101, S011, S300, S210, S201, H200, beta)

Check and correct realizability of S210 and S201 moments.

Arguments

• S110, S101, S011: Second-order standardized moments
• S300: Third-order moment along first axis
• S210, S201: Third-order cross-moments to correct
• H200: Variance-related quantity (H200 = max(eps, S400 - S300^2 - 1))
• beta: Scaling factor (typically 1.0)

Returns

• S210r, S201r: Corrected moments

Algorithm

Uses matrix square root and quadratic form to determine feasible region. Scales the deviation from the mean by xr in [0,1] to satisfy realizability.

realizability_S211(e11, e22, e33, e12, e13, d23, S211, beta)

Check and correct realizability of S211 moment.

Arguments

• e11, e22, e33, e12, e13: Elements from realizability matrix
• d23: Offset term
• S211: Current value of S211
• beta: Scaling factor (typically 1.0)

Returns

• S211r: Corrected value of S211

Algorithm

Computes bounds based on quadratic form involving e-matrix elements. Clamps S211 to feasible range [s211min, s211max].

realizability_S310(S110, S101, S011, S300, S210, S201, S120, S111, S310, S220, H200, beta)

Check and correct realizability of S310 and S220 moments.

Arguments

• S110, S101, S011: Second-order standardized moments
• S300, S210, S201, S120, S111: Third-order moments
• S310, S220: Fourth-order moments to correct
• H200: Variance-related quantity
• beta: Scaling factor (typically 1.0)

Returns

• S310r, S220r: Corrected moments

Algorithm

Computes bounds for S310 based on complex quadratic forms involving both 2nd and 3rd order moments. The bound depends on S220 and H200.

realizability_S310_220(S110, S101, S011, S210, S120, S111, S220)

Check and correct realizability of S220 for S310 constraint.

Arguments

• S110, S101, S011: Second-order standardized moments
• S210, S120, S111: Third-order moments
• S220: Fourth-order moment to correct

Returns

• S220r: Corrected S220 (lower bound)

Algorithm

Computes a lower bound for S220 based on the constraint that certain quadratic forms must be non-negative.

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

realizablity_S220(S110, S220, A220)

Check maximum bounds and correct S220.

Arguments

• S110: Second-order standardized moment
• S220: Fourth-order moment to correct
• A220: Bound parameter (typically sqrt((H200+S300^2)*(H020+S030^2)))

Returns

• S220r: Corrected S220

Algorithm

Clamps S220 to the range [s220min, s220max] where:

• s220min = max(S110^2, 1 - A220)
• s220max = 1 + A220

region_to_moments(region::CubicRegion, r110, r101, r011)

Convert a CubicRegion specification to a moment vector.

rootsR(Del1, H02, H20, s03, s11, s12, s13, s21, s30, s31)

Find roots of degree 3 polynomial for realizability checking. Returns a 3-element vector of roots (may be complex).

rootsR_X_Y(Y, e11, e12, e13, e23, e24, e34, e44, ex)

Find roots for realizability checking with specific matrix elements.

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

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)

send_M(M_interior::Array{T,3}, i0i1::Tuple{Int,Int}, j0j1::Tuple{Int,Int},
       dest_rank::Int, comm::MPI.Comm) where T

Send moment array to destination rank.

Arguments

• M_interior: Local interior moment array
• i0i1: Global index range for x
• j0j1: Global index range for y
• dest_rank: Destination rank
• comm: MPI communicator

Notes

This is a lower-level function. Typically use gather_M instead.

setup_mpi_cartesian_2d(np::Int, halo::Int, comm::MPI.Comm)

Create a 2D Cartesian domain decomposition.

Arguments

• np: Global grid size along x and y (np x np)
• halo: Halo width in cells (typically 1)
• comm: MPI communicator

Returns

A named tuple with fields:

• np_global: Scalar, global grid size
• halo: Scalar, halo width
• dims: (Px, Py) process grid dimensions
• coords: (rx, ry) coordinates of this rank in process grid (0-based)
• rank: This rank index (0-based)
• size: Total number of ranks
• neighbors: Named tuple with fields: left, right, down, up (rank indices or -1)
• local_size: (nxlocal, nylocal) 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

Algorithm

1. Choose approximately square process grid [Px, Py] with Px*Py = nprocs
2. Use block decomposition with remainder cells assigned to lower coordinates
3. Compute neighbor ranks using Cartesian topology

Examples

using MPI
MPI.Init()
comm = MPI.COMM_WORLD
decomp = setup_mpi_cartesian_2d(120, 1, comm)
println("Rank $(decomp.rank) has local size $(decomp.local_size)")

Notes

• Uses 0-based rank indexing (Julia convention for MPI)
• Neighbor ranks are -1 at physical boundaries
• Process grid is chosen to be as square as possible

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

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.