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: Sort eigenvalues by (:R, :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)

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

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

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)

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

λ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