API Reference

ConvergenceHistory

Tracks adaptive-shift convergence information during an eigenvalue solve.

Cache of pre-discretized sparse matrix components, separated by k-power.

Domain(; x=FourierTransformed(), y=Fourier(N=60, L=1.0), z=Chebyshev(N=30, lower=0, upper=1))

A multi-dimensional domain with named coordinates. Each coordinate is one of:

• FourierTransformed() — wavenumber parameter direction
• Fourier(N=..., L=...) — resolved periodic direction
• Chebyshev(N=..., lower=..., upper=...) — resolved non-periodic direction

EVP(domain; variables, eigenvalue)

Eigenvalue problem container. Accumulates equations, BCs, parameters, and substitutions before discretization.

Usage:

prob = EVP(domain, variables=[:psi, :b], eigenvalue=:sigma)
prob[:U] = U_field
prob[:E] = 1e-12

Marker type for a Fourier-transformed (wavenumber) direction.

SolverResults

Computed eigenpairs plus metadata for one eigenproblem.

Fields

• eigenvalues::Vector{ComplexF64}
• eigenvectors::Matrix{ComplexF64}
• converged::Bool
• method_used::Symbol
• final_shift::Float64
• iterations::Int
• solve_time::Float64
• history::ConvergenceHistory

@bc left(psi) = 0
@bc right(dz(psi)) = 0
@bc left(psi) == 0        # == also works

Add a boundary condition. prob is optional.

@compute v = -dy(dz(w)) + dx(zeta)

Evaluate a DSL expression on the solved eigenvector fields and assign the result. Requires cache, eigvec, and k to be in scope (set via @compute_setup).

cache = discretize(prob)
results = solve(cache, [0.1]; sigma_0=0.02)

@compute_setup cache results[1].eigenvectors[:, 1] 0.1

@compute v = -dy(dz(w)) + dx(zeta)
@compute u = -dx(dz(w)) - dy(zeta)
@compute vort = dx(dx(psi)) + dy(dy(psi))
@compute flux = U * dy(b) - E * D2(zeta)

@compute_setup cache eigvec k

Set up the context for @compute. Call once after solving.

@compute_setup cache results[1].eigenvectors[:, 1] 0.1
@compute v = -dy(dz(w)) + dx(zeta)

@derive prob Dh2_v = rhs                              # Form 1: named operator
@derive prob v where dx(dx(v)) + dy(dy(v)) = rhs     # Form 2: inline operator

Define a derived (auxiliary) variable implicitly. The code computes Op^{-1} and substitutes v = Op^{-1} * rhs_expr into any equation that references v.

Form 1 (Op_var = rhs): The name is split on _ — everything before the last _ is a @substitution name, the part after is the variable. Dh2_v means "Dh2(v) = rhs".

Form 2 (v where lhs_expr = rhs): The operator is defined inline as any expression involving v. No separate @substitution needed.

v is NOT added to the variable list — it's eliminated automatically.

Examples:

# Form 1 (requires @substitution Dh2):
@derive prob Dh2_v = -dy(dz(w)) + dx(zeta)

# Form 2 (inline, no @substitution needed):
@derive prob v where dx(dx(v)) + dy(dy(v)) = -dy(dz(w)) + dx(zeta)

# General operator:
@derive prob v where 3*dz(dz(v)) + alpha*v = dz(w)

@derive_bc v left(v) == 0
@derive_bc v right(dz(v)) == 0
@derive_bc prob v left(v) == 0   # explicit prob

Add BCs for a derived variable's inversion. prob is optional.

@equation sigma * psi == U * dx(psi)
@equation sigma * psi = U * dx(psi)     # = also works

Add an equation LHS = RHS (or LHS == RHS). prob is optional.

@substitution Lap(A) = dx(dx(A)) + dz(dz(A))
@substitution prob Lap(A) = dx(dx(A)) + dz(dz(A))   # explicit prob

Define a substitution template. prob is optional — uses the active EVP if omitted.

assemble(cache; k_x=0.0, k_y=0.0, ...) -> (A, B)

Assemble for multiple FourierTransformed directions using keyword wavenumbers. Each keyword name must match a FourierTransformed coordinate name.

# Domain with two FourierTransformed directions:
domain = Domain(x=FourierTransformed(), y=FourierTransformed(), z=Chebyshev(N=30, lower=0, upper=1))
cache = discretize(prob)
A, B = assemble(cache; k_x=1.0, k_y=0.5)

chebyshev_evaluate(c, x) -> Vector

Evaluate a Chebyshev expansion with coefficients c at points x ∈ [-1, 1]. Uses Clenshaw's algorithm for numerical stability. The output element type is promoted from c (so complex coefficients — e.g. reconstructed eigenfunctions — give a complex result rather than throwing).

chebyshev_points(N, a=-1.0, b=1.0) -> Vector{Float64}

Compute N Chebyshev-Gauss-Lobatto points on the interval [a, b].

On [-1,1]: x_j = cos(j * pi / (N-1)) for j = 0, ..., N-1

These are then linearly mapped to [a, b].

differentiate(f, domain, coord; order=1, filter=:none) -> Vector

Spectrally differentiate a function given as physical-space values on the grid.

Chebyshev directions: converts to T coefficients, applies the differentiation operator chain D{order-1} * ⋯ * D0, converts the C^(order) result back to physical space. Chebyshev methods do not suffer from Gibbs ringing for smooth non-periodic functions.

Fourier directions: applies (ik)^order in Fourier space. An optional spectral filter can suppress Gibbs ringing from under-resolved gradients:

• :none — no filter (default; exact for smooth periodic functions)
• :exp — exponential filter exp(-α(|k|/k_max)^p) with α=36, p=36 (sharp cutoff that preserves well-resolved modes)
• :2/3 — 2/3 dealiasing rule (zeros out the top 1/3 of modes)
• A Function — custom filter σ(k, k_max) applied as multiplier

Examples

domain = Domain(z = Chebyshev(N=64, lower=0.0, upper=1.0))
z = gridpoints(domain, :z)
U = @. z^2
dUdz = differentiate(U, domain, :z)            # ≈ 2z
d2Udz2 = differentiate(U, domain, :z; order=2) # ≈ 2

domain_y = Domain(y = Fourier(N=64, L=2π))
y = gridpoints(domain_y, :y)
g = sin.(y)
dgdy = differentiate(g, domain_y, :y)                  # cos(y), no filter
dgdy = differentiate(g, domain_y, :y; filter=:exp)     # cos(y), filtered
dgdy = differentiate(g, domain_y, :y; filter=:2/3)     # cos(y), dealiased

discretize(prob::EVP; augment_derived=true) -> DiscretizationCache

Validate the problem, expand substitutions, lower derivatives, separate by k-power, and build sparse matrix components. Returns a cache for fast wavenumber assembly.

Derived variables (@derive) are augmented into the system by default: each becomes a regular unknown governed by its defining equation Op(v) = rhs, which keeps the assembled operator sparse (no dense Op⁻¹) and handles operators whose inverse the legacy path cannot form. Pass augment_derived=false to use the legacy path that eliminates derived variables via an explicit operator inverse (a smaller but denser system). The augmented form is a singular-B descriptor system; solve filters the resulting infinite/spurious modes, but shift-target the physical region.

discretize_distributed(prob; ngroups=1, kwargs...) -> DiscretizationCache

MPI-coupled discretize: builds the full cache, then returns a per-rank row-restricted cache (every rank keeps only its owned rows of the k-components; the singular-B mass filter is computed distributedly). Note: legacy derived-variable caches (augment_derived=false) are still kept full on every rank, since H(k) is rebuilt per wavenumber — so that path does not get the row-distribution memory saving. Provided by the extension BiGSTARSMPIExt (requires MPI, PetscWrap, SlepcWrap + a complex-scalar PETSc/SLEPc). Use the SAME ngroups in the subsequent solve. kwargs are forwarded to discretize (e.g. augment_derived). Run under mpiexec -n P julia ….

evaluate_field(cache, prob, eigvec, k, expr_fn) -> Vector{ComplexF64}

Evaluate a general expression on the eigenvector fields.

expr_fn is a function that takes a Dict of variable coefficient vectors and returns the result. Operators from the cache can be applied.

# Manual field computation:
v_coeffs = evaluate_field(cache, prob, eigvec, 0.1) do fields
    # fields[:w], fields[:zeta], fields[:b] are coefficient vectors
    # Apply operators manually:
    D_y = get_diff_operator(prob.domain, :y, 1)
    D_z_T = ...  # build T-basis z-derivative
    return -D_y * D_z_T * fields[:w] + im * 0.1 * fields[:zeta]
end

For derived variables, use reconstruct(cache, prob, eigvec, k, :v) instead.

fourier_points(N, L) -> Vector{Float64}

Return N equally-spaced points on [0, L), excluding the right endpoint.

meshgrid(domain, dim1, dim2) -> (Matrix, Matrix)

Return 2D meshgrid arrays for two resolved dimensions. Convention: dim2 (typically z) varies along rows, dim1 (typically y) along columns. Result matrices are size (Ndim2, Ndim1) — z-fastest, matching state vector ordering.

Y, Z = meshgrid(domain, :y, :z)
U = @. Z - 0.5 + 0.1 * sin(2π * Y)  # 2D field, N_z × N_y
prob[:U] = vec(U)  # flatten for DSL

print_evals(λs::Vector{<:Complex})

Pretty-print every eigenvalue in λs, highest index first, in the order given — this does NOT sort. Call sort_evals first if you want a particular ordering.

print_summary(r::SolverResults)

Print a compact summary of one solve result.

reconstruct(cache, prob, eigvec, k, field_name::Symbol) -> Vector{ComplexF64}

Reconstruct a derived variable's coefficient vector from an eigenvector.

Given the solved eigenvector eigvec (containing coefficients for all declared variables) and wavenumber k, computes the derived field using its stored operator inverse: v = H(k) * rhs(eigvec).

cache = discretize(prob)
results = solve(cache, [0.1]; sigma_0=0.02)
eigvec = results[1].eigenvectors[:, 1]

# Reconstruct derived variables:
v_coeffs = reconstruct(cache, prob, eigvec, 0.1, :v)
u_coeffs = reconstruct(cache, prob, eigvec, 0.1, :u)

reconstruct_all(cache, prob, eigvec, k) -> Dict{Symbol, Vector{ComplexF64}}

Reconstruct ALL fields (declared variables + derived variables) from an eigenvector.

fields = reconstruct_all(cache, prob, eigvec, 0.1)
# fields[:w], fields[:zeta], fields[:b] — from eigenvector
# fields[:v], fields[:u] — reconstructed from @derive definitions

Accept Symbol for remove_evals as well (convenience).

solve(cache, k_values; sigma_0, nev=1, which=:LM, tol=1e-10, maxiter=300,
      ncv=0, mat_solver=:mumps, eps_type=:krylovschur,
      n_tries=8, Δσ₀=0.2, incre=1.2, ϵ=1e-5,
      manage_init=true, verbose=false) -> Vector{SolverResults}

Distributed eigensolve via SLEPc over PETSc, one generalized pencil per wavenumber spread across all MPI ranks, with adaptive-σ shift-and-invert. Provided by the package extension BiGSTARSMPIExt, which loads only when MPI, PetscWrap, and SlepcWrap are all imported. Returns Vector{SolverResults}, fully populated on rank 0; other ranks return empty markers.

Requires a complex-scalar PETSc/SLEPc build (--with-scalar-type=complex, with PETSC_DIR/PETSC_ARCH/SLEPC_DIR exported). Run under mpiexec -n P julia --project=. script.jl. See docs/src/mpi.md.

λs_sorted, χ_sorted = sort_evals(λs, χ, "R"; sorting="lm")
Sort the eigenvalues `λs` and corresponding eigenvectors `χ` based on the specified criterion 
(`"M"` for magnitude, `"I"` for imaginary part, or `"R"` for real part).
The `sorting` argument determines the order: `"lm"` for descending order.

sort_evals(λ, Χ, by::Symbol; rev::Bool=true)

Sort eigenvalues λ and corresponding eigenvectors Χ by:

• :R → real part
• :I → imaginary part
• :M → magnitude (abs)

Set rev=true for descending (default), false for ascending.

to_coefficients(f_physical, coord_type, N) -> Vector

Transform physical-space values to spectral coefficients.

• :chebyshev → Chebyshev T coefficients via DCT
• :fourier → Fourier coefficients via FFT (scaled by 1/N)

to_physical(c, domain, coord; x=nothing) -> Vector

Domain-aware evaluation of spectral coefficients. Unlike the coord_type method (whose x must be reference coordinates ∈ [-1, 1]), this overload maps the supplied physical points to reference coordinates using domain's bounds before evaluating — so you can pass x = gridpoints(domain, coord) directly.

• Chebyshev coord: maps x ∈ [lower, upper] → [-1, 1], then evaluates the expansion.
• Fourier coord: inverse-FFTs (the x argument is ignored; the FFT grid is implied).

to_physical(c, coord_type, x) -> Vector

Transform spectral coefficients back to physical-space values.

• :chebyshev → evaluate Chebyshev expansion at points x
• :fourier → inverse FFT