API Reference

Pre-allocated workspace for in-place assembly. One per thread.

Create with allocate_workspace(cache), reuse across wavenumbers.

ConvergenceHistory

Tracks convergence information during eigenvalue computation.

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

EigenSolver

Main solver object for generalized eigenvalue problems Ax = λBx.

Fields

• A::AbstractMatrix: Left-hand side matrix
• B::AbstractMatrix: Right-hand side matrix
• config::SolverConfig: Solver configuration
• results::Union{SolverResults, Nothing}: Latest results (if available)

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

SolverConfig

Configuration parameters for eigenvalue solvers.

Fields

• method::Symbol: Solver method (:Arnoldi, :Arpack, :Krylov)
• σ₀::Float64: Initial shift value
• which::Symbol: Which eigenvalues to find (:LM, :LR, :SR, etc.)
• nev::Int: Number of eigenvalues to compute
• maxiter::Int: Maximum iterations
• tol::Float64: Convergence tolerance
• sortby::Symbol: Order of returned eigenvalues — :nearest (default, distance from the shift σ₀, i.e. the mode shift-and-invert targets), :R (real part), :I, :M
• n_tries::Int: Number of retry attempts
• Δσ₀::Float64: Initial shift increment
• incre::Float64: Increment growth factor
• ϵ::Float64: Successive eigenvalue tolerance
• krylovdim::Int: Krylov subspace dimension (Krylov method only). Default 30; raise it only for hard/clustered spectra (cost scales steeply with this). It is clamped at runtime to [nev+2, n].

SolverResults

Contains results from eigenvalue computation.

Fields

• eigenvalues::Vector{ComplexF64}: Computed eigenvalues
• eigenvectors::Matrix{ComplexF64}: Corresponding eigenvectors
• converged::Bool: Whether computation converged
• method_used::Symbol: Method that produced the result
• final_shift::Float64: Final shift value used
• iterations::Int: Total iterations
• solve_time::Float64: Wall-clock time for solution
• history::ConvergenceHistory: Detailed convergence information

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

allocate_workspace(cache) -> AssemblyWorkspace

Allocate a dense (A, B, temp) workspace matching the cache dimensions. Dense because eigenvalue solvers need dense matrices for factorization anyway, and in-place overwrites require fixed memory layout.

assemble!(ws, cache, k)

Assemble A and B into the pre-allocated workspace ws for wavenumber k. Overwrites ws.A and ws.B in-place. Derived-variable terms still allocate while rebuilding their k-dependent inverse operators.

assemble(cache, k) -> (A, B)

Assemble the full A and B matrices for a given wavenumber k (allocating). A(k) = Σp k^p * Ap + derived variable contributions with H(k).

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_coefficients(f_values) -> Vector{Float64}

Compute Chebyshev expansion coefficients from function values at Chebyshev-Gauss-Lobatto points (in the ordering from chebyshev_points, i.e., from x=1 down to x=-1).

Uses the DCT-I (type 1) relation. Given N function values at the CGL points x_j = cos(j*pi/(N-1)), j=0,...,N-1, the Chebyshev coefficients are:

ck = (2 / (N-1)) * sum{j=0}^{N-1} '' f(x_j) * cos(kjpi/(N-1))

where '' means the first and last terms are halved.

The returned vector has length N, where c[k] is the coefficient of T_{k-1}.

chebyshev_evaluate(c, x) -> Vector{Float64}

Evaluate a Chebyshev expansion with coefficients c at points x ∈ [-1, 1]. Uses Clenshaw's algorithm for numerical stability.

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

compare_methods!(solver::EigenSolver; methods=[:Arnoldi, :Arpack, :Krylov], verbose=true)

Compare multiple methods on the same problem.

Returns

• Dict{Symbol, SolverResults}: Results for each method

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.

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.

get_results(solver::EigenSolver)

Extract eigenvalues and eigenvectors from solved system.

Returns

• (λ, Χ): Eigenvalues and eigenvectors

Throws

• ArgumentError: If solver hasn't been solved or didn't converge

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::AbstractVector{<:Number}; n::Int=length(λs), by::Symbol=:abs)

Pretty-print the top n eigenvalues from the list λs, optionally sorted by :real, :imag, or :abs.

print_summary(solver::EigenSolver)

Print a summary of solver results.

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, method=:Krylov, parallel=false, kwargs...)

Solve the eigenvalue problem for each wavenumber in k_values. Returns a Vector{SolverResults}.

Uses in-place assembly to reuse the main dense matrix workspace. When parallel=true, one pre-allocated workspace per thread is created. Problems with derived variables still allocate per wavenumber while rebuilding H(k).

Keyword arguments

• sigma_0::Float64 — initial shift for the eigenvalue solver (required)
• method::Symbol — solver method: :Krylov (default), :Arnoldi, or :Arpack
• parallel::Bool — if true, solve wavenumbers in parallel using Threads.@threads. Start Julia with multiple threads (julia -t auto) to benefit from this.
• verbose::Bool — print solver progress (default false)

Example

cache = discretize(prob)
k_values = range(0.1, 5.0, length=50)

# Sequential (1 workspace, reused for all wavenumbers):
results = solve(cache, k_values; sigma_0=0.02)

# Parallel (1 workspace per thread, reused across wavenumbers):
results = solve(cache, k_values; sigma_0=0.02, parallel=true)

solve!(solver::EigenSolver; verbose::Bool=true, A_buf=nothing)

Solve the generalized eigenvalue problem with adaptive shift selection.

Returns

• solver.results: Updated with solution results

Example

solver = EigenSolver(A, B; σ₀=1.0, method=:Arnoldi, nev=5)
solve!(solver)
λ, Χ = get_results(solver)

solve(cache; sigma_0, kwargs...) -> Vector{SolverResults}

Solve a problem with no FourierTransformed direction (pure 1D/2D, no wavenumber sweep).

solve_mpi(cache, k_values; sigma_0, nev=1, which=:LM, tol=1e-10, maxiter=300, kwargs...)

Distributed-memory eigensolve via SLEPc over PETSc, one eigenproblem per wavenumber spread across all MPI ranks. 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.

Run under mpiexec -n P julia script.jl with SlepcInitialize() / SlepcFinalize() bracketing the work, and a complex-scalar system PETSc/SLEPc.

λ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, coord_type, x) -> Vector

Transform spectral coefficients back to physical-space values.

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