linstabpy 1.0.0

Documentation: https://linstabpy.vercel.app/

linstabpy performs sparse linear stability and resolvent analysis of compressible viscous flows. It builds the linearized compressible Navier–Stokes operators for two-dimensional structured base flows and solves the eigenvalue and resolvent problems through PETSc/SLEPc with sparse direct factorization (MUMPS).

It began as a Python port of the MATLAB LinStab2D code (Martini and Schmidt, 2024) and has grown into a standalone framework: it adopts the operator abstractions of resolvent4py (Padovan et al., 2025), extends the discretization to fd2/fd4/fd6/fd8 summation-by-parts operators, and verifies the assembled operator against independent derivations rather than against MATLAB output alone.

Memory was the original motivation. Stability and resolvent analysis of 2D base flows produces large linear systems, and MATLAB's dense backslash path scales as O(n^2) in storage. A 150k-DOF operator needs about 180 GB of RAM that way, which puts most LES/DNS-derived base flows out of reach on an ordinary workstation. Keeping the operators sparse and sending the linear algebra through PETSc/SLEPc and MUMPS brings those same problems back inside a workstation's budget.

Contents

• Overview
• Installation
• Quick start
• Repository layout
• Numerical setup
• Nondimensionalization
• Boundary conditions
• Parity with LinStab2D
• Case studies
• Development
• References
• License and citation

Overview

The package targets two-dimensional structured base flows and builds operators for the linearized compressible Navier–Stokes equations. The public project name is LinStabPy; the repository, distribution, and import package name is linstabpy.

It currently supports:

• Resolvent analysis: leading gain, forcing modes, and response modes
• Temporal stability: direct and adjoint eigenvalue problems
• Spatial stability: streamwise-wavenumber eigenvalue problems
• Coordinate systems: Cartesian, cylindrical/axisymmetric, and curvilinear structured meshes
• Boundary treatment: Dirichlet rows, Neumann rows, sponge damping fields, and masking for structured grids with excluded points
• Energy norms: Chu compressible energy norm with quadrature weighting

Conventions used in this repo. Temporal eigenvalues are omega = eig(i * L0) with exp(-i omega t) time dependence, so Im(omega) > 0 is unstable. The resolvent is H = C F ((gamma - i omega) I - L0)^{-1} F B with gains from the weighted normal problem (the sign is explained in wiki/resolvent-sign-convention.md; solver.resolvent_sign_convention: matlab_parity restores the LinStab2D sign). SBP orders are 2, 4, 6, 8 only. The default shifted solve is KSPPREONLY with PCLU through MUMPS. If you come from math and HPC, the new vocabulary is mostly fluids side: base flow, sponge, Chu norm, wavemaker. The solvers are standard shift-invert Krylov-Schur.

The compressible formulation follows Martini and Schmidt (2024), and the operator-handling design draws on resolvent4py (Padovan et al., 2025). Where the original code was validated mainly against its own MATLAB output, LinStabPy adds a layer of provenance-independent checks: a symbolic re-derivation of the coefficients, a complex-step Jacobian, and direct verification of the curvilinear transform, global assembly, boundary rows, and energy norm.

Installation

One supported numerical stack: CPU, complex-scalar PETSc/SLEPc, MUMPS. On Linux:

git clone https://github.com/openfluids/linstabpy.git
cd linstabpy

bash install.sh --diagnose   # read-only checks, installs nothing
bash install.sh              # builds .venv with PETSc/SLEPc, MUMPS and this package
bash install.sh --check      # confirms the solver stack works

The full guide covers conda, macOS, clusters, Intel MPI, GPU stacks and the MPI rules: Installation. When a build fails, see Troubleshooting.

Reference base flows

The Martini example loaders (example_1_blasius, example_2_jet, example_4_cavity, example_5_paraboloid) read four .mat base flows from the upstream LinStab2D project (GPL-3.0): M09jet.mat, BL_BF.mat, Step_BF.mat, Paraboloid_BF.mat. They ship in examples/linstab2d/baseflows/, so the examples run from a fresh clone. To read them from a LinStab2D clone instead, set:

git clone https://github.com/LinStab2D/LinStab2D.git ~/linstab2d-reference
export LINSTABPY_BASEFLOW_DATA_DIR=$HOME/linstab2d-reference

The Song example ships its 32° base flow, reference eigenvalues and leading mode under examples/song/reference/ (redistributed with the authors' permission); Song's raw finite-volume files are not shipped. The Savarino example ships no data from its Zenodo record: examples/savarino/fetch_zenodo.py downloads the base flow and modes when needed.

Quick start

git clone https://github.com/openfluids/linstabpy.git
cd linstabpy

bash install.sh --diagnose
bash install.sh
bash install.sh --check
bash benchmark.sh

This repository has no required submodules. The data each example needs is listed in Reference base flows.

Run the smallest Blasius example:

.venv/bin/python -m examples.linstab2d.example_1_blasius

Run a full example:

.venv/bin/python -m examples.linstab2d.example_1_blasius

Do not use uv run on macOS for PETSc/SLEPc examples. Use the interpreter installed and patched by install.sh:

.venv/bin/python -m examples.linstab2d.example_1_blasius

Collect the test suite:

.venv/bin/python -m pytest --collect-only tests

Run the suite with:

.venv/bin/python -m pytest tests/

MATLAB-parity tests are opt-in

This repository does not ship MATLAB reference data — neither the raw .mat exports nor the derived .npz bundles. The tests that compare against that data are therefore deselected by default via the matlab_parity marker, and run only when you ask for them:

LINSTABPY_MATLAB_PARITY=1 .venv/bin/python -m pytest -m matlab_parity

Two consequences are deliberate:

• A default run reports no skips it cannot justify. These tests are deselected, not skipped, so they never pad the skip count. Every skip in a default run is an absent optional backend (STRUMPACK, CUDA, superlu_dist) — a true statement about the machine you are on.
• Requesting them means they must run. With LINSTABPY_MATLAB_PARITY=1, a missing reference is a failure, not a skip. A requested comparison that quietly does not execute is indistinguishable from one that passed.

Supply the reference data through LINSTABPY_BASEFLOW_DATA_DIR, the martini_examples/ export directory, LINSTABPY_TEST_REFERENCE_SPECTRUM, and the examples/linstab2d/ bundles. The failure message names the exact paths it wanted.

Canonical example runner

Use run_all_examples.sh as the single example runner for serial acceptance, cache cleanup, parity checks, and MPI launches. The runner computes each selected Martini example, then runs the matching example_*_plots module. Direct .venv/bin/python -m examples.linstab2d.example_1_blasius writes the NPZ only. For a figure-by-figure walkthrough of the Martini examples, see docs/tutorial/tutorial.pdf, Section 5. That section also explains how to read the generated baseflow, spectrum, and mosaic figures, and records the plotting improvements proposed for the next figure pass.

Task	Command
List available examples	bash run_all_examples.sh --list
Run all examples, plots, and Martini parity	bash run_all_examples.sh
Re-plot one example from an existing NPZ	.venv/bin/python -m examples.linstab2d.example_1_plots
Recompute selected examples from cached PNGs/NPZs	bash run_all_examples.sh --fresh --example 1 --skip-parity
Remove selected generated NPZ/PNG outputs first	bash run_all_examples.sh --clean --example 1 --skip-parity
Run one example under MPI	bash run_all_examples.sh --mpi -n 10 -m examples.linstab2d.example_3_streaky_bl
Pass PETSc/SLEPc options to one example	bash run_all_examples.sh --mpi -n 1 -m examples.linstab2d.example_1_blasius -- -eps_monitor

--mpi folds the old one-module MPI wrapper behavior into the same script: the runner probes a working MPI launcher, defaults to half of the detected physical cores when -n/--ranks is omitted, and caps OpenMP/BLAS thread counts to one thread per rank unless the caller already set them. Override the runtime with LINSTABPY_MPIEXEC, LINSTABPY_MPI_RANKS, or LINSTABPY_PYTHON when a machine needs explicit paths.

Why YAML configs?

Each example is a numerical case definition: Reynolds number, Mach number, grid size, boundary conditions, sponge settings, forcing masks, solver targets, and PETSc/MPI choices should stay visible and diffable, with an explanation sitting next to each value.

Plain JSON would remove those explanations. Long command-line flags would hide too much state in shell history or runner scripts. Python config files would mix data with executable code. YAML keeps the comments next to the knobs.

Repository layout

For a maintainer map of where numerical policy belongs, see docs/architecture.md.

.
|-- src/linstabpy/                 Importable Python package
|   |-- baseflow/                  Packed-field loader, BaseFlow, viscosity
|   |-- boundary_conditions/       Dirichlet/Neumann row replacement and recipes
|   |-- config/                    YAML schema and typed solver policy
|   |-- mesh/                      SBP finite differences, grids, filters, deformation
|   |-- operators/                 LNSE coefficients, assembly, spatial operators, Chu norm
|   |-- solvers/                   PETSc/SLEPc resolvent and eigenvalue solvers
|   |-- driver/                    Shared run path, payload, and NPZ write
|   `-- utils/                     YAML case loader, PETSc helpers
|-- examples/                      Public cases, figure scripts, and tracked PNGs
|   |-- linstab2d/                 Six Martini examples and MATLAB field constructors
|   |-- odd_even/                  Odd-even decoupling on a wavy-wall mesh (paper figure)
|   |-- song/                      Song / Robinet SWBLI EVP case
|   |-- savarino/                  Savarino 30.8 deg resolvent case
|   `-- sawtooth/                  Lui 2022 suction-side SWBLI, stride-2 field
|-- docs/                          Architecture, backends, and tutorial
|-- tests/                         Pytest suite, file-size and provenance checks
|-- install.sh                     Installer and host preflight (`--diagnose`)
|-- environment.yml                CPU conda-forge complex PETSc + pip from pyproject
|-- benchmark.sh                   Install benchmark
`-- run_all_examples.sh            Serial, parity, and MPI example runner

Keep reusable code under src/linstabpy/ and the public cases under examples/linstab2d, examples/song, examples/savarino, and examples/sawtooth. Martini .mat constructors live in examples/linstab2d/matlab_baseflow.py. Set LINSTABPY_BASEFLOW_DATA_DIR if those files sit outside the example directory.

Data and file-size policy

Tracked GitHub-facing files must stay at or below 60 MiB. Large generated outputs, raw MATLAB exports, and machine-local benchmark products are kept out of Git and either regenerated or stored externally. The guardrail is:

.venv/bin/python -m pytest tests/test_tracked_file_size.py -q --tb=line

The same guard runs in the GitHub guardrails workflow. Ignored local outputs above the report threshold are listed separately so they can be cleaned or archived without blocking normal source commits.

Compact Martini fixtures belong under examples/linstab2d/. Raw MATLAB exports such as operator_convention.mat and coefficients_raw.mat are intentionally not tracked when they exceed the 60 MiB policy.

Numerical setup

Temporal stability

The current temporal stability path uses:

• SLEPc EPS
• KRYLOVSCHUR as the outer eigensolver
• shift-invert spectral transformation, ST=SINVERT
• direct LU factorization with MUMPS for shifted linear systems in the default production path

The temporal eigenproblem is formed on i * L0, so reported temporal eigenvalues are omega = eig(i * L0).

Residual bookkeeping

The code records several residual-like quantities:

• SLEPc EPS residuals from EPS.computeError(...)
• the configured EPS convergence test: abs, rel, or norm
• the selected reported residual type: absolute, relative, or backward
• an explicit postcheck residual computed outside SLEPc on returned modes
• boundary-condition leakage diagnostics

These quantities are related but not interchangeable. EPSSetTolerances() sets the tolerance and maximum iteration count used by the EPS convergence test. EPSSetConvergenceTest() selects how the EPS convergence test is measured. EPSComputeError() reports an error estimate for a returned eigenpair, but that estimate is not always the same quantity used internally during the solve.

For shift-invert runs, true_residual=true is useful because cheap transformed residual estimates can be optimistic for the original problem.

PETSc KSP tolerances

KSPSetTolerances() controls a linear Krylov solve: relative tolerance, absolute tolerance, divergence tolerance, and maximum iteration count. It is not an eigenpair tolerance.

In the default direct path, shifted systems use:

KSPPREONLY + PCLU + MUMPS

In that configuration, iterative KSP tolerances are not the main accuracy knobs. The dominant knobs are EPS settings, Krylov subspace size ncv, target placement, and the conditioning of the shifted operator.

When an iterative linear-solve path is selected, ksp_rtol, ksp_atol, ksp_max_it, pc_type, subdomain settings, and ksp_error_if_not_converged become meaningful.

Compute backends

The validated sparse-direct default is PETSc LU with MUMPS. Backend selection is centralized in the factor-backend seam and applies only to PETSc LU paths, not to iterative Krylov-only solves.

Selection precedence is:

1. YAML solver subsection factor_backend
2. per-solver environment override, such as LINSTABPY_EIGEN_PC_FACTOR_SOLVER_TYPE or LINSTABPY_RESOLVENT_PC_FACTOR_SOLVER_TYPE
3. LINSTABPY_FACTOR_BACKEND
4. the call-site default

Examples:

LINSTABPY_FACTOR_BACKEND=mumps .venv/bin/python -m examples.linstab2d.example_1_blasius
LINSTABPY_FACTOR_BACKEND=strumpack .venv/bin/python -m examples.linstab2d.example_1_blasius

Unavailable backends fail during PETSc LU setup with a direct error message instead of silently falling back to MUMPS or PETSc native LU.

Backend	Status	Notes
mumps	Validated default	CPU sparse-direct backend used by the current test suite.
strumpack	Experimental, run it on one host thread	Not built by default; GPU claims are unvalidated for LinStabPy matrices. Multithreaded it returns a wrong solve now and then: see the note below.
superlu_dist	Experimental	Not built by default; PETSc requires a SuperLU_DIST build with GPU offload support for GPU use.
cudss	Validated GPU backend	Not a PETSc factor type: the cuDSS LU (through nvmath-python, pip install "linstabpy[cudss]") preconditions the right-preconditioned FGMRES path. One rank, one GPU; the adjoint uses a second factorization of the conjugate transpose. Solutions match MUMPS to 1e-12; see docs/COMPUTE_BACKENDS.md.

See docs/COMPUTE_BACKENDS.md for the backend audit, version floors, and the current GPU go/no-go notes.

STRUMPACK is not thread safe

LinStabPy sets the OpenMP thread count to one before it uses STRUMPACK, because the library returns a wrong solve now and then on several threads, without any error. The line in the log is STRUMPACK: host OpenMP threads set to 1. MUMPS, cuDSS and the BLAS keep their threads, which use a different OpenMP runtime. On a system with no libgomp the package cannot set the count; it warns, and there you must set OMP_NUM_THREADS=1 yourself.

What we measured (September 2026, one workstation, one NVIDIA A6000): factor K = L0 - sigma*I for the Blasius temporal example, 750 unknowns, with the shift close to the spectrum, then solve it for a fixed right-hand side and record the relative residual ||Kx-b||/||b||. Sixteen separate processes per setting, identical input every time. For scale, LAPACK's dense LU on the same system gives 9.4e-13, and STRUMPACK on one thread gives 2.102e-10; a sparse solver with its own ordering is not expected to match a dense one, and both are far below the failures below.

	PETSc 3.24.2, STRUMPACK 8.0.0	PETSc 3.25.5, STRUMPACK master
one thread, GPU on	2.102e-10 sixteen times	2.102e-10 sixteen times
ten threads, GPU on	8e-11 to 2.5e-3, 2 of 16 bad	8e-11 to 7.0e-7, 1 of 16 bad
ten threads, GPU off	1e-11 to 1.1e-5, 1 of 16 bad	9e-12 to 1.0e+01, 2 of 16 bad

On one thread the answer is the same to the last bit, sixteen times, on both builds. On ten threads it scatters, with the card and without it. The newer STRUMPACK is master of May 2026, twenty-four commits after the 8.0.0 tag, and it behaves the same way. STRUMPACK's own TODO file states that the library "is not completely thread safe at the moment".

This is a property of the library, not of LinStabPy: the residual is measured against the same matrix the solver was given, so a wrong matrix would cancel out. We have not found the faulty code inside STRUMPACK.

The error reaches the answer, not only the probe. Solving the same example for its five leading eigenvalues through FGMRES with a STRUMPACK preconditioner, eight fresh processes per setting: on ten host threads the largest error against the MUMPS reference ran from 1.1e-11 to 6.9e-9, with two runs of eight above 1e-9; on one thread it was 2.209e-11 eight times out of eight. MUMPS on the same ten threads stayed at 1e-11. The flexible outer iteration does not repair it, because the error enters through the shift-invert solves that the eigensolver reads.

Two checks still stand behind the thread setting. The backend accuracy probe (eigen.backend_probe) solves one system after the factorization and stops the run when the relative residual exceeds its gate, 1e-8 by default; that catches the large failures above but not a draw that lands between 1e-10 and the gate. The eigenvalue and resolvent paths then refine with FGMRES and check the residual of the answer itself. The mumps and cudss backends do not show this behaviour.

One host thread costs wall time. On the small example above a whole run took 4.10 s against 3.64 s, three runs each. The cost on a production-size operator is not measured.

Resolvent region masks

resolvent_map accepts optional input_region and output_region dictionaries that build sharp diagonal forcing/response masks from physical coordinates and state variables. Omitted keys mean the full extent and all variables. Example configuration fragment, assuming the usual example variables base_flow and case_cfg are already loaded:

from linstabpy import resolvent_map

cfg = {
    **case_cfg,
    "input_region": {"x": [0.02, 0.06], "y": [0.0, 0.01], "vars": ["u", "v", "w", "T"]},
    "output_region": {"x": [0.04, 0.08], "vars": ["rho", "u", "v", "w", "T"]},
}
result = resolvent_map(base_flow, beta_range=[20.0], f_range=[0.15], cfg=cfg)

A windowed/component-restricted gain changes the forcing/response question; it is not an approximation of the global unmasked gain.

Boundary-condition acceptance

Returned stability modes are postchecked against:

• full-operator residual
• boundary-condition leakage
• constrained-Dirichlet leakage
• target distance
• optional spectral boxes on Re(omega) and Im(omega)

An interior eigenvalue solve can report converged Ritz pairs and still hand back modes that are unacceptable under the full operator and boundary treatment. The postcheck is what catches that.

Convergence principles

eps_tol alone does not define mode quality. For interior non-Hermitian problems, tightening eps_tol does not guarantee proportional improvement in accepted modes. Mode quality can be limited by ncv, target placement, operator conditioning, non-normality, and spectral crowding.

ncv is a quality-control parameter. If the target neighborhood is crowded or non-normal, too small an ncv can make accepted-mode quality plateau even when the solve is active.

Eigenvalues can stabilize before eigenvectors. A spectrum that appears stable is not enough; explicit postcheck residuals and boundary leakage still matter.

Shift-invert quality depends on the shifted operator. If iL0 - sigma I is poorly conditioned, residuals can plateau even with direct LU.

Target choice matters. A target on the real axis can sit in a crowded neighborhood. Moving the target slightly into the complex plane can isolate a branch more cleanly. Treat this as a numerical search choice, not a change in physics.

Branch-localize before scaling nev. Start with a small solve, test a few targets around the branch of interest, judge success from postcheck residual and target distance, then increase nev and ncv.

What to check on every run

Before trusting a spectrum, check:

• requested nev
• actual nconv
• ncv
• EPS tolerance and EPS convergence test
• whether true residuals were enabled
• whether shifted solves were direct or iterative
• explicit postcheck residual counts
• boundary-condition leakage diagnostics
• target distance for accepted and rejected modes
• wall time
• peak RSS
• MPI rank count and thread policy
• the exact YAML config and git commit that produced the artifact

What not to assume

Do not assume:

• nconv > 0 means the run succeeded
• passing the EPS convergence test implies a small full-operator residual
• lowering eps_tol is always the cheapest way to improve results
• a larger nev gives more trustworthy leading modes
• the fastest setup on one machine is the fastest setup on another
• a mean-flow linearization is physically meaningful just because the operator was assembled correctly

Nondimensionalization

All flow fields stored in BaseFlow are nondimensional. The governing equations are the compressible Navier–Stokes equations in the primitive variables [rho, u, v, w, T], linearized around the base flow. The nondimensionalization follows the standard convective scaling described below.

Reference scales

Symbol	Definition	Units	Source in code
rho_ref	freestream density	kg/m³	implicit: RHO = 1 at freestream (set in linstabpy.baseflow.grid_fields)
U_ref	freestream velocity	m/s	implicit: U = 1 at freestream
L_ref	problem reference length	m	user-supplied (e.g. displacement thickness delta* for Blasius, shock-impingement length X_sh for SWBLI)
T_ref	freestream static temperature	K	BaseFlow.T_0 (linstabpy.baseflow.baseflow)
p_ref	rho_ref * U_ref**2	Pa	confirmed by Chu norm weights in energy_norm.py
mu_ref	mu(T_0) via Sutherland	Pa·s	sutherland_air in linstabpy.baseflow.viscosity, variable mu_0
t_ref	L_ref / U_ref	s	derived

L_ref is not stored inside BaseFlow; it is the length scale used when computing Re and when building the mesh coordinate arrays passed to the loaders. Its physical interpretation depends on the case.

Nondimensional variables

Variable	Stored as	Nondim form	Freestream value
density	BaseFlow.RHO	rho* = rho / rho_ref	1
streamwise velocity	BaseFlow.U	u* = u / U_ref	1
wall-normal velocity	BaseFlow.V	v* = v / U_ref	0
spanwise velocity	BaseFlow.W	w* = w / U_ref	0
temperature	BaseFlow.T	T* = T / T_ref	1
viscosity	BaseFlow.MU	mu* = mu / mu(T_0)	1
streamwise coordinate	mesh X	x* = x / L_ref	—
wall-normal coordinate	mesh Y	y* = y / L_ref	—
time	operator eigenvalue omega	t* = t * U_ref / L_ref	—

Pressure is not stored explicitly. From the ideal gas law nondimensionalized with p_ref = rho_ref * U_ref**2:

p* = RHO * T / (kappa * Ma**2)

This Ma**2 factor is absent in acoustic-scaled codes; see LES vs linstabpy convention notes.

Dimensionless numbers

Reynolds number

Re = rho_ref * U_ref * L_ref / mu(T_0)

mu(T_0) is Sutherland's law evaluated at the freestream temperature T_0. Enters viscous terms as 1/Re; for example, coefficients.py has c_Arr_44 = -MU / Re (streamwise momentum wall-normal diffusion).

Mach number

Ma = U_ref / a_ref,  a_ref = sqrt(kappa * R * T_ref)

Enters the operator via c2 = kappa * Ma**2 (pressure-velocity coupling) and cv = 1 / (kappa * (kappa - 1) * Ma**2) (energy equation). Both are derived in sutherland_air (viscosity.py) and used throughout coefficients.py.

Prandtl number

Pr = mu(T_0) * cp / k_ref,  cp = kappa * R / (kappa - 1)

Enters thermal diffusion as 1 / (Re * Pr); for example, c_Arr_55 = -kappa * MU / Re / Pr (temperature Laplacian in energy equation).

Heat capacity ratio

kappa = cp / cv  (gamma)

Stored as BaseFlow.kappa. The code uses kappa rather than gamma; the BaseFlow docstring notes "Heat capacity ratio (gamma)" to keep the physical meaning explicit. Enters the ideal gas relation, the Chu norm weights, and all compressibility terms.

Re, Ma, Pr, and kappa are user-supplied in the YAML config under physics. The code uses them as given without internal consistency checks.

Sutherland viscosity

Dynamic viscosity follows Sutherland's law for air. The normalization reference is mu_0 = mu(T_0), the Sutherland value at the actual freestream temperature. The relevant section of viscosity.py:

    S1: float = 110.4  # Sutherland temperature [K]
    T_ref: float = 280.0  # reference temperature [K]
    mu_ref: float = 1.735e-5  # reference dynamic viscosity [Pa.s]

    Tstar = bf.T * T_0  # dimensional temperature

    # Reference viscosity at T_0
    mu_0 = mu_ref * (T_ref + S1) / (T_0 + S1) * (T_0 / T_ref) ** 1.5

    # Viscosity and its first two derivatives w.r.t. dimensional T
    mu = mu_ref * (T_ref + S1) / (Tstar + S1) * (Tstar / T_ref) ** 1.5

T_ref = 280.0 K and mu_ref = 1.735e-5 Pa.s are anchoring constants in the Sutherland formula, not the nondimensionalization scale. The operative nondim scale is mu_0, which depends on T_0. MU = mu / mu_0, so MU = 1 at freestream regardless of what T_0 is.

The temperature derivatives are nondimensionalized consistently:

dmudT*  = d(mu/mu_0) / d(T*)  =  (dmu/dT_dim) * T_0 / mu_0
d2mudT2* = (d2mu/dT_dim**2) * T_0**2 / mu_0

LES vs linstabpy convention notes

Three specific points where LES codes commonly differ.

Pressure reference. Some LES codes use p_ref = rho_ref * R * T_ref (acoustic scaling). linstabpy uses p_ref = rho_ref * U_ref**2 (dynamic pressure). The two coincide only at Ma = 1. For Ma != 1, the nondim ideal gas law differs by a factor of kappa * Ma**2:

Quantity	linstabpy	acoustic-scaled LES
p*	p / (rho_ref U_ref**2)	p / (rho_ref R T_ref)
ideal gas	p* = RHO * T / (kappa * Ma**2)	p* = RHO * T
conversion	multiply acoustic p* by kappa * Ma**2	divide linstabpy p* by kappa * Ma**2

Energy equation scaling. In linstabpy cv* = 1 / (kappa * (kappa-1) * Ma**2). An LES code that writes cv* = 1 / (kappa - 1) is using acoustic-based nondim. To convert: cv_LES = cv_linstabpy * kappa * Ma**2.

Temperature. linstabpy stores T* = T / T_ref with T* = 1 in the freestream. Some LES post-processing outputs temperature as excess (T - T_inf) / T_inf or as a ratio (T - T_w) / (T_0 - T_w). If an external base-flow file provides temperature already divided by T_ref, load it directly into BaseFlow.T. If it provides temperature in Kelvin, divide by T_0 before assigning.

Boundary conditions

Boundary setup is YAML-driven. Row-replacement boundary conditions are listed under bc.recipe; sponge damping is configured in YAML as an operator damping field or loaded from reference case data.

Example:

physics:
  sponge:
    enabled: true
    preset: tanh
    amplitude: 0.1
bc:
  penalty: 1.0e8
  recipe:
    - type: dirichlet
      borders: lt
      variables: ruvwT
    - type: dirichlet
      borders: b
      variables: ruvw
    - type: neumann
      borders: b
      variables: T
      exclude_corners: true
      curvilinear_normals: true
    - type: neumann
      borders: r
      variables: ruvwT
      exclude_corners: true
      curvilinear_normals: true
analysis:
  stability_acceptance:
    postcheck_rel_max: 5.0e-5
    target_distance_max: 0.15
    omega_real_abs_max: 1.5
    omega_imag_abs_max: 2.0
    conjugate_aware: true

Recipe fields:

Field	Meaning
type	dirichlet or neumann for row-replacement BCs. Sponge damping is configured separately under physics.sponge or loaded from case input data.
borders	Boundary letters. Common values are l, r, t, b, or combinations such as lt.
variables	Variable block selector, such as u, T, ruvw, or ruvwT.
penalty	Sits beside recipe under bc:, not inside a rule. Dirichlet: rows and columns are zeroed, then the diagonal is set to this value, pinning the DOF to zero (constrained eigenvalues land near i * penalty, far from physical targets). Neumann: the stencil row is scaled by this value. Driver default 1e8; assembly, router, and Martini examples default 1e4. Keep the default unless constrained-leakage diagnostics complain, then raise toward 1e6 before reaching higher.
exclude_corners	For Neumann rules, remove corner points that would overlap with other boundary rules.
curvilinear_normals	For Neumann rules, use per-face curvilinear normals instead of axis-aligned normals.

Use curvilinear_normals: true for deformed or body-fitted curvilinear meshes. The flag makes Neumann rows use per-face normals from the stored mesh Jacobian entries. Set it to false when the base-flow mesh is Cartesian and the curvilinear normal path is not needed; that preserves the axis-aligned path and can save a small amount of work. Older parity configurations may spell this flag explicitly to preserve historical behavior. New configs should spell the intended value instead of relying on omission.

analysis.stability_acceptance.conjugate_aware defaults to true, except when beta != 0 and the base flow has a nonzero W: the driver then sets it to false. With omega = eig(i*L0) and a real L0, the spectrum is closed under omega -> -conj(omega), so the distance to the target is

min(|omega - target|, |omega + conj(target)|)

This is the right default for real-valued operators, where oscillatory modes come in pairs omega and -conj(omega). Set conjugate_aware: false when you need strict non-redundant mode selection around one target branch, or when the operator is not real-valued and conjugate symmetry should not be admitted.

Parity with LinStab2D

Example-level parity

Reference MATLAB examples:

• examples/linstab2d/linstab2D_scripts/example_0_TestCurvilinearDerivative.m
• examples/linstab2d/linstab2D_scripts/example_1_Blasius.m
• examples/linstab2d/linstab2D_scripts/example_2_Jet.m
• examples/linstab2d/linstab2D_scripts/example_3_Streaky_BL.m
• examples/linstab2d/linstab2D_scripts/example_4_Cavity.m
• examples/linstab2d/linstab2D_scripts/example_5_Paraboloid.m

Python counterparts:

• examples/linstab2d/example_0_curvilinear_validation.py
• examples/linstab2d/example_1_blasius.py
• examples/linstab2d/example_2_jet.py
• examples/linstab2d/example_3_streaky_bl.py
• examples/linstab2d/example_4_cavity.py
• examples/linstab2d/example_5_paraboloid.py

MATLAB example	Python example	Status	Required artifacts
example_0_TestCurvilinearDerivative.m	example_0_curvilinear_validation.py	implemented	example_0_curvilinear_validation.npz (local), example_0_curvilinear_validation.png, example_0_curvilinear_mesh.png, example_0_field_df_triptych.png
example_1_Blasius.m	example_1_blasius.py	implemented	example_1_blasius.npz (local), example_1_blasius_baseflow.png, example_1_blasius_modes_forcing_response_{u,v}[_imag]_mosaic.png
example_2_Jet.m	example_2_jet.py	implemented	example_2_jet.npz (local), example_2_jet_baseflow.png, example_2_jet_modes_forcing_response_{u,v}[_imag]_mosaic.png
example_3_Streaky_BL.m	example_3_streaky_bl.py	implemented	example_3_*.npz (local), example_3_streaky_bl_baseflow.png, example_3_streaky_bl_spectrum_*.png, example_3_streaky_bl_modes_*_mosaic.png
example_4_Cavity.m	example_4_cavity.py	implemented	example_4_cavity.npz (local), example_4_cavity_baseflow.png, example_4_cavity_modes_{direct,adjoint}_*_mosaic.png
example_5_Paraboloid.m	example_5_paraboloid.py	implemented	example_5_paraboloid.npz, example_5_paraboloid_baseflow.png, modes_forcing_response_{u,v,w}[_imag]_mosaic.png (see example_5_D2y_asymm_fix.md)

Operator parity matrix

Status semantics:

• pending: not implemented yet
• implemented: code exists but parity artifact/signoff is still pending
• validated: implemented and parity-validated by artifacts
• deferred: intentionally postponed
• out_of_scope: intentionally not targeted
• blocked: cannot proceed due to an external blocker

Every required: true entry must reach status = validated.

ID	MATLAB	Python target	Status
construct_get_lin_problem	ConstructOperators/GetLinProblem.m	linstabpy/operators/assembly.py::build_lnse_operator	validated
construct_get_spatial_lin_problem	ConstructOperators/GetSpatialLinProblem.m	linstabpy/operators/spatial.py::build_spatial_operator	validated
construct_cartesian_coefficients	ConstructOperators/getCoeffsLaminar_Cartesian.m	linstabpy/operators/coefficients.py::compute_cartesian_coefficients	validated
construct_axisymmetric_coefficients	ConstructOperators/getCoeffsLaminar.m	linstabpy/operators/coefficients.py::compute_cylindrical_coefficients	validated
construct_create_diff_matrices_axy	ConstructOperators/CreateDiffMatrices_Axy.m	linstabpy/operators/assembly.py::_create_axy_diff_matrices	validated
construct_lnseoperators_cyl	ConstructOperators/LNSEoperators_cyl.m	linstabpy/operators/assembly.py::build_lnse_operator(model="axy")	validated
construct_get_comp_energy_norm	ConstructOperators/GetCompEnergyNorm.m	linstabpy/operators/energy_norm.py::compute_chu_norm	validated
bc_dirichlet	BoundaryConditions/BC_Dirichlet.m	linstabpy/boundary_conditions/dirichlet.py::apply_dirichlet	validated
bc_dirichlet_spatial	BoundaryConditions/BC_Dirichlet_SpatialStability.m	linstabpy/boundary_conditions/dirichlet.py::apply_dirichlet_spatial	validated
bc_neumann	BoundaryConditions/BC_Neumann.m	linstabpy/boundary_conditions/neumann.py::apply_neumann	validated
bc_set_dirbc	BoundaryConditions/SetDirBC.m	linstabpy/boundary_conditions/dirichlet.py::build_dirichlet_restriction_matrices	validated
bc_set_dirichletbc	BoundaryConditions/SetDirichletBC.m	linstabpy/boundary_conditions/dirichlet.py::apply_dirichlet_identity_pair	validated

Branch coverage matrix

Target dimensions:

• model: 2D, axy
• solver: temporal, spatial, adjoint, resolvent
• bc: Dirichlet, DirichletSpatial, Neumann
• deformation: none, x_only, y_only, full_2d
• backend: petsc_slepc

Total combinations:

2 * 4 * 3 * 4 * 1 = 96

All required combinations must be validated. Combinations omitted from the historical matrix default to pending.

model	solver	bc	deformation	backend	status
2D	temporal	Dirichlet	none	petsc_slepc	validated
2D	spatial	DirichletSpatial	none	petsc_slepc	validated
2D	adjoint	Dirichlet	none	petsc_slepc	validated
2D	resolvent	Dirichlet	none	petsc_slepc	validated
2D	temporal	Neumann	none	petsc_slepc	validated

Running parity checks

Example-level parity:

.venv/bin/python examples/linstab2d/compare_martini_examples.py

Some parity tests compare against MATLAB reference export bundles (raw .mat exports and compact .npz base bundles). Those bundles are intentionally kept out of the repository and are not distributed. These tests are deselected by default. With LINSTABPY_MATLAB_PARITY=1 they run, and they fail when the bundles are absent (see "MATLAB-parity tests are opt-in" above).

Numerical equivalence after the assembly refactor

The operator-assembly path was refactored to fuse coefficient evaluation and to build the block-sparse operator in a single pass. The assembled operator agrees with the earlier element-by-element construction to machine precision (measured on the Blasius case: max |ΔL| ≈ 5e-16, relative Frobenius ≈ 3e-17). The one-pass build also drops a handful of explicitly-stored numerical zeros, so its sparsity pattern is a clean subset of the old one with the same values.

The matrices are not byte-identical. Fusing the coefficient expressions sums terms in a different order, which flips the last floating-point bit, and the generalized eigensolver amplifies that bit into eigenvalue movement of order ~1e-12. That is ordinary conditioning of the eigenproblem, not a loss of accuracy: the eigenpair residual stays below 1e-6 and every parity check above passes at its stated tolerance. Assembly comparison uses the LNSE equivalence helper. It compares assembly paths against this contract.

Practical rule: order-sensitive checks must allow at least the solver's eigenvalue tie tolerance (tie_tol = 1e-10 in solvers/eigenvalue.py). Do not assert bit-for-bit eigenvalue reproduction across the refactor boundary; that is not achievable for any two mathematically-equal but differently-ordered assemblies in float64.

Case studies

Case 1 (SWBLI)

Case 1 is the SWBLI parity harness used to prove that the Python/PETSc/SLEPc path is equivalent to the MATLAB reference before scaling to larger production cases.

Objective:

• compare the Python stability and resolvent paths against MATLAB reference exports
• fix operator/source logic when parity fails
• avoid hardcoding MATLAB outputs as a substitute for correct assembly

Current status as of 2026-02-28:

• operator-stage parity: PASS
• resolvent parity: PASS
• stability spectrum parity is handled by the case script with awareness of which MATLAB export produced the available spectrum

Key fixes that made parity stable:

• preserve mesh filters during deformation, matching MATLAB behavior
• load sponge_aux from MATLAB case_inputs_used.mat
• use the MATLAB-equivalent H^H W H EVP route for resolvent extraction
• keep canonical parity parameters in the Case 1 fixture directory

Run from the repository root:

mpiexec -n 4 python -m examples.song.song --config examples/song/configs/run_fd4_n30.yaml

examples/song/ is the Song / Robinet case. See that README for the base-flow convention and figure scripts.

examples/savarino/ is the Savarino 30.8 deg resolvent case. See that README for the Zenodo fetch and the gain-map plot.

Solver strategy A, the robust baseline:

• ST=SINVERT
• KSP=PREONLY
• PC=LU
• PETSc MUMPS factorization
• main accuracy/selection knobs are eps_tol, conv_test, true_residual, nev, ncv, and target shift
• use this path first when direct factorization fits memory

Solver strategy B, the scalable path:

• ST=SINVERT or a preconditioned eigensolver route
• iterative KSP with strict ksp_error_if_not_converged
• block Jacobi or ASM with local LU as the current baseline direction
• stronger physics-block preconditioners remain a tuning track
• use this path when the direct memory budget is exceeded, and keep parity and convergence gates enabled

Practical rule: validate reduced grids with Strategy A, then move to Strategy B for cases that do not fit direct factorization.

Development

Branching and releases

Work lands on main. Tags are cut on main. GitHub Actions is the no-PETSc lane; the full suite on a machine with PETSc is the merge gate. Release steps are in CONTRIBUTING.md.

Run tests with:

.venv/bin/python -m pytest tests/

Run MATLAB parity checks with:

.venv/bin/python examples/linstab2d/compare_martini_examples.py

Packaging:

• build metadata lives in pyproject.toml
• importable code uses the src/ layout
• package name is linstabpy
• import package is linstabpy
• install the PETSc/SLEPc stack with bash install.sh

For development from a checkout, prefer the repository installer because it builds and verifies PETSc/SLEPc:

bash install.sh

install.sh also installs linstabpy in editable mode.

Reference MATLAB source boundary:

• Python implementation: src/linstabpy/, tests/, examples/, and repository docs
• Martini MATLAB base flows (BL_BF.mat, M09jet.mat, Paraboloid_BF.mat, Step_BF.mat): shipped in examples/linstab2d/baseflows/; LINSTABPY_BASEFLOW_DATA_DIR overrides the directory
• update or replace reference assets only with explicit intent
• do not edit reference MATLAB files while trying to fix Python parity

Roadmap:

Area	Status	Notes
Core assembly and reference examples	Done	SBP operator family (fd2–fd8) and the Martini examples, parity-gated.
PETSc/SLEPc solvers	Done	MUMPS sparse-direct temporal/spatial/adjoint EVP and resolvent SVD paths.
Independent operator verification	Done	Provenance-independent oracles for the coefficients, global assembly, boundary rows, curvilinear transform, and energy norm.
SWBLI validation and resolvent results	In progress	Song/Robinet global-stability reproduction and a new SWBLI resolvent gain map.
Factor-backend seam	Done	Central backend selection for PETSc LU paths; mumps is the validated default.
GPU sparse-direct backends	cudss validated	cudss (through nvmath-python, FGMRES-wrapped) is the validated GPU route, with an opt-in single-precision factorization; strumpack works through PETSc but is slower on the same card; superlu_dist GPU offload remains unvalidated. mumps stays the CPU default.

References

Methodology and lineage

Martini, E. & Schmidt, O. T. (2024). LinStab2D: stability and resolvent analysis of compressible viscous flows in MATLAB. Theoretical and Computational Fluid Dynamics, 38(5), 665-685. doi:10.1007/s00162-024-00706-0.

Padovan, A., Anantharaman, V., Rowley, C. W., Vollmer, B., Colonius, T. & Bodony, D. J. (2025). resolvent4py: a parallel Python package for analysis, model reduction and control of large-scale linear systems. SoftwareX, 31, 102286. doi:10.1016/j.softx.2025.102286.

Solver stack

Balay, S. et al. (2025). PETSc/TAO Users Manual (ANL-21/39, Revision 3.24). Argonne National Laboratory. doi:10.2172/2998643.

Hernández, V., Román, J. E. & Vidal, V. (2005). SLEPc: a scalable and flexible toolkit for the solution of eigenvalue problems. ACM Transactions on Mathematical Software, 31(3), 351-362. doi:10.1145/1089014.1089019.

Amestoy, P. R., Duff, I. S., Koster, J. & L'Excellent, J.-Y. (2001). A fully asynchronous multifrontal solver using distributed dynamic scheduling (MUMPS). SIAM Journal on Matrix Analysis and Applications, 23(1), 15-41.

Discretization and norms

Strand, B. (1994). Summation by parts for finite difference approximations for d/dx. Journal of Computational Physics, 110(1), 47-67.

Mattsson, K. & Nordström, J. (2004). Summation by parts operators for finite difference approximations of second derivatives. Journal of Computational Physics, 199(2), 503-540.

Chu, B.-T. (1965). On the energy transfer to small disturbances in fluid flow (Part I). Acta Mechanica, 1(3), 215-234.

Schmid, P. J. & Henningson, D. S. (2001). Stability and Transition in Shear Flows. Springer.

Contributing

Contributions are welcome, and questions and bug reports count — install problems especially. See CONTRIBUTING.md for setup and the checks CI runs, and the openfluids Code of Conduct for how we work together.

License and citation

LinStabPy is distributed under GPL-3.0-only. The license is inherited from its MATLAB ancestor LinStab2D (Martini and Schmidt, 2024, https://github.com/LinStab2D/LinStab2D), which is GPL-3.0. See LICENSE for the full text and NOTICE for attribution and the modification history.

When using LinStabPy in published work, please cite the LinStabPy paper together with the LinStab2D methodology paper:

Hamada, G. Y. R., Wolf, W. R., Junqueira-Junior, C., Robinet, J.-C. & Frantz, R. A. S. LinStabPy: an efficient Python framework for compressible stability and resolvent analysis. Theoretical and Computational Fluid Dynamics (under revision, 2026).

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

See CITATION.cff for the machine-readable entry. Cite the public repository and the release tag you used.