BiGSTARS.jl Documentation

Bi-global stability analysis in Julia

Sparse spectral eigenvalue problems from readable equations.

BiGSTARS.jl turns physical-space notation such as dx, dy, and dz into sparse generalized eigenvalue matrices for rotating, stratified, and geophysical flows.

Start with the DSL Open an example

Overview

Bi-global analysis bridges the gap between 1D stability tools (often too idealized) and fully 3D tri-global solvers (often too costly). BiGSTARS.jl offers a practical middle ground for rotating, stratified flows.

Equation DSL

Write governing equations directly with dx, dy, dz, substitutions, and boundary-condition macros.

Sparse spectral operators

Chebyshev directions use ultraspherical operators; Fourier directions stay diagonal in coefficient space.

Wavenumber caching

Discretize once, then assemble quickly for scalar wavenumbers or named directions such as k_x and k_y.

Solver workflow

Use Arnoldi, ARPACK, or KrylovKit with adaptive shift-and-invert and thread-parallel wavenumber sweeps.

What To Read First

Equation DSL

Define domains, variables, equations, boundary conditions, and derived variables.

Eigenvalue Solver

Choose a solver, shift, tolerance, and number of eigenvalues.

Post-processing

Reconstruct fields from eigenvectors and convert them to physical space.

Examples

Compare your workflow with complete Eady, Stone, and rotating RBC examples.

Quick Start

using BiGSTARS

# Define the domain
domain = Domain(
    x = FourierTransformed(),       # wavenumber direction
    y = Fourier(60, [0, 1]),        # periodic, N=60, domain [0, 1)
    z = Chebyshev(30, [0, 1])       # bounded, N=30, domain [0, 1]
)

# Set up eigenvalue problem
prob = EVP(domain, variables=[:psi], eigenvalue=:sigma)

# Background state and parameters
Z = gridpoints(domain, :z)
prob[:U] = Z .- 0.5
prob[:E] = 1e-12

# Define operators and equations
@substitution Lap(A) = dx(dx(A)) + dy(dy(A)) + dz(dz(A))
@equation sigma * Lap(psi) = U * dx(Lap(psi)) - E * Lap(Lap(psi))

# Boundary conditions
@bc left(psi) = 0
@bc right(psi) = 0

# Discretize once, solve over wavenumbers
cache = discretize(prob)
k_values = range(0.1, 5.0, length=50)
results = solve(cache, k_values; sigma_0=0.02, method=:Krylov)

# For problems without wavenumber direction:
results = solve(cache; sigma_0=0.02)

For domains with more than one FourierTransformed() direction, assemble with named wavenumbers:

domain = Domain(
    x = FourierTransformed(),
    y = FourierTransformed(),
    z = Chebyshev(40, [0, 1])
)

cache = discretize(prob)
A, B = assemble(cache; k_x=1.0, k_y=0.5)

Keyword names are checked. Use either the coordinate name (x=1.0) or the explicit wavenumber name (k_x=1.0), but do not provide conflicting aliases.

Examples

Eady example — QG PV baroclinic instability with dynamic BCs
Stone example — Non-hydrostatic Boussinesq with derived variables
rRBC example — Rotating Rayleigh-Benard with constraint equations