Skip to content

LinStabPy

LinStabPy builds the discrete linearized compressible Navier–Stokes operator about a given base flow and solves the eigenvalue and resolvent problems on it at HPC scale.

Intuition first. You bring a steady flow field computed elsewhere. LinStabPy answers two questions. Stability: does a small disturbance grow on its own, at what frequency, and with what spatial shape? Resolvent: at a given frequency, which steady harmonic forcing produces the largest response, and how large is that amplification?

The discrete state is the primitive perturbation vector \( \mathbf{q}' \) at each grid point, giving 5 unknowns per point before boundary elimination (\( \rho', u', v', w', T' \)). The code assembles the sparse matrix \( \mathcal{L}_0 \) such that \( \partial_t \mathbf{q}' = \mathcal{L}_0 \mathbf{q}' \) in the undisturbed interior, then applies sponge, dissipation, and boundary recipes in a fixed order. Everything downstream consumes that one operator family.

Governing Equations & Mathematical Formulation

Perturbation Decomposition

A two-dimensional steady compressible base flow \( \bar{\mathbf{q}}(x,y) \) is perturbed by small three-dimensional unsteady fluctuations:

\[ \mathbf{q}(x, y, z, t) = \bar{\mathbf{q}}(x, y) + \mathbf{q}'(x, y, z, t) \]

Assuming periodicity in the spanwise direction \( z \) and harmonic time dependence with spanwise wavenumber \( \beta \) and frequency \( \omega \):

\[ \mathbf{q}'(x, y, z, t) = \hat{\mathbf{q}}(x, y) \exp\left( i (\beta z - \omega t) \right) \]

where the state vector of primitive perturbation variables is:

\[ \hat{\mathbf{q}} = \begin{bmatrix} \hat{\rho} & \hat{u} & \hat{v} & \hat{w} & \hat{T} \end{bmatrix}^T \]

Linearized Dynamical System (LNSE)

Linearizing the compressible Navier–Stokes equations around \( \bar{\mathbf{q}} \) yields the semi-discrete linear dynamical system:

\[ \frac{\partial \mathbf{q}'}{\partial t} = \mathcal{L}_0(\beta) \mathbf{q}' + \mathbf{B} \mathbf{f}' \]

where \( \mathcal{L}_0 \) incorporates summation-by-parts (SBP) finite differences (orders 2, 4, 6, 8), coordinate metrics, characteristic sponge layers, artificial dissipation, and boundary conditions (with mass/descriptor inversion applied during assembly).

Under harmonic time dependence \( \mathbf{q}'(t) = \hat{\mathbf{q}} \exp(-i\omega t) \), the unforced homogeneous system reduces to:

\[ i \mathcal{L}_0(\beta) \hat{\mathbf{q}} = \omega \hat{\mathbf{q}} \]
  • Angular frequency: \( \operatorname{Re}(\omega) \)
  • Temporal growth rate: \( \operatorname{Im}(\omega) > 0 \) denotes exponential instability (\( \propto \exp(\operatorname{Im}(\omega) t) \))
  • Solved via SLEPc EPS using Krylov–Schur with spectral shift-invert transformation.
  • For spatial stability at real \( \omega \), the system rearranges into the generalized pencil \( (L_c + \omega L_w) \hat{\mathbf{q}} = \alpha R \hat{\mathbf{q}} \) for complex streamwise wavenumber \( \alpha \).

Input–Output Resolvent Analysis

Under external harmonic forcing \( \mathbf{f}'(t) = \mathbf{B} \hat{\mathbf{f}} \exp(-i\omega t) \), the resolvent operator \( \mathcal{H}(\beta, \omega) \) maps inputs to linear responses \( \hat{\mathbf{q}} = \mathcal{H}(\beta, \omega) \hat{\mathbf{f}} \):

\[ \mathcal{H}(\beta, \omega) = \mathbf{C} \mathbf{F} \left( \mathcal{L}_0(\beta) - (\gamma + i\omega) \mathbf{I} \right)^{-1} \mathbf{F} \mathbf{B} \]

where \( \mathbf{B} \) and \( \mathbf{C} \) are input and output restriction masks, \( \mathbf{F} \) is an optional spatial filter placed inside the masks, and \( \gamma \ge 0 \) is an optional real discount.

The optimal energy gain \( \sigma_1^2 \) is defined as the supremum over all forcing shapes under the Chu compressible energy norm:

\[ \sigma_1^2(\beta, \omega) = \max_{\hat{\mathbf{f}} \neq 0} \frac{\| \mathcal{H} \hat{\mathbf{f}} \|_E^2}{\| \hat{\mathbf{f}} \|_F^2} = \max_{\hat{\mathbf{f}} \neq 0} \frac{\langle \hat{\mathbf{q}}, \mathbf{W}_q \hat{\mathbf{q}} \rangle}{\langle \hat{\mathbf{f}}, \mathbf{W}_f \hat{\mathbf{f}} \rangle} \]

Computed as \( \sigma = \sqrt{\operatorname{Re}(\lambda)} \) from the weighted normal eigenproblem on a PETSc MatShell. Large gains when all modes are stable (\( \operatorname{Im}(\omega) < 0 \)) quantify non-normal convective amplification, receptivity to environmental disturbances, and flow-control sensitivity.

Reading guide

This wiki keeps solver language exact and explains the fluids side alongside:

  • Base flow: A steady solution of the compressible Navier–Stokes equations, loaded from .npz files. We linearize about it; the code never computes it.
  • Temporal stability: Eigenvalues \( \omega = \operatorname{eig}(i \mathcal{L}_0) \). Growth in time is \( \operatorname{Im}(\omega) > 0 \). The shift-invert target you pass in YAML lives in \( \omega \) space.
  • Spatial stability: Same physics, swapped role of axes. Frequency is prescribed real, streamwise wavenumber is the eigenvalue.
  • Resolvent: For fixed \( (\beta, \omega) \), the gain is the leading weighted singular value of \( \mathcal{R} \). High gain with all \( \operatorname{Im}(\omega) < 0 \) indicates transient or convective amplification.
  • Chu norm: The compressible energy inner product. Diagonal in our discretization, built from quadrature weights, base density and temperature, \( \gamma \), and Mach number.
  • Solver objects: KSP, PC, EPS, ST, and MatShell are defined in the solver glossary, with links to the official PETSc and SLEPc docs.

  • Spectra & Sensitivity

    Growth rates, frequencies, adjoint modes, and where the flow is most sensitive to a structural perturbation.

    Read Overview

  • LNSE Discretization

    Grids, SBP difference stencils (orders 2, 4, 6, 8), metric invariants, and energy-stable boundary condition recipes.

    Explore Architecture

  • Orchestrated Workflow

    End-to-end staged pipeline from mean-flow ingest and YAML config to MPI execution and contract-stable NPZ payloads.

    Study Workflow

  • Resolvent & Solvers

    Cost-aware parametric sweeps across \( (\beta, \omega) \) with anchor factorizations, preconditioner reuse, and residual gates.

    Module Guide

Quick Example

You can run a case from the terminal with MPI, or import the package in Python. Import stays lightweight and does not start MPI or PETSc until a solver stage is invoked.

# Run one case on 4 MPI ranks. Config and output travel by env vars (see Getting started).
export LINSTABPY_CASE_CONFIG=case.yaml LINSTABPY_CASE_OUTPUT_NPZ=results/out.npz
mpirun -n 4 linstabpy

# Load and validate a case configuration without pulling PETSc
import linstabpy
from linstabpy.config.typed import load_case_config

config = load_case_config("case.yaml")
print(f"Case: {config.meta.case_name}, Mode: {config.solver.mode}")

The YAML file is the case definition. It names the base flow files, the Reynolds and Mach numbers, the grid and boundary choices, and what to compute. Keep it next to the results for reproducible scientific provenance.

Results land in contract-stable .npz files. Load them with numpy.load():

import numpy as np

data = np.load("results/case_temporal.npz")
omega = data["omega"]          # Complex eigenvalues
modes = data["modes"]          # Eigenvectors [n_points, 5, n_modes]
print(f"Computed {len(omega)} eigenvalues. Leading growth rate: {np.max(np.imag(omega)):.4e}")

Documentation Directory

Section Core Focus Target Audience
Getting started Install, first runs, examples, your first case Start here
Overview What LinStabPy computes and where it stops First read, scope, vocabulary
Architecture Component relationships & C4 models Contributors, system developers
Workflow Steps of a run in order Running and orchestrating cases
Module guide Per-package derivations and code notes Advanced customization, algorithms
Boundaries YAML schema, CLI flags, env vars, NPZ contract Reference when building cases
Solver glossary Every KSP, PC, EPS, ST, and MatShell term defined and linked Open beside any solver page
Core concepts Base flow, growth rate, gain, Chu norm, wavemaker, sponge, and friends Open beside any science page

Cite linstabpy

If this software contributed to published work, please cite the linstabpy paper together with the LinStab2D methodology paper:

Hamada, G. Y. R. and Frantz, R. A. S. linstabpy: open-source compressible stability and resolvent analysis. Theor. Comput. Fluid Dyn. (in preparation, 2026).

Martini, E. and Schmidt, O. T. (2024). LinStab2D: stability and resolvent analysis of compressible viscous flows in MATLAB. Theor. Comput. Fluid Dyn. 38, 665–685. doi:10.1007/s00162-024-00706-0.

The full bibliography lives on the references page, and CITATION.cff carries the machine-readable entry (DOI and version land with the first tagged Zenodo release).