Metadata-Version: 2.4
Name: linstabpy
Version: 1.0.0
Summary: Stability and resolvent analysis of compressible viscous flows
Author: Gabriel Hamada, William Wolf, Carlos Junqueira-Junior, Jean-Christophe Robinet, Ricardo Frantz
License-Expression: GPL-3.0-only
Project-URL: Homepage, https://linstabpy.vercel.app/
Project-URL: Documentation, https://linstabpy.vercel.app/
Project-URL: Repository, https://github.com/openfluids/linstabpy
Project-URL: Changelog, https://github.com/openfluids/linstabpy/blob/main/CHANGELOG.md
Project-URL: Issues, https://github.com/openfluids/linstabpy/issues
Keywords: stability,resolvent,compressible,Navier-Stokes,PETSc,SLEPc
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
License-File: NOTICE
Requires-Dist: numpy>=1.23
Requires-Dist: scipy>=1.10
Requires-Dist: matplotlib>=3.7
Requires-Dist: pyyaml>=6.0
Provides-Extra: hdf5
Requires-Dist: h5py; extra == "hdf5"
Provides-Extra: petsc
Requires-Dist: mpi4py>=3.1; extra == "petsc"
Requires-Dist: petsc4py>=3.18; extra == "petsc"
Requires-Dist: slepc4py>=3.18; extra == "petsc"
Provides-Extra: cudss
Requires-Dist: nvmath-python[cu12]<2,>=1.0; extra == "cudss"
Dynamic: license-file

![linstabpy banner](docs/assets/readme-banner-v3.jpg)

[![Python](https://img.shields.io/badge/python-3.10%2B-blue.svg)](https://www.python.org/downloads/)
[![PETSc](https://img.shields.io/badge/solvers-PETSc%20%7C%20SLEPc%20%7C%20MUMPS-blue.svg)](https://petsc.org/)
[![License](https://img.shields.io/badge/license-GPL--3.0-blue.svg)](LICENSE)

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

`linstabpy` performs [sparse](https://grokipedia.com/page/Sparse_matrix) linear stability and [resolvent](https://grokipedia.com/page/Resolvent_formalism) analysis of
[compressible](https://grokipedia.com/page/Compressible_flow) viscous flows. It builds the linearized compressible [Navier–Stokes](https://grokipedia.com/page/Navier–Stokes_equations)
operators for two-dimensional structured base flows and solves the [eigenvalue](https://grokipedia.com/page/Eigenvalues_and_eigenvectors) and
resolvent problems through [PETSc](https://en.wikipedia.org/wiki/PETSc)/[SLEPc](https://en.wikipedia.org/wiki/SLEPc) with sparse direct factorization ([MUMPS](https://en.wikipedia.org/wiki/MUMPS_(software))).

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](https://grokipedia.com/page/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](#overview)
- [Installation](#installation)
- [Quick start](#quick-start)
- [Repository layout](#repository-layout)
- [Numerical setup](#numerical-setup)
- [Nondimensionalization](#nondimensionalization)
- [Boundary conditions](#boundary-conditions)
- [Parity with LinStab2D](#parity-with-linstab2d)
- [Case studies](#case-studies)
- [Development](#development)
- [References](#references)
- [License and citation](#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](wiki/core-concepts.md), forcing modes, and response modes
- **Temporal stability:** direct and [adjoint](https://grokipedia.com/page/Adjoint) eigenvalue problems
- **Spatial stability:** streamwise-[wavenumber](https://grokipedia.com/page/Wavenumber) eigenvalue problems
- **Coordinate systems:** Cartesian, cylindrical/axisymmetric, and curvilinear structured meshes
- **Boundary treatment:** [Dirichlet](https://grokipedia.com/page/Dirichlet_boundary_condition) rows, [Neumann](https://grokipedia.com/page/Neumann_boundary_condition) rows, sponge damping fields, and masking for structured grids with excluded points
- **Energy norms:** [Chu](wiki/core-concepts.md) 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](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](wiki/core-concepts.md). The solvers are standard [shift-invert](wiki/solver-glossary.md) [Krylov-Schur](wiki/solver-glossary.md).

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:

```bash
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](https://linstabpy.vercel.app/installation/). When a build
fails, see [Troubleshooting](https://linstabpy.vercel.app/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:

```bash
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

```bash
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](#reference-base-flows).

Run the smallest [Blasius](https://grokipedia.com/page/Blasius_boundary_layer) example:

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

Run a full example:

```bash
.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`:

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

Collect the test suite:

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

Run the suite with:

```bash
.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:

```bash
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](https://portal.nersc.gov/project/sparse/strumpack/), [CUDA](https://en.wikipedia.org/wiki/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](https://grokipedia.com/page/Reynolds_number), [Mach number](https://grokipedia.com/page/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`.

```text
.
|-- 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:

```bash
.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](https://grokipedia.com/page/LU_decomposition) 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](https://grokipedia.com/page/Krylov_subspace) 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:

```text
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](https://en.wikipedia.org/wiki/Condition_number) 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:

```bash
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:

```python
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](https://grokipedia.com/page/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](#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](https://grokipedia.com/page/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](https://grokipedia.com/page/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`:

```python
    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:

```yaml
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

```text
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:

```text
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:

```bash
.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:

```bash
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:

```bash
.venv/bin/python -m pytest tests/
```

Run MATLAB parity checks with:

```bash
.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
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](https://en.wikipedia.org/wiki/Multifrontal_method) 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](https://grokipedia.com/page/Finite_difference_method) 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](http://www.scholarpedia.org/article/Transition_to_turbulence) in Shear
Flows*. Springer.

## Contributing

Contributions are welcome, and questions and bug reports count — install
problems especially. See [CONTRIBUTING.md](CONTRIBUTING.md) for setup and the
checks CI runs, and the
[openfluids Code of Conduct](https://github.com/openfluids/.github/blob/main/CODE_OF_CONDUCT.md)
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.
